Codimension-one bifurcation analysis and chaos of a discrete prey–predator system

Abstract

In this paper, we explore local dynamics at equilibrium points, the existence of bifurcation sets and codimension-one bifurcation analysis of a discrete prey–predator system with Holling type-II functional response. Further, OGY and Hybrid control strategies are utilized to control chaos in the under study discrete model due to the occurrence of Neimark–Sacker and flip bifurcations. Finally, numerical simulations are given to verify theoretical results.

Keywords:

• Predator–prey model
• codimension-one bifurcations
• chaos
• non-standard scheme
• numerical simulation

MSC classifications:

• 70K50
• 92D25
• 40A05

1. Introduction

1.1. Motivation and mathematical formulation

The Lotka–Volterra model, often known as the prey–predator model, is a mathematical representation of the dynamics of two species, one being the predator and the other being the prey in a given population. The model explains the long-term interconnections among the populace of the predator and prey species. It considers variables like the rate of birth and death of the predator and prey, as well as the rate at which the predator hunts and kills the prey. The prey–predator model has been broadly used to explore the dynamics of population interactions in various fields such as designing conservation procedure to safeguard endangered species and understanding the interconnections between several species in an ecosystem. Although these dynamical concerns resulting from the mathematical simulation of prey–predator systems may initially appear to be straightforward, a deep study of these models frequently give rise to incredibly challenging problems. The main goal of modelling a population ecosystem is to ensure that the mathematical model in question can represent the noteworthy system behaviours for the system being observed. Dynamic modelling of ecological systems is often an evolving technique. Prey–predator relationships modelled mathematically using the famous Lotka–Volterra type prey–predator model with Holling type-II functional response [1]: (1) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{ax}(1-\frac{x}{K})-\frac{\mathrm{cxy}}{m+x},{\displaystyle \dot{y}=y(\frac{\mathrm{fx}}{m+x}-d),}}\end{array}\end{array}$(1) where x is prey's density and y is predator's density, while parameters c, a, d, f, K, and m show the capturing rate, intrinsic growth rate of prey, predator death rate, maximal predator growth rate, carrying capacity and half saturation constant, respectively. The so-called ratio-dependent theory, which holds that the rate of per capita predator growth should rely on the ratio of prey to predator abundance, should be the foundation of a more suitable general theory of predator–prey relationships. This is especially true when the predator has to look for food, and it is supported by numerous field, lab studies and observations [2–5]. The basic form of ratio-dependent prey–predator system is: (2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{xf}\left(x\right)-\mathrm{yp}\left(\frac{x}{y}\right),{\displaystyle \dot{y}=(\mathrm{cq}\left(\frac{x}{y}\right)-d)y,}}\end{array}\end{array}$(2) where c is conversion rate, predator functional response is represented by $p(x)$, $q(x)$ is supplanted by $p(x)$. Moreover, Kuang and Beretta [6] have extended the prey–predator model (2(2) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{xf}\left(x\right)-\mathrm{yp}\left(\frac{x}{y}\right),{\displaystyle \dot{y}=(\mathrm{cq}\left(\frac{x}{y}\right)-d)y,}}\end{array}\end{array}$(2) ) by the effect of Holling type-II functional response as follows: (3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{ax}(1-\frac{x}{K})-\frac{\mathrm{cxy}}{\mathrm{my}+x},{\displaystyle \dot{y}=y(-d+\frac{\mathrm{fx}}{\mathrm{my}+x}).}}\end{array}\end{array}$(3) It is anticipated that model (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{ax}(1-\frac{x}{K})-\frac{\mathrm{cxy}}{\mathrm{my}+x},{\displaystyle \dot{y}=y(-d+\frac{\mathrm{fx}}{\mathrm{my}+x}).}}\end{array}\end{array}$(3) ) becomes the form [6, 7]: (4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=x(1-x)-\frac{\mathrm{sxy}}{x+y},{\displaystyle \dot{y}=\mathrm{\delta y}(-r+\frac{x}{x+y}),}}\end{array}\end{array}$(4) by (5) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x=\frac{x}{K},}\\ t=\mathrm{at},\\ {\displaystyle y=\frac{\mathrm{my}}{K},}\end{array}\end{array}$(5) where (6) $\begin{array}{r}\{\begin{array}{l}{\displaystyle s=\frac{c}{\mathrm{ma}},{\displaystyle \delta =\frac{f}{a},{\displaystyle r=\frac{d}{f}.}}}\end{array}\end{array}$(6)

1.2. Literature survey

Many investigators investigated the dynamics of prey–predator models. For instance, Singh and Malik [8] have investigated the dynamics at equilibrium points and exploration the bifurcation of the following model: (7) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=x_t+x_t(1-x_t-\frac{p{y}_t}{\gamma +x_t}-\frac{\alpha }{\delta +x_t}),{\displaystyle y_{t+1}=y_t+\beta y_t(1-\frac{q{y}_t}{\gamma +x_t}),}}\end{array}\end{array}$(7) where $y_t$ and $x_t$ denote the predator and prey populace densities. Ma et al. [9] have studied the local dynamics and bifurcation of following model: (8) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=x_t+\frac{x_t}{1+k{y}_t}(r-a{x}_t)-\frac{(1-m)x_t{y}_t}{b+(1-m)x_t},{\displaystyle y_{t+1}=y_t(1+\mu -\frac{c{y}_t}{b+(1-m)x_t}),}}\end{array}\end{array}$(8) where $r,a,m,b,c,r,\mu ,k$ are positive constants. Santra et al. [10] have explored behaviour of the following model involving prey refuge: (9) $\begin{array}{r}\{\begin{array}{l}{x}_{t+1}=a{x}_t(1-x_t)-c(x_t-b{y}_t)y_t,\\ y_{t+1}=d(x_t-b{y}_t)y_t,\end{array}\end{array}$(9) where a, c, b, r, d are positive parameters. Chen et al. [11] have studied local stability and bifurcation at fixed points of the following model: (10) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=x_t+x_t(1-x_t)\frac{x_t}{u+x_t}-b{x}_t{y}_t,}\\ y_{t+1}=y_t+c{x}_t{y}_t+m{y}_t-y_t^2,\end{array}\end{array}$(10) where u, b, c, m are positive constants. In contrast to discrete models, many investigators investigated the dynamics of continuous-time systems designated by differential equations. For instance, Tunç and Tunç [12] have explored dynamics of the following system of second-order differential equations with delay: (11) $\begin{array}{r}\ddot{x}+F(x,\dot{x})\dot{x}+H(x(t-\tau (t)))=P(t,x,\dot{x}).\end{array}$(11) Tunç [13] has investigated the dynamics of the following vector lienard equation: (12) $\begin{array}{r}\ddot{x}+F(x,\dot{x})\dot{x}+H(x(t-\tau ))=P\left(t\right).\end{array}$(12) Tunç [14] has investigated the boundedness of solutions of a class of non-autonomous differential equations.

1.3. The novelty of the proposed work

It is important here to note that the earlier work is about the global qualitative study on system (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{ax}(1-\frac{x}{K})-\frac{\mathrm{cxy}}{\mathrm{my}+x},{\displaystyle \dot{y}=y(-d+\frac{\mathrm{fx}}{\mathrm{my}+x}).}}\end{array}\end{array}$(3) ) where the discussion of authors focuses on the qualitative study of the ratio-dependent predator–prey systems. Their work builds upon previous work but aims to address the open questions regarding the global qualitative behaviour of the model. By transforming the Michaelis-Menten-type ratio-dependent model into a Gause-type predator–prey system, the authors obtain a complete classification of the asymptotic behaviour of the solutions. They resolve open questions regarding the global stability of equilibria and the uniqueness of limit cycles. Their work provides insights into how the outcomes of the model depend on initial conditions. Furthermore, they present the biological implications of the findings, showcasing the practical relevance of the study. Overall, Kuang and Beretta [6] discussed the dynamics of continuous-time model (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{ax}(1-\frac{x}{K})-\frac{\mathrm{cxy}}{\mathrm{my}+x},{\displaystyle \dot{y}=y(-d+\frac{\mathrm{fx}}{\mathrm{my}+x}).}}\end{array}\end{array}$(3) ), and demonstrate that if the system's positive steady state (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{ax}(1-\frac{x}{K})-\frac{\mathrm{cxy}}{\mathrm{my}+x},{\displaystyle \dot{y}=y(-d+\frac{\mathrm{fx}}{\mathrm{my}+x}).}}\end{array}\end{array}$(3) ) is locally asymptotically stable, then the system cannot have non-trivial positive periodic solutions. It also includes some results regarding the global stability of the positive steady state. Moreover, they non-dimensionalize the system (3(3) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=\mathrm{ax}(1-\frac{x}{K})-\frac{\mathrm{cxy}}{\mathrm{my}+x},{\displaystyle \dot{y}=y(-d+\frac{\mathrm{fx}}{\mathrm{my}+x}).}}\end{array}\end{array}$(3) ) into (4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=x(1-x)-\frac{\mathrm{sxy}}{x+y},{\displaystyle \dot{y}=\mathrm{\delta y}(-r+\frac{x}{x+y}),}}\end{array}\end{array}$(4) ), then, Hsu et al. [7] investigated the full classification for the asymptotic behaviour of (4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=x(1-x)-\frac{\mathrm{sxy}}{x+y},{\displaystyle \dot{y}=\mathrm{\delta y}(-r+\frac{x}{x+y}),}}\end{array}\end{array}$(4) ) by transforming it into a Gause-type predator–prey system.

In contrast to the continuous-time model, our aim in this paper is to explore the dynamical characteristics of the following prey–predator system: (13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) which is discrete version of (4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=x(1-x)-\frac{\mathrm{sxy}}{x+y},{\displaystyle \dot{y}=\mathrm{\delta y}(-r+\frac{x}{x+y}),}}\end{array}\end{array}$(4) ) by non-standard finite difference scheme, where h is a step size. Furthermore, s denotes the ratio of capturing rate to the half saturation constant times intrinsic growth rate of prey, δ is the ratio of maximal predator growth rate to the intrinsic growth rate of prey, and finally r is the ratio of predator death rate to the maximal predator growth rate. The reason for studying the dynamics of discrete-time system (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ), instead of of continuous-time system (4(4) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \dot{x}=x(1-x)-\frac{\mathrm{sxy}}{x+y},{\displaystyle \dot{y}=\mathrm{\delta y}(-r+\frac{x}{x+y}),}}\end{array}\end{array}$(4) ), is that the discrete model simplifies the complexity of the system, making it more accessible for analysis and interpretation. Continuous models often involve differential equations, which can be challenging to solve analytically or numerically. By converting the model into a set of difference equations, the computational complexity is reduced, making it easier to simulate and analyse the dynamics of the system. Additionally, the discretization of the model allows for direct comparison with experimental data. Experimental data is often collected at discrete-time points, and by converting the model into a discrete version, the model's predictions can be directly compared with the observed data. This facilitates experimental validation and enhances the credibility of the model. Furthermore, stability analysis becomes more feasible with discrete models. By studying the behaviour of the system over discrete-time steps, it becomes possible to analyse the stability of the fixed points and determine if the system will converge to an equilibrium or exhibit oscillatory behaviour. This stability analysis can provide valuable insights for system design and control. Lastly, parameter estimation is often easier with discrete models. Estimating parameters in continuous models can be challenging due to the complexity of the equations. Discrete models, with their simpler forms, allow for more straightforward parameter estimation techniques, enhancing the accuracy and reliability of the model's predictions. In summary, converting a continuous-time model into a discrete version offers advantages in terms of computational efficiency, simplicity, compatibility with experimental data, stability analysis, and parameter estimation. These advantages make the discrete model a valuable tool in the literature of ecology and biology for studying and analysing the local behaviour of fixed points. In the context of ecology and biology, calculating the local behaviour for fixed points in discrete models offers several advantages. Firstly, it allows for the identification of bifurcation sets, which are critical points where the qualitative behaviour of the system changes. By analysing the local behaviour of fixed points, researchers can identify bifurcations such as the emergence of multiple stable states or the occurrence of oscillatory behaviour. Understanding these bifurcation sets provides insights into the underlying mechanisms driving system dynamics. Secondly, discrete-time models enable bifurcation analysis, which involves studying how the system's behaviour changes as model parameters are varied. This analysis can reveal the existence of different types of bifurcations, such as saddle-node, transcritical, Neimark–Sacker, flip, or pitchfork bifurcations. Bifurcation analysis helps in understanding the stability and dynamics of the system under different conditions. Additionally, discrete models are useful for studying chaos in ecological and biological systems. Techniques like the OGY (Ott-Grebogi-Yorke) method can be employed to analyse the local behaviour of fixed points and identify chaotic regions in the parameter space. This study of chaos can provide insights into the underlying mechanisms driving chaotic behaviour. Furthermore, discrete-time models facilitate the development and analysis of hybrid control strategies. Hybrid control combines continuous and discrete dynamics to design control strategies that effectively manage complex systems. By calculating the local behaviour for fixed points, researchers can determine control strategies that stabilize the system, suppress undesirable behaviour, or guide the system towards desired states. Lastly, discrete-time models offer opportunities for numerical validation of theoretical results. Researchers can simulate the discrete model and compare the results with theoretical predictions, validating the accuracy and reliability of their theoretical findings. This numerical validation helps build confidence in the model and its ability to capture the essential dynamics of the system. In summary, calculating the local behaviour for fixed points in discrete models provides advantages in identifying bifurcation sets, conducting bifurcation analysis, studying chaos, developing hybrid control strategies, and numerically validating theoretical results. These advantages enhance our understanding of ecological and biological systems and enable more effective decision-making and control strategies. Due to aforementioned fact, in this paper we explore dynamics of discrete system model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) where our key investigations include:

• Local behaviour for fixed points.

• Identification of bifurcation sets for fixed points.

• Bifurcation analysis of discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ).

• Study of chaos by OGY and Hybrid control strategies.

• Numerical validation of theoretical results.

1.4. The advantages of the proposed work

In the setting of a discrete prey–predator system, codimension-one bifurcation analysis has numerous advantages for understanding the dynamics, stability, and potential for chaos in ecological systems:

• Codimension-one bifurcation analysis aids in identifying key parameter value at which qualitative changes in the system dynamics occur. It enables researchers to identify the precise parameter value that cause bifurcations or the formation of chaotic attractors.

• This type of analysis gives information about the stability of the system's equilibrium points and periodic orbits. It helps identify whether prey and predator populations tend to stabilize or oscillate under varied parametric value.

• Researchers can forecast system transitions by examining codimension-one bifurcations. They can, for example, forecast when a stable equilibrium will become unstable, resulting in the formation of chaotic behaviour.

• Codimension-one bifurcation analysis can aid in detecting the emergence of chaos in a discrete prey–predator system.

• Codimension-one bifurcation analysis can be used as a teaching tool for students and researchers to learn about complex ecological dynamics.

• The calculation of local behaviour for fixed points in discrete-time models offers advantages in identifying bifurcation sets, conducting bifurcation analysis, studying chaos, developing hybrid control strategies, and numerically validating theoretical results. These advantages enhance our understanding of ecological systems, aid in predicting and managing ecological dynamics, and inform conservation and ecosystem management practices.

1.5. Paper layout

The organization of rest of the paper is as follows: stability analysis and bifurcation sets are studied in Section 2. In Section 3, we studied flip bifurcation for boundary fixed point, whereas Neimark–Sacker and flip bifurcations at interior fixed points are briefly investigated in Sections 4 and 5, respectively. The chaos control by OGY and Hybrid control strategies is presented in Section 6. Finally, numerical simulations and conclusion are given in Sections 7 and 8, respectively.

2. Bifurcation sets and stability analysis

In the present section, we explore the stability analysis at equilibrium points and bifurcation sets for the discrete prey–predator model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ). If $\phi =(x,y)$ is an equilibrium point of (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ), then one has (14) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x=\frac{x(1+h)}{1+\mathrm{hx}}-\frac{\mathrm{hsxy}}{(x+y)(1+\mathrm{hx})},{\displaystyle y=y+\mathrm{h\delta y}(-r+\frac{x}{x+y}).}}\end{array}\end{array}$(14) Since ${\phi }_1=(1,0)$ satisfied system (14(14) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x=\frac{x(1+h)}{1+\mathrm{hx}}-\frac{\mathrm{hsxy}}{(x+y)(1+\mathrm{hx})},{\displaystyle y=y+\mathrm{h\delta y}(-r+\frac{x}{x+y}).}}\end{array}\end{array}$(14) ) obviously, and therefore for all h, s, δ and r discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) has semitrivial equilibrium point ${\phi }_1=(1,0)$. In order to get the model's interior equilibrium point, one need to solve the following system (15) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{1+h}{1+\mathrm{hx}}-\frac{\mathrm{hsy}}{(x+y)(1+\mathrm{hx})}=1,{\displaystyle \delta (-r+\frac{x}{x+y})=0.}}\end{array}\end{array}$(15) Equation (15(15) $\begin{array}{r}\{\begin{array}{l}{\displaystyle \frac{1+h}{1+\mathrm{hx}}-\frac{\mathrm{hsy}}{(x+y)(1+\mathrm{hx})}=1,{\displaystyle \delta (-r+\frac{x}{x+y})=0.}}\end{array}\end{array}$(15) ) can also be rewritten as (16) $\begin{array}{r}y=\frac{x^2-x}{1-x-s},\end{array}$(16) and (17) $\begin{array}{r}y=\frac{(1-r)x}{r}.\end{array}$(17) Using (17(17) $\begin{array}{r}y=\frac{(1-r)x}{r}.\end{array}$(17) ) into (16(16) $\begin{array}{r}y=\frac{x^2-x}{1-x-s},\end{array}$(16) ), one gets (18) $\begin{array}{r}x=1-s+\mathrm{sr}.\end{array}$(18) Further, using (18(18) $\begin{array}{r}x=1-s+\mathrm{sr}.\end{array}$(18) ) into (17(17) $\begin{array}{r}y=\frac{(1-r)x}{r}.\end{array}$(17) ), one gets (19) $\begin{array}{r}y=\frac{(1-r)(1-s+\mathrm{sr})}{r}.\end{array}$(19) So, from Equations (18(18) $\begin{array}{r}x=1-s+\mathrm{sr}.\end{array}$(18) ) and (19(19) $\begin{array}{r}y=\frac{(1-r)(1-s+\mathrm{sr})}{r}.\end{array}$(19) ), one can obtained that if $\frac{s-1}{s}<r<1$ then in ${\mathbb{R}}_{+}^2$ model's interior equilibrium point is ${\phi }_2=(1-s+\mathrm{sr},\frac{(1-r)(1-s+\mathrm{sr})}{r})$. Now the variation matrix $V{|}_{\phi }$ of the linearized system of (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) under the map $(f_1,f_2)\mapsto (x_{t+1},y_{t+1})$ is (20) $\begin{array}{r}V{|}_{\phi }:=\left(\begin{array}{ccc}{\displaystyle \frac{\begin{array}{c}(1+h)(x+y{)}^2\\ +h^{2}s{x}^{2}y-hs{y}^2\end{array}}{(x+y{)}^2(1+\mathrm{hx}{)}^2}}& {\displaystyle \frac{-\mathrm{hs}{x}^2}{(x+y{)}^2(1+\mathrm{hx})}{\displaystyle \frac{\mathrm{h\delta }{y}^2}{(x+y{)}^2}}}& {\displaystyle 1-\mathrm{h\delta r}+\frac{\mathrm{h\delta }{x}^2}{(x+y{)}^2}}\end{array}\right),\end{array}$(20) where (21) $\begin{array}{r}\{\begin{array}{l}{\displaystyle f_1:=\frac{x(1+h)}{1+\mathrm{hx}}-\frac{\mathrm{hsxy}}{(x+y)(1+\mathrm{hx})},{\displaystyle f_2:=y+\mathrm{h\delta y}(-r+\frac{x}{x+y}).}}\end{array}\end{array}$(21) Hereafter, we will give local dynamics for the equilibria ${\phi }_{1,2}$ by stability theory [15–17]. So for ${\phi }_1$, (20(20) $\begin{array}{r}V{|}_{\phi }:=\left(\begin{array}{ccc}{\displaystyle \frac{\begin{array}{c}(1+h)(x+y{)}^2\\ +h^{2}s{x}^{2}y-hs{y}^2\end{array}}{(x+y{)}^2(1+\mathrm{hx}{)}^2}}& {\displaystyle \frac{-\mathrm{hs}{x}^2}{(x+y{)}^2(1+\mathrm{hx})}{\displaystyle \frac{\mathrm{h\delta }{y}^2}{(x+y{)}^2}}}& {\displaystyle 1-\mathrm{h\delta r}+\frac{\mathrm{h\delta }{x}^2}{(x+y{)}^2}}\end{array}\right),\end{array}$(20) ) becomes (22) $\begin{array}{r}V{|}_{{\phi }_1}:=\left(\begin{array}{ccc}{\displaystyle \frac{1}{1+h}}& {\displaystyle \frac{-\mathrm{hs}}{1+h}0}& 1-\mathrm{h\delta r}+\mathrm{h\delta }\end{array}\right),\end{array}$(22) with (23) $\begin{array}{r}{\lambda }_1=\frac{1}{1+h},{\lambda }_2=1-\mathrm{h\delta r}+\mathrm{h\delta }.\end{array}$(23)

Theorem 2.1

${\phi }_1$ of model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) is

1. a sink if $h>\frac{2}{\delta (r-1)}$;

2. never source;

3. a saddle if $0<h<\frac{2}{\delta (r-1)}$;

4. non-hyperbolic if $h=\frac{2}{\delta (r-1)}$.

Proof.

By stability theory, ${\phi }_1$ of discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) is a sink if $|{\lambda }_{1,2}|<1$. Therefore, from (23(23) $\begin{array}{r}{\lambda }_1=\frac{1}{1+h},{\lambda }_2=1-\mathrm{h\delta r}+\mathrm{h\delta }.\end{array}$(23) ) if $|{\lambda }_1|=\frac{1}{1+h}<1$ and $|{\lambda }_2|=|1-\mathrm{h\delta r}+\mathrm{h\delta }|<1$, that is, $h>\frac{2}{\delta (r-1)}$ then ${\phi }_1$ is a sink. Similarly, it is easy to obtain that ${\phi }_1$ is never source, saddle if $0<h<\frac{2}{\delta (r-1)}$, and non-hyperbolic if $h=\frac{2}{\delta (r-1)}$.

Now in the following, one has the following result regarding bifurcation set at ${\phi }_1$ if $h=\frac{2}{\delta (r-1)}$ holds.

Theorem 2.2

For the model's semitrivial fixed point ${\phi }_1$, the flip bifurcation set is written as $\mathcal{F}{|}_{{\phi }_1}:=\{(h,s,\delta ,r),h=\frac{2}{\delta (r-1)}\}$.

Proof.

Since ${\phi }_1$ is non-hyperbolic if $h=\frac{2}{\delta (r-1)}$. Therefore, if $h=\frac{2}{\delta (r-1)}$ holds then ${\lambda }_2{|}_{h=\frac{2}{\delta (r-1)}}=-1$ but ${\lambda }_1{|}_{h=\frac{2}{\delta (r-1)}}=\frac{\delta (r-1)}{\delta (r-1)+2}\ne 1\mathrm{or}-1$ which implies that at ${\phi }_1$ eigenvalues criterion for the existence of flip bifurcation holds if $(h,s,\delta ,r)$ passes $\mathcal{F}{|}_{{\phi }_1}:=\{(h,s,\delta ,r),h=\frac{2}{\delta (r-1)}\}$.

Now for ${\phi }_2$, (20(20) $\begin{array}{r}V{|}_{\phi }:=\left(\begin{array}{ccc}{\displaystyle \frac{\begin{array}{c}(1+h)(x+y{)}^2\\ +h^{2}s{x}^{2}y-hs{y}^2\end{array}}{(x+y{)}^2(1+\mathrm{hx}{)}^2}}& {\displaystyle \frac{-\mathrm{hs}{x}^2}{(x+y{)}^2(1+\mathrm{hx})}{\displaystyle \frac{\mathrm{h\delta }{y}^2}{(x+y{)}^2}}}& {\displaystyle 1-\mathrm{h\delta r}+\frac{\mathrm{h\delta }{x}^2}{(x+y{)}^2}}\end{array}\right),\end{array}$(20) ) becomes (24) $\begin{array}{r}V{|}_{{\phi }_2}:=\left(\begin{array}{ccc}{\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}& {\displaystyle \frac{-h{r}^{2}s}{1+h+\mathrm{hs}(r-1)}\mathrm{h\delta }(r-1{)}^2}& 1+\mathrm{hr\delta }(r-1)\end{array}\right),\end{array}$(24) with the following characteristic equation:

(25) $\begin{array}{r}det\left(\begin{array}{ccc}{\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}-\lambda }& {\displaystyle \frac{-h{r}^{2}s}{1+h+\mathrm{hs}(r-1)}\mathrm{h\delta }(r-1{)}^2}& 1+\mathrm{hr\delta }(r-1)-\lambda \end{array}\right)=0,\end{array}$(25)

that is, (26) $\begin{array}{r}{\lambda }^2-{\mathrm{\Lambda }}_1\lambda +{\mathrm{\Lambda }}_2=0,\end{array}$(26) where (27) $\begin{array}{rl}{\mathrm{\Lambda }}_1& =\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr\delta }(r-1),\\ {\mathrm{\Lambda }}_2& =\left(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}\right)(1+\mathrm{hr\delta }(r-1))\\ & +\frac{h^2{r}^2\mathrm{s\delta }(r-1{)}^2}{1+h+\mathrm{hs}(r-1)}.\end{array}$(27) From (26(26) $\begin{array}{r}{\lambda }^2-{\mathrm{\Lambda }}_1\lambda +{\mathrm{\Lambda }}_2=0,\end{array}$(26) ), one has (28) $\begin{array}{r}{\lambda }_{1,2}=\frac{{\mathrm{\Lambda }}_1\pm \sqrt{\mathrm{\Delta }}}{2},\end{array}$(28) where (29) $\begin{array}{rl}\mathrm{\Delta }& ={\mathrm{\Lambda }}_1^2-4{\mathrm{\Lambda }}_2,\\ & ={(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr\delta }(r-1))}^2\\ & -4(\left(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}\right)(1+\mathrm{hr\delta }(r-1))\\ & +\frac{h^2{r}^2\mathrm{s\delta }(r-1{)}^2}{1+h+\mathrm{hs}(r-1)}).\end{array}$(29)

Theorem 2.3

If $\mathrm{\Delta }<0$, then the model's interior equilibrium point ${\phi }_2$ is

1. a stable focus if $0<\delta <\frac{1-s+s{r}^2}{r^2-r}$ with $r>max\{\pm \sqrt{\frac{s-1}{s}},1\}$;

2. an unstable focus if $\delta >\frac{1-s+s{r}^2}{r^2-r}$ with $r>max\{\pm \sqrt{\frac{s-1}{s}},1\}$;

3. non-hyperbolic if $\delta =\frac{1-s+s{r}^2}{r^2-r}$.

Proof.

If $\mathrm{\Delta }<0$ then from (28(28) $\begin{array}{r}{\lambda }_{1,2}=\frac{{\mathrm{\Lambda }}_1\pm \sqrt{\mathrm{\Delta }}}{2},\end{array}$(28) ) one gets $|{\lambda }_{1,2}|=\sqrt{(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)})(1+\mathrm{hr\delta }(r-1))+\frac{h^2{r}^2\mathrm{s\delta }(r-1{)}^2}{1+h+\mathrm{hs}(r-1)}}<1$. This implies that if $0<\delta <\frac{1-s+s{r}^2}{r^2-r}$ with $r>max\{\pm \sqrt{\frac{s-1}{s}},1\}$ then ${\phi }_2$ is a stable focus. Similarly, easy calculation yields the fact that it an unstable focus if $\delta >\frac{1-s+s{r}^2}{r^2-r}$ with $r>max\{\pm \sqrt{\frac{s-1}{s}},1\}$, and non-hyperbolic if $\delta =\frac{1-s+s{r}^2}{r^2-r}$.

Theorem 2.4

If $\mathrm{\Delta }\ge 0$ then the model's interior equilibrium point ${\phi }_2$ is

1. a stable node if $\delta >\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$ with $s>max\{\frac{2+h}{h(r-1{)}^2},\frac{2+h}{h(1-r)}\}$;

2. an unstable node if $0<\delta <\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$ and $s>max\{\frac{2+h}{h(r-1{)}^2},\frac{2+h}{h(1-r)}\}$;

3. non-hyperbolic if $\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$.

Proof.

It is noted that ${\phi }_2$ is a stable node if real roots of (28(28) $\begin{array}{r}{\lambda }_{1,2}=\frac{{\mathrm{\Lambda }}_1\pm \sqrt{\mathrm{\Delta }}}{2},\end{array}$(28) ) satisfying $|{\lambda }_{1,2}|<1$, that is, if $\delta >\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$ with $s>max\{\frac{2+h}{h(r-1{)}^2},\frac{2+h}{h(1-r)}\}$ then ${\phi }_2$ is a stable node. Similarly, easy calculation also implies that it is an unstable node if $0<\delta <\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$, and non-hyperbolic if $\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$.

Now in the following, one has the following result regarding bifurcation set at ${\phi }_2$ if $\delta =\frac{1-s+s{r}^2}{r^2-r}$ and $\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$ hold, respectively.

Theorem 2.5

For the model's interior equilibrium point ${\phi }_2$, the Neimark–Sacker and flip bifurcation sets, respectively are

1. $\mathcal{N}{|}_{{\phi }_2}:=\{(h,s,\delta ,r),\delta =\frac{1-s+s{r}^2}{r^2-r}\}$;

2. $\mathcal{F}{|}_{{\phi }_2}:=\{(h,s,\delta ,r),\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}\}$.

Proof.

(i). It is noted that if $\delta =\frac{1-s+s{r}^2}{r^2-r}$ holds, then $V{|}_{{\phi }_2}$ has complex eigenvalues with $|{\lambda }_{1,2}{|}_{\delta =\frac{1-s+s{r}^2}{r^2-r}}=1$ which implies that at ${\phi }_2$ eigenvalues criterion for the existence of N-S bifurcation holds if $(h,s,\delta ,r)$ passes $\mathcal{N}{|}_{{\phi }_2}:=\{(h,s,\delta ,r),\delta =\frac{1-s+s{r}^2}{r^2-r}\}$.

(ii). On the other hand, if $\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$ holds then one has ${\lambda }_1{|}_{\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}}=-1$ but ${\lambda }_2{|}_{\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}}=\frac{2+h(1+s(r-1))(1+\mathrm{hrs}(r-1))}{(1+h+\mathrm{hs}(r-1))(2+h+\mathrm{hs}(r-1))}\ne 1\mathrm{or}-1$ which implies that at ${\phi }_2$ eigenvalues criterion for the existence of flip bifurcation holds if $(h,s,\delta ,r)$ passes through the curve $\mathcal{F}{|}_{{\phi }_2}:=\{(h,s,\delta ,r),\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}\}$.

3. Analysis of a flip bifurcation for semitrivial fixed point

Based on the local dynamic study for semitrivial fixed point ${\phi }_1$ in Section 2, here we investigate flip bifurcation by bifurcation theory [18–27].

Theorem 3.1

If $(h,s,\delta ,r)\in \mathcal{F}{|}_{{\phi }_1}$ then for ${\phi }_1$ the flip bifurcation does not take place in the discrete-time model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ).

Proof.

Since model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) is invariant under y = 0, and so one restrict it on the line y = 0 where it becomes (30) $\begin{array}{r}{x}_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}.\end{array}$(30) From (30(30) $\begin{array}{r}{x}_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}.\end{array}$(30) ), one has (31) $\begin{array}{r}f(h,x):=\frac{x(1+h)}{1+\mathrm{hx}}.\end{array}$(31) Finally, if $x=x^{\ast }=1$ and $h=h^{\ast }=\frac{2}{\delta (r-1)}$, from (31(31) $\begin{array}{r}f(h,x):=\frac{x(1+h)}{1+\mathrm{hx}}.\end{array}$(31) ), one has (32) $\begin{array}{r}{\frac{\mathrm{\partial }f}{\mathrm{\partial }x}|}_{h^{\ast }=\frac{2}{\delta (r-1)},x^{\ast }=1}:=\frac{\delta (r-1)}{\delta (r-1)+2}\ne -1,\end{array}$(32) (33) $\begin{array}{r}{\frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{x}^2}|}_{h^{\ast }\frac{2}{\delta (r-1)},x^{\ast }=1}:=\frac{-4\delta (r-1)}{{(\delta (r-1)+2)}^2}\ne 0,\end{array}$(33) and (34) $\begin{array}{r}{\frac{\mathrm{\partial }f}{\mathrm{\partial }h}|}_{h^{\ast }=\frac{2}{\delta (r-1)},x^{\ast }=1}:=0.\end{array}$(34) From (32(32) $\begin{array}{r}{\frac{\mathrm{\partial }f}{\mathrm{\partial }x}|}_{h^{\ast }=\frac{2}{\delta (r-1)},x^{\ast }=1}:=\frac{\delta (r-1)}{\delta (r-1)+2}\ne -1,\end{array}$(32) ) and (34(34) $\begin{array}{r}{\frac{\mathrm{\partial }f}{\mathrm{\partial }h}|}_{h^{\ast }=\frac{2}{\delta (r-1)},x^{\ast }=1}:=0.\end{array}$(34) ), one can be concluded that the flip bifurcation does not take place if $(h,s,\delta ,r)\in \mathcal{F}{|}_{{\phi }_1}$.

4. Analysis of a Neimark–Sacker bifurcation for interior fixed point

Based on the local dynamic study for the equilibrium ${\phi }_2$ in Section 2, we investigate the Neimark–Sacker bifurcation if $(h,s,\delta ,r)\in \mathcal{N}{|}_{{\phi }_2}$ where δ is regarded as a bifurcation parameter.

Theorem 4.1

At ${\phi }_2$, discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) undergoes N–S bifurcation if $(h,s,\delta ,r)\in \mathcal{N}{|}_{{\phi }_2}$.

Proof.

Recall that if δ is a bifurcation parameter then perturbation system of discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) takes the form (35) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+h({\delta }^{\ast }+ϵ)y_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(35) where $\delta ={\delta }^{\ast }+ϵ$ and $0<|{\delta }^{\ast }|\ll 1$. For the perturbation system (35(35) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+h({\delta }^{\ast }+ϵ)y_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(35) ), the roots of $V{|}_{{\phi }_2}$ are (36) $\begin{array}{r}{\lambda }_{1,2}:=\frac{{\mathrm{\Lambda }}_1(ϵ)\pm \iota \sqrt{4{\mathrm{\Lambda }}_2(ϵ)-{\mathrm{\Lambda }}_1^2(ϵ)}}{2},\end{array}$(36) where ${\mathrm{\Lambda }}_1(ϵ)=\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr}({\delta }^{\ast }+ϵ)(r-1)$ and ${\mathrm{\Lambda }}_2(ϵ)=(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)})(1+\mathrm{hr}({\delta }^{\ast }+ϵ)(r-1))+\frac{h^2{r}^{2}s({\delta }^{\ast }+ϵ)(r-1{)}^2}{1+h+\mathrm{hs}(r-1)}$. Now due to the fact that $(h,s,\delta ,r)\in \mathcal{N}{|}_{{\phi }_2}$, the characteristic equation (36(36) $\begin{array}{r}{\lambda }_{1,2}:=\frac{{\mathrm{\Lambda }}_1(ϵ)\pm \iota \sqrt{4{\mathrm{\Lambda }}_2(ϵ)-{\mathrm{\Lambda }}_1^2(ϵ)}}{2},\end{array}$(36) ) of $V{|}_{{\phi }_2}$ has two conjugate complex root with $|{\lambda }_{1,2}|=\sqrt{\begin{array}{c}\left(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}\right)(1+\mathrm{hr}({\delta }^{\ast }+ϵ)(r-1))\\ +\frac{h^2{r}^{2}s({\delta }^{\ast }+ϵ)(r-1{)}^2}{1+h+\mathrm{hs}(r-1)}\end{array}}$, and subsequently one gets the following quantity by noticing that $\delta =\frac{1-s+s{r}^2}{r^2-r}$: (37) $\begin{array}{r}\frac{d|{\lambda }_{1,2}|}{\mathrm{dϵ}}{|}_{ϵ=0}=\frac{\mathrm{hr}(r-1)}{2(1+h+\mathrm{hs}(r-1))}\ne 0.\end{array}$(37) Moreover for $ϵ=0$, the characteristics roots of (36(36) $\begin{array}{r}{\lambda }_{1,2}:=\frac{{\mathrm{\Lambda }}_1(ϵ)\pm \iota \sqrt{4{\mathrm{\Lambda }}_2(ϵ)-{\mathrm{\Lambda }}_1^2(ϵ)}}{2},\end{array}$(36) ) must not be found at the intersections of unit circle with coordinate axes. Therefore, one need to suppose that ${\lambda }_1^{\ell }\ne 1$ and ${\lambda }_2^{\ell }\ne 1(\ell =1,2,3,4)$ for $ϵ=0$, which is equivalent to ${\mathrm{\Lambda }}_1(0)=\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr\delta }(r-1)\ne -2,0,1,2$. But if $\delta =\frac{1-s+s{r}^2}{r^2-r}$ holds, then ${\mathrm{\Lambda }}_2(0)=1$, and so, ${\mathrm{\Lambda }}_1(0)\ne -2,2$. Therefore, ${\mathrm{\Lambda }}_1(0)=\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr\delta }(r-1)\ne 0,1$, and so, by straightforward calculation one obtains (38) $\begin{array}{rl}s& \ne \frac{-2-h(2+r)\pm \sqrt{-4-8r+\mathrm{hr}(-4+\mathrm{hr})}}{2h(r^2-1)},\\ & \frac{-1-h(2+r)\pm \sqrt{-3-4r+\mathrm{hr}(-2+\mathrm{hr})}}{2h(r^2-1)}.\end{array}$(38) Now to transform ${\phi }_2$ one use the following transformations: (39) $\begin{array}{r}{u}_t=x_t-x^{\ast },v_t=y_t-y^{\ast },\end{array}$(39) with $x^{\ast }=1-s+\mathrm{sr}$ and $y^{\ast }=\frac{(1-r)(1-s+\mathrm{sr})}{r}$. From (39(39) $\begin{array}{r}{u}_t=x_t-x^{\ast },v_t=y_t-y^{\ast },\end{array}$(39) ) and (35(35) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+h({\delta }^{\ast }+ϵ)y_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(35) ), we get (40) $\begin{array}{r}\{\begin{array}{l}{\displaystyle u_{t+1}={\displaystyle \frac{(u_t+x^{\ast })(1+h)}{1+h(u_t+x^{\ast })}}-{\displaystyle \frac{\mathrm{hs}(u_t+x^{\ast })(v_t+y^{\ast })}{((u_t+x^{\ast })+(v_t+y^{\ast }))(1+h(u_t+x^{\ast }))}}-x^{\ast },{\displaystyle v_{t+1}=(v_t+y^{\ast })+h({\delta }^{\ast }+ϵ)(v_t+y^{\ast })(-r+{\displaystyle \frac{(u_t+x^{\ast })}{(u_t+x^{\ast })+(v_t+y^{\ast })}})-y^{\ast }.}}\end{array}\end{array}$(40) If $ϵ=0$, then we will derive normal form of (40(40) $\begin{array}{r}\{\begin{array}{l}{\displaystyle u_{t+1}={\displaystyle \frac{(u_t+x^{\ast })(1+h)}{1+h(u_t+x^{\ast })}}-{\displaystyle \frac{\mathrm{hs}(u_t+x^{\ast })(v_t+y^{\ast })}{((u_t+x^{\ast })+(v_t+y^{\ast }))(1+h(u_t+x^{\ast }))}}-x^{\ast },{\displaystyle v_{t+1}=(v_t+y^{\ast })+h({\delta }^{\ast }+ϵ)(v_t+y^{\ast })(-r+{\displaystyle \frac{(u_t+x^{\ast })}{(u_t+x^{\ast })+(v_t+y^{\ast })}})-y^{\ast }.}}\end{array}\end{array}$(40) ). On expanding (40(40) $\begin{array}{r}\{\begin{array}{l}{\displaystyle u_{t+1}={\displaystyle \frac{(u_t+x^{\ast })(1+h)}{1+h(u_t+x^{\ast })}}-{\displaystyle \frac{\mathrm{hs}(u_t+x^{\ast })(v_t+y^{\ast })}{((u_t+x^{\ast })+(v_t+y^{\ast }))(1+h(u_t+x^{\ast }))}}-x^{\ast },{\displaystyle v_{t+1}=(v_t+y^{\ast })+h({\delta }^{\ast }+ϵ)(v_t+y^{\ast })(-r+{\displaystyle \frac{(u_t+x^{\ast })}{(u_t+x^{\ast })+(v_t+y^{\ast })}})-y^{\ast }.}}\end{array}\end{array}$(40) ) by Taylor series up to second-order at origin, one gets (41) $\begin{array}{r}\{\begin{array}{l}{u}_{t+1}={\mathrm{\Omega }}_{11}^1{u}_t+{\mathrm{\Omega }}_{12}^1{v}_t+{\mathrm{\Omega }}_{13}^1{u}_t^2+{\mathrm{\Omega }}_{14}^1{u}_t{v}_t+{\mathrm{\Omega }}_{15}^1{v}_t^2,\\ v_{t+1}={\mathrm{\Omega }}_{11}^2{u}_t+{\mathrm{\Omega }}_{12}^2{v}_t+{\mathrm{\Omega }}_{13}^2{u}_t^2+{\mathrm{\Omega }}_{14}^2{u}_t{v}_t+{\mathrm{\Omega }}_{15}^2{v}_t^2,\end{array}\end{array}$(41)

where its coefficients are mentioned in (A1(A1) $\begin{array}{rl}{\mathrm{\Omega }}_{11}^1& =\frac{(1+h){(x^{\ast }+y^{\ast })}^2+h^{2}s{x}^{\ast 2}{y}^{\ast }-\mathrm{hs}{y}^{\ast 2}}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{12}^1& =\frac{-\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^1& =\frac{\begin{array}{c}-h{(x^{\ast }+y^{\ast })}^3-h^2{(x^{\ast }+y^{\ast })}^3-h^{3}s{x}^{\ast 3}{y}^{\ast }\\ +\mathrm{hs}{y}^{\ast 2}+3h^{2}s{x}^{\ast }{y}^{\ast 2}+h^{2}s{y}^{\ast 3}\end{array}}{{(1+h{x}^{\ast })}^3{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^1& =\frac{h^{2}s{x}^2(x^{\ast }-y^{\ast })-2\mathrm{hs}{x}^{\ast }{y}^{\ast }}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{15}^1& =\frac{\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{11}^2& =\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},{\mathrm{\Omega }}_{12}^2=1-\mathrm{hr}{\delta }^{\ast }+\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^2& =-\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^2& =\frac{2h{\delta }^{\ast }{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},{\mathrm{\Omega }}_{15}^2=-\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3}.\end{array}$(A1) ). Now (41(41) $\begin{array}{r}\{\begin{array}{l}{u}_{t+1}={\mathrm{\Omega }}_{11}^1{u}_t+{\mathrm{\Omega }}_{12}^1{v}_t+{\mathrm{\Omega }}_{13}^1{u}_t^2+{\mathrm{\Omega }}_{14}^1{u}_t{v}_t+{\mathrm{\Omega }}_{15}^1{v}_t^2,\\ v_{t+1}={\mathrm{\Omega }}_{11}^2{u}_t+{\mathrm{\Omega }}_{12}^2{v}_t+{\mathrm{\Omega }}_{13}^2{u}_t^2+{\mathrm{\Omega }}_{14}^2{u}_t{v}_t+{\mathrm{\Omega }}_{15}^2{v}_t^2,\end{array}\end{array}$(41) ) takes the form: (42) $\begin{array}{r}\{\begin{array}{l}{x}_{t+1}=\eta x_t-\zeta y_t+\overline{P}(x_t,y_t),\\ y_{t+1}=\zeta x_t+\eta y_t+\overline{Q}(x_t,y_t),\end{array}\end{array}$(42) by following matrix transformation (43) $\begin{array}{r}\left(\begin{array}{c}{u}_t\\ v_t\end{array}\right):=\left(\begin{array}{ccc}{\mathrm{\Omega }}_{12}^1& 0\eta -{\mathrm{\Omega }}_{11}^1& -\zeta \end{array}\right)\left(\begin{array}{c}{x}_t\\ y_t\end{array}\right),\end{array}$(43) where (44) $\begin{array}{r}\{\begin{array}{l}\overline{P}(x_t,y_t)=r_{11}{x}_t^2+r_{12}{x}_t{y}_t+r_{13}{y}_t^2,\\ \overline{Q}(x_t,y_t)=r_{21}{x}_t^2+r_{22}{x}_t{y}_t+r_{23}{y}_t^2,\end{array}\end{array}$(44) and the quantities η, ζ and $r_{\mathrm{ij}}(i,j=1,2,3)$ are depicted in (A2(A2) $\begin{array}{rl}\eta & =\frac{1}{2}(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr\delta }(r-1)),\\ \zeta & =\frac{1}{2}\sqrt{\begin{array}{c}4(\left({\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}\right)(1+\mathrm{hr\delta }(r-1))+{\displaystyle \frac{h^2{r}^2\mathrm{s\delta }{(r-1)}^2}{1+h+\mathrm{hs}(r-1)}})-{({\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}+1+\mathrm{hr\delta }(r-1))}^2\end{array}},\\ r_{11}& ={\mathrm{\Omega }}_{13}^1{\mathrm{\Omega }}_{12}^1+\eta {\mathrm{\Omega }}_{14}^1-{\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{14}^1-{\eta }^2\frac{{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}+\frac{{\left({\mathrm{\Omega }}_{11}^1\right)}^2{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}-\frac{2\eta {\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{12}& =-\zeta {\mathrm{\Omega }}_{14}^1-\frac{2\mathrm{\zeta \eta }{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}+\frac{2\zeta {\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{13}& ={\zeta }^2\frac{{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{21}& =\frac{1}{\zeta }((\eta -{\mathrm{\Omega }}_{11}^1){\mathrm{\Omega }}_{13}^1{\mathrm{\Omega }}_{12}^1+{(\eta -{\mathrm{\Omega }}_{11}^1)}^2{\mathrm{\Omega }}_{14}^1\frac{{(\eta -{\mathrm{\Omega }}_{11}^1)}^3{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}-{\mathrm{\Omega }}_{13}^2{\left({\mathrm{\Omega }}_{12}^1\right)}^2{\mathrm{\Omega }}_{14}^2{\mathrm{\Omega }}_{12}^1(\eta -{\mathrm{\Omega }}_{11}^1)-{\mathrm{\Omega }}_{15}^2{(\eta -{\mathrm{\Omega }}_{11}^1)}^2),\\ r_{22}& =-(\eta -{\mathrm{\Omega }}_{11}^1){\mathrm{\Omega }}_{14}^1-\frac{2(\eta -{\mathrm{\Omega }}_{11}^1)}{{\mathrm{\Omega }}_{12}^1}{\mathrm{\Omega }}_{15}^1+{\mathrm{\Omega }}_{14}^2{\mathrm{\Omega }}_{12}^1\\ & +2{\mathrm{\Omega }}_{15}^2(\eta -{\mathrm{\Omega }}_{11}^1),\\ r_{23}& =\zeta \frac{{\mathrm{\Omega }}_{15}^1(\eta -{\mathrm{\Omega }}_{11}^1)}{{\mathrm{\Omega }}_{12}^1}-\zeta {\mathrm{\Omega }}_{15}^2.\end{array}$(A2) ). Furthermore, from (44(44) $\begin{array}{r}\{\begin{array}{l}\overline{P}(x_t,y_t)=r_{11}{x}_t^2+r_{12}{x}_t{y}_t+r_{13}{y}_t^2,\\ \overline{Q}(x_t,y_t)=r_{21}{x}_t^2+r_{22}{x}_t{y}_t+r_{23}{y}_t^2,\end{array}\end{array}$(44) ) we get $\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t^2}{|}_{{\phi }_2=(0,0)}=2r_{11}$, $\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}{|}_{{\phi }_2=(0,0)}=r_{12}$, $\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{y}_t^2}{|}_{{\phi }_2=(0,0)}=2r_{13}$, $\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{x}_t^3}{|}_{{\phi }_2=(0,0)}=0$, $\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{x}_t^2\mathrm{\partial }{y}_t}{|}_{{\phi }_2=(0,0)}=0$, $\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t^2}{|}_{{\phi }_2=(0,0)}=0$, $\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{y}_t^3}{|}_{{\phi }_2=(0,0)}=0$, $\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t^2}{|}_{{\phi }_2=(0,0)}=2r_{21}$, $\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}{|}_{{\phi }_2=(0,0)}=r_{22}$, $\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{y}_t^2}{|}_{{\phi }_2=(0,0)}=2r_{23}$, $\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t^3}{|}_{{\phi }_2=(0,0)}=\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t^2\mathrm{\partial }{y}_t}{|}_{{\phi }_2=(0,0)}=\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t^2}{|}_{{\phi }_2=(0,0)}=\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{y}_t^3}{|}_{{\phi }_2=(0,0)}=0$. Now the following discriminatory quantity should not be zero in order to guarantee the occurrence of Neimark–Sacker bifurcation at ${\phi }_2$ of system (42(42) $\begin{array}{r}\{\begin{array}{l}{x}_{t+1}=\eta x_t-\zeta y_t+\overline{P}(x_t,y_t),\\ y_{t+1}=\zeta x_t+\eta y_t+\overline{Q}(x_t,y_t),\end{array}\end{array}$(42) ): (45) $\begin{array}{rl}\ell & :=-\mathrm{\Re }\left(\frac{(1-2\overline{\lambda }){\overline{\lambda }}^2}{1-\lambda }{\varrho }_{11}{\varrho }_{20}\right)-\frac{1}{2}\Vert {\varrho }_{11}{\Vert }^2\\ & -\Vert {\varrho }_{02}{\Vert }^2+\mathrm{\Re }\left(\overline{\lambda }{\varrho }_{21}\right),\end{array}$(45) where (46) $\begin{array}{rl}{\varrho }_{02}& =\frac{1}{8}(\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t^2}-\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{y}_t^2}+2\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}\\ & {+\iota (\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t^2}-\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{y}_t^2}+2\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}))|}_{{\phi }_2=(0,0)},\\ {\varrho }_{11}& =\frac{1}{4}{(\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t^2}+\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{y}_t^2}+\iota (\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t^2}+\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{y}_t^2}))|}_{{\phi }_2=(0,0)},\\ {\varrho }_{20}& =\frac{1}{8}(\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t^2}-\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{y}_t^2}+2\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}\\ & {+\iota (\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t^2}-\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{y}_t^2}-2\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}))|}_{{\phi }_2=(0,0)},\\ {\varrho }_{21}& =\frac{1}{16}(\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{x}_t^3}+\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{y}_t^3}+\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t^2\mathrm{\partial }{y}_t}+\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{y}_t^3}\\ & {+\iota (\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t^3}+\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t^2}-\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{x}_t^2\mathrm{\partial }{y}_t}-\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{y}_t^3}))|}_{{\phi }_2=(0,0)}.\end{array}$(46) Using partial derivatives, (46(46) $\begin{array}{rl}{\varrho }_{02}& =\frac{1}{8}(\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t^2}-\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{y}_t^2}+2\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}\\ & {+\iota (\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t^2}-\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{y}_t^2}+2\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}))|}_{{\phi }_2=(0,0)},\\ {\varrho }_{11}& =\frac{1}{4}{(\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t^2}+\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{y}_t^2}+\iota (\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t^2}+\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{y}_t^2}))|}_{{\phi }_2=(0,0)},\\ {\varrho }_{20}& =\frac{1}{8}(\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t^2}-\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{y}_t^2}+2\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}\\ & {+\iota (\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{x}_t^2}-\frac{{\mathrm{\partial }}^2\overline{Q}}{\mathrm{\partial }{y}_t^2}-2\frac{{\mathrm{\partial }}^2\overline{P}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t}))|}_{{\phi }_2=(0,0)},\\ {\varrho }_{21}& =\frac{1}{16}(\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{x}_t^3}+\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{y}_t^3}+\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t^2\mathrm{\partial }{y}_t}+\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{y}_t^3}\\ & {+\iota (\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t^3}+\frac{{\mathrm{\partial }}^3\overline{Q}}{\mathrm{\partial }{x}_t\mathrm{\partial }{y}_t^2}-\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{x}_t^2\mathrm{\partial }{y}_t}-\frac{{\mathrm{\partial }}^3\overline{P}}{\mathrm{\partial }{y}_t^3}))|}_{{\phi }_2=(0,0)}.\end{array}$(46) ) becomes (47) $\begin{array}{rl}{\varrho }_{02}& =\frac{1}{4}(r_{11}-r_{13}+r_{22}+\iota (r_{21}-r_{23}+r_{12})),\\ {\varrho }_{11}& =\frac{1}{2}(r_{11}+r_{13}+\iota (r_{21}+r_{23})),\\ {\varrho }_{20}& =\frac{1}{4}(r_{11}-r_{13}+r_{22}+\iota (r_{21}-r_{23}-r_{12})),\\ {\varrho }_{21}& =0.\end{array}$(47) From (45(45) $\begin{array}{rl}\ell & :=-\mathrm{\Re }\left(\frac{(1-2\overline{\lambda }){\overline{\lambda }}^2}{1-\lambda }{\varrho }_{11}{\varrho }_{20}\right)-\frac{1}{2}\Vert {\varrho }_{11}{\Vert }^2\\ & -\Vert {\varrho }_{02}{\Vert }^2+\mathrm{\Re }\left(\overline{\lambda }{\varrho }_{21}\right),\end{array}$(45) ) the condition $\ell \ne 0$ holds as $(h,s,\delta ,r)\in \mathcal{N}{|}_{{\phi }_2}$ then Neimark–Sacker bifurcation takes place, and additionally supercritical (subcritical) Neimark–Sacker bifurcation take place if $\ell <0$ ($\ell >0$).

5. Analysis of a flip bifurcation for interior fixed point

Based on the local dynamic study for the equilibrium ${\phi }_2$ in Section 2, we investigate the flip bifurcation if $(h,s,\delta ,r)\in \mathcal{F}{|}_{{\phi }_2}$ where δ is regarded as a bifurcation parameter.

Theorem 5.1

At ${\phi }_2$, discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) undergoes flip bifurcation if $(h,s,\delta ,r)\in \mathcal{F}{|}_{{\phi }_2}$.

Proof.

Recall that if δ is a bifurcation parameter then the discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) takes the form (35(35) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+h({\delta }^{\ast }+ϵ)y_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(35) ), and by (39(39) $\begin{array}{r}{u}_t=x_t-x^{\ast },v_t=y_t-y^{\ast },\end{array}$(39) ) it further becomes (48) $\begin{array}{rl}{u}_{t+1}& ={\mathrm{\Omega }}_{11}^1{u}_t+{\mathrm{\Omega }}_{12}^1{v}_t+{\mathrm{\Omega }}_{13}^1{u}_t^2+{\mathrm{\Omega }}_{14}^1{u}_t{v}_t+{\mathrm{\Omega }}_{15}^1{v}_t^2,\\ v_{t+1}& ={\mathrm{\Omega }}_{11}^2{u}_t+{\mathrm{\Omega }}_{12}^2{v}_t+{\mathrm{\Omega }}_{13}^2{u}_t^2+{\mathrm{\Omega }}_{14}^2{u}_t{v}_t+{\mathrm{\Omega }}_{15}^2{v}_t^2\\ & +k_1{u}_tϵ+k_2{v}_tϵ+k_3{u}_t^2ϵ\\ & +k_4{u}_t{v}_tϵ+k_5{v}_t^2ϵ,\end{array}$(48) where it coefficients are depicted in (A3(A3) $\begin{array}{rl}{\mathrm{\Omega }}_{11}^1& =\frac{(1+h){(x^{\ast }+y^{\ast })}^2+h^{2}s{x}^{\ast 2}{y}^{\ast }-\mathrm{hs}{y}^{\ast 2}}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{12}^1& =\frac{-\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^1& =\frac{\begin{array}{c}-h{(x^{\ast }+y^{\ast })}^3-h^2{(x^{\ast }+y^{\ast })}^3-h^{3}s{x}^{\ast 3}{y}^{\ast }+\mathrm{hs}{y}^{\ast 2}\\ +3h^{2}s{x}^{\ast }{y}^{\ast 2}+h^{2}s{y}^{\ast 3}\end{array}}{{(1+h{x}^{\ast })}^3{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^1& =\frac{h^{2}s{x}^2(x^{\ast }-y^{\ast })-2\mathrm{hs}{x}^{\ast }{y}^{\ast }}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{15}^1& =\frac{\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{11}^2& =\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},{\mathrm{\Omega }}_{12}^2=1-\mathrm{hr}{\delta }^{\ast }+\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^2& =-\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^2& =\frac{2h{\delta }^{\ast }{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},{\mathrm{\Omega }}_{15}^2=-\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},k_1=\frac{h{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ k_2& =-\mathrm{hr}+\frac{h{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},k_3=\frac{-h{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ k_4& =\frac{2h{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},k_5=\frac{-h{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3}.\end{array}$(A3) ). Now using following transformation: (49) $\begin{array}{rl}\left(\begin{array}{c}{u}_t\\ v_t\end{array}\right)& :=(\begin{array}{c}1\\ {\displaystyle \frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}}\end{array}\\ & \begin{array}{c}1\\ {\displaystyle \frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}}\end{array})\left(\begin{array}{c}{x}_t\\ y_t\end{array}\right),\end{array}$(49) Equation (48(48) $\begin{array}{rl}{u}_{t+1}& ={\mathrm{\Omega }}_{11}^1{u}_t+{\mathrm{\Omega }}_{12}^1{v}_t+{\mathrm{\Omega }}_{13}^1{u}_t^2+{\mathrm{\Omega }}_{14}^1{u}_t{v}_t+{\mathrm{\Omega }}_{15}^1{v}_t^2,\\ v_{t+1}& ={\mathrm{\Omega }}_{11}^2{u}_t+{\mathrm{\Omega }}_{12}^2{v}_t+{\mathrm{\Omega }}_{13}^2{u}_t^2+{\mathrm{\Omega }}_{14}^2{u}_t{v}_t+{\mathrm{\Omega }}_{15}^2{v}_t^2\\ & +k_1{u}_tϵ+k_2{v}_tϵ+k_3{u}_t^2ϵ\\ & +k_4{u}_t{v}_tϵ+k_5{v}_t^2ϵ,\end{array}$(48) ) gives (50) $\begin{array}{r}\left(\begin{array}{c}{x}_{t+1}\\ y_{t+1}\end{array}\right)=\left(\begin{array}{cc}-1& 0\\ 0& {\lambda }_2\end{array}\right)\left(\begin{array}{c}{x}_t\\ y_t\end{array}\right)+\left(\begin{array}{c}\hat{P}(u_t,v_t,ϵ)\\ \hat{Q}(u_t,v_t,ϵ)\end{array}\right),\end{array}$(50) where (51) $\begin{array}{rl}\hat{P}& =\left\{\begin{array}{rl}& \frac{2\mathrm{hrs}(r-1)(1+h+\mathrm{hsr}-\mathrm{hs})}{4+h(1+s(r-1))(4+h+\mathrm{hs}{r}^2-\mathrm{hs})}\\ & ({\mathrm{\Omega }}_{13}^1{u}_t^2+{\mathrm{\Omega }}_{14}^1{u}_t{v}_t+{\mathrm{\Omega }}_{15}^1{v}_t^2)\\ & +\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{\begin{array}{c}4+h(1+s(-1+r))\\ (4+h+hs(r^2-1))\end{array}}\\ & \times ({\mathrm{\Omega }}_{13}^2{u}_t^2+{\mathrm{\Omega }}_{14}^2{u}_t{v}_t+{\mathrm{\Omega }}_{15}^2{v}_t^2+k_1ϵu_t\\ & +k_2ϵv_t+k_3ϵu_t^2+k_4ϵu_t{v}_t+k_5ϵv_t^2)\end{array}\right\},\\ \hat{Q}& =\left\{\begin{array}{rl}& \frac{\begin{array}{c}-(2+h(1+\mathrm{sr}-s))\\ (-2+h(-1+s(r-1{)}^2))\end{array}}{4+h(1+s(r-1))(4+h+\mathrm{hs}(r^2-1))}\\ & ({\mathrm{\Omega }}_{13}^1{u}_t^2+{\mathrm{\Omega }}_{14}^1{u}_t{v}_t+{\mathrm{\Omega }}_{15}^1{v}_t^2)\\ & -\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{\begin{array}{c}4+h(1+s(-1+r))\\ (4+h+hs(r^2-1))\end{array}}\\ & \times ({\mathrm{\Omega }}_{13}^2{u}_t^2+{\mathrm{\Omega }}_{14}^2{u}_t{v}_t+{\mathrm{\Omega }}_{15}^2{v}_t^2+k_1ϵu_t\\ & +k_2ϵv_t+k_3ϵu_t^2+k_4ϵu_t{v}_t+k_5ϵv_t^2)\end{array}\right\},\end{array}$(51)

Now we calculate the centre manifold $F^{C}O$ at O in a small neighborhood of $ϵ=0$ where mathematically it can be expressed as (53) $\begin{array}{rl}{F}^{C}O& =\{(x_t,y_t):y_t=v_0ϵ+v_1{x}_t^2+v_2ϵx_t+v_3{ϵ}^3\\ & +O\left({(|x_t|+|ϵ|)}^3\right)\},\end{array}$(53) with (54) $\begin{array}{r}\{\begin{array}{l}{v}_1=\left\{\begin{array}{rl}& {\displaystyle \frac{1}{1-{\lambda }_2}(\frac{\begin{array}{c}-(2+h(1+\mathrm{sr}-s))\\ (-2+h(-1+s(r-1{)}^2))\end{array}}{\begin{array}{c}4+h(1+s(r-1))\\ (4+h+hs(r^2-1))\end{array}}}\\ & {\displaystyle \times ({\mathrm{\Omega }}_{13}^1+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}{\mathrm{\Omega }}_{14}^1}\\ & +{\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{\mathrm{\Omega }}_{15}^1)\\ & {\displaystyle -\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{\begin{array}{c}4+h(1+s(r-1))\\ (4+h+\mathrm{hs}(r^2-1))\end{array}}}\\ & {\displaystyle \times ({\mathrm{\Omega }}_{13}^2+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}{\mathrm{\Omega }}_{14}^2}\\ & +{\left(\frac{2-\mathrm{hs}(-1+h{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{\mathrm{\Omega }}_{15}^2))\end{array}\right\},\\ v_2=\left\{\begin{array}{r}{\displaystyle \frac{1}{1-{\lambda }_2}}\\ (\frac{-h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{\begin{array}{c}4+h(1+s(r-1))\\ (4+h+\mathrm{hs}(r^2-1))\end{array}}\\ {\displaystyle \times (k_1+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}{k}_2))}\end{array}\right\},\\ v_0=v_3=0.\end{array}\end{array}$(54) Finally, one write (50(50) $\begin{array}{r}\left(\begin{array}{c}{x}_{t+1}\\ y_{t+1}\end{array}\right)=\left(\begin{array}{cc}-1& 0\\ 0& {\lambda }_2\end{array}\right)\left(\begin{array}{c}{x}_t\\ y_t\end{array}\right)+\left(\begin{array}{c}\hat{P}(u_t,v_t,ϵ)\\ \hat{Q}(u_t,v_t,ϵ)\end{array}\right),\end{array}$(50) ) restrict to $F^{c}O$ as (55) $\begin{array}{rl}{f}_1(x_t)& =-x_t+m_1{x}_t^2+m_2{x}_tϵ+m_3{x}_t^2ϵ+m_4{x}_t{ϵ}^2\\ & +m_5{x}_t^3+O\left({(|x_t|+|ϵ|)}^4\right),\end{array}$(55) where

(56) $\begin{array}{r}\{\begin{array}{l}{m}_1=\left\{\begin{array}{rl}& {\displaystyle \frac{2\mathrm{hrs}(r-1)(1+h+\mathrm{hs}(r-1))}{4+h(1+s(r-1))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle \times ({\mathrm{\Omega }}_{13}^1+\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right){\mathrm{\Omega }}_{14}^1{\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{\mathrm{\Omega }}_{15}^1)}\\ & {\displaystyle +\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle \times ({\mathrm{\Omega }}_{13}^2+\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right){\mathrm{\Omega }}_{14}^2{\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{\mathrm{\Omega }}_{15}^2)}\end{array}\right\},\\ m_2={\displaystyle \frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}}(k_1+{\displaystyle \frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}}{k}_2),\\ m_3=\left\{\begin{array}{rl}& {\displaystyle \frac{2\mathrm{hrs}(r-1)(1+h+\mathrm{hs}(r-1))}{4+h(1+s(r-1))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle \times (2{\mathrm{\Omega }}_{13}^1{v}_2+{\mathrm{\Omega }}_{14}^1{v}_2(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2})}\\ & {\displaystyle +{\mathrm{\Omega }}_{15}^1{v}_2\left(2\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)\left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right)\right))}\\ & {\displaystyle +\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}(2{\mathrm{\Omega }}_{13}^2{v}_2{\displaystyle +{\mathrm{\Omega }}_{14}^2{v}_2(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}}}\\ \\ & +\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2})\\ & {\displaystyle +{\mathrm{\Omega }}_{15}^2{v}_2(2\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}\\ & \left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right))\\ & {\displaystyle +k_1{v}_1+\left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right)k_2{v}_1}\\ & {\displaystyle +k_3+\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)k_4}\\ & +{\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{k}_5)\end{array}\right\},\\ m_4=\left\{\begin{array}{rl}& {\displaystyle \frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle (k_1{v}_2+\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}{k}_2{v}_2)}\end{array}\right\},\\ m_5=\left\{\begin{array}{rl}& {\displaystyle \frac{2\mathrm{hrs}(r-1)(1+h+\mathrm{hs}(r-1))}{4+h(1+s(r-1))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle \times (2{\mathrm{\Omega }}_{13}^1{v}_1+{\mathrm{\Omega }}_{14}^1{v}_1(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2})}\\ & {\displaystyle +{\mathrm{\Omega }}_{15}^1{v}_1\left(2\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)\left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right))\right)}\\ & {\displaystyle +\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}(2{\mathrm{\Omega }}_{13}^2{v}_1\phantom{\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}}}\\ & {\displaystyle +{\mathrm{\Omega }}_{14}^2{v}_1(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2})}\\ & {\displaystyle +{\mathrm{\Omega }}_{15}^2{v}_1\left(2\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)\left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right))\right)}\end{array}\right\}.\end{array}\end{array}$(56)

Finally, in order for the existence of flip bifurcation at ${\phi }_2$, the following two discriminatory quantities ${\ell }_1$ and ${\ell }_2$ should not be zero [18, 19]: (57) $\begin{array}{rl}{\ell }_1& :=(\frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_t\mathrm{\partial }ϵ}+\frac{1}{2}\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }ϵ}\frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_t^2}){|}_{(0,0)}=m_2,\\ {\ell }_2& :=(\frac{1}{6}\frac{{\mathrm{\partial }}^3{f}_1}{\mathrm{\partial }{x}_t^3}+{\left(\frac{1}{2}\frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_t^2}\right)}^2){|}_{(0,0)}=m_5+m_1^2.\end{array}$(57) Consequently, based on above calculation, we can say that if $(h,s,\delta ,r)\in \mathcal{F}{|}_{{\phi }_2}$ and ${\ell }_2\ne 0$ then at ${\phi }_2$ flip bifurcation takes place. Besides, the period-2 points bifurcate from ${\phi }_2$ of model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) are stable (unstable) if ${\ell }_2>0$ (${\ell }_2<0$).

6. Chaos control

In the present section, we will apply the following two feedback control strategies to the model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) to get a stable trajectory:

6.1. By OGY method

In order to study chaos, we fist use OGY method, which was proposed by Ott et al. [28], for the discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ). By existing chaos theory [29, 30], one can write the model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) as follows: (58) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)}=f_1(x_t,y_t,s),{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t})=f_2(x_t,y_t,s),}}\end{array}\end{array}$(58) where s represents the control parameter under which one can acquire the desired chaos control by small perturbations. To do this, one restrict it where $s\in (s_0-\sigma ,s_0+\sigma )$ with $\sigma >0$, and $s_0$ indicates nominal value corresponds to chaotic region. The trajectory is moved in the direction of the target orbit using the stabilizing feedback control strategy. Assuming that ${\phi }_2$ be unstable equilibrium point of the model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ), which is in the chaotic region created by the occurrence of N-S bifurcation, then the given below linear map can be used to approximate the model (58(58) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)}=f_1(x_t,y_t,s),{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t})=f_2(x_t,y_t,s),}}\end{array}\end{array}$(58) ): (59) $\begin{array}{r}\left(\begin{array}{c}{x}_{t+1}-x\\ y_{t+1}-y\end{array}\right)\approx V(x,y,s_0)\left(\begin{array}{c}{x}_t-x\\ y_t-y\end{array}\right)+P(s-s_0),\end{array}$(59) where (60) $\begin{array}{rlr}V{|}_{(x,y,s_0)}& =\left(\begin{array}{ccc}{\displaystyle \frac{\mathrm{\partial }{f}_1(x,y,s_0)}{\mathrm{\partial }x}}& {\displaystyle \frac{\mathrm{\partial }{f}_1(x,y,s_0)}{\mathrm{\partial }y}{\displaystyle \frac{\mathrm{\partial }{f}_2(x,y,s_0)}{\mathrm{\partial }x}}}& {\displaystyle \frac{\mathrm{\partial }{f}_2(x,y,s_0)}{\mathrm{\partial }y}}\end{array}\right)& =\left(\begin{array}{ccc}{\displaystyle \frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}}& {\displaystyle \frac{-h{r}^2{s}_0}{1+h+h{s}_0(r-1)}\mathrm{h\delta }(r-1{)}^2}& 1+\mathrm{hr\delta }(r-1)\end{array}\right),\end{array}$(60) and (61) $\begin{array}{r}P=\left(\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }{f}_1(x,y,s_0)}{\mathrm{\partial }s}{\displaystyle \frac{\mathrm{\partial }{f}_2(x,y,s_0)}{\mathrm{\partial }s}}}\end{array}\right)=\left(\begin{array}{c}{\displaystyle \frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}0}\end{array}\right).\end{array}$(61) Moreover, the model (58(58) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)}=f_1(x_t,y_t,s),{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t})=f_2(x_t,y_t,s),}}\end{array}\end{array}$(58) ) is controlled provided that the following matrix (62) $\begin{array}{rl}C& =\left(\begin{array}{c}P:\mathrm{VP}\end{array}\right)=(\begin{array}{c}{\displaystyle \frac{h(-1+r)(1-s_0+s_{0}r)}{1+h+h{s}_0(-1+r)}0}\end{array}\\ & \begin{array}{c}{\displaystyle \frac{h(-1+r)(1-\mathrm{hr}{s}_0(-1+r))(1-s_0+s_{0}r)}{(1+h+h{s}_0(-1+r){)}^2}{\displaystyle \frac{h^2\delta (-1+r{)}^3(1-s_0+s_{0}r)}{1+h+h{s}_0(-1+r)}}}\end{array}),\end{array}$(62) is of rank 2. Additionally, it must be understood that δ cannot be used in this situation as a control parameter to utilize the OGY method to the model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ). As, if we take $f_1(x_t,y_t,\delta )=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)}$ and $f_2(x_t,y_t,\delta )=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t})$, then $\mathrm{Rank}(C)$ becomes zero. So, to apply the OGY method to the model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ), s is taken as a control parameter. Now, if one take $s-s_0=-T(\begin{array}{c}{x}_t-x\\ y_t-y\end{array})$, where $T=(k_1{k}_2)$, then the model (59(59) $\begin{array}{r}\left(\begin{array}{c}{x}_{t+1}-x\\ y_{t+1}-y\end{array}\right)\approx V(x,y,s_0)\left(\begin{array}{c}{x}_t-x\\ y_t-y\end{array}\right)+P(s-s_0),\end{array}$(59) ) can be written as (63) $\begin{array}{r}\left(\begin{array}{c}{x}_{t+1}-x\\ y_{t+1}-y\end{array}\right)\approx (V-\mathrm{PT})\left(\begin{array}{c}{x}_t-x\\ y_t-y\end{array}\right).\end{array}$(63) Furthermore, the corresponding controlled model of (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) is (64) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-{\displaystyle \frac{h(s_0-k_1(x_t-x)-k_2(y_t-y))x_t{y}_t}{(x_t+y_t)(1+h{x}_t)}},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}).}}\end{array}\end{array}$(64) By stability theory ${\phi }_2$ is a sink iff roots of V −PT satisfying $|{\lambda }_{1,2}|<1$ where its variational matrix is (65) $\begin{array}{rl}V-\mathrm{PT}& =(\begin{array}{c}{\displaystyle \frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}-{\displaystyle \frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}}{k}_1\mathrm{h\delta }(r-1{)}^2}\end{array}\\ & \begin{array}{c}{\displaystyle \frac{-h{r}^2{s}_0}{1+h+h{s}_0(r-1)}-{\displaystyle \frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}}{k}_{2}1+\mathrm{hr\delta }(r-1)}\end{array}).\end{array}$(65) The characteristic equation of V −PT is (66) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\lambda }^2-(\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}-\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1+1+\mathrm{hr\delta }(r-1))\lambda }\\ {\displaystyle +(\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}-{\displaystyle \frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}}{k}_1)(1+\mathrm{hr\delta }(r-1))}\\ {\displaystyle +\mathrm{h\delta }(r-1{)}^2(\frac{h{r}^2{s}_0}{1+h+h{s}_0(r-1)}+{\displaystyle \frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}}{k}_2)=0.}\end{array}\end{array}$(66) If ${\lambda }_{1,2}$ denotes the roots of (66(66) $\begin{array}{r}\{\begin{array}{l}{\displaystyle {\lambda }^2-(\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}-\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1+1+\mathrm{hr\delta }(r-1))\lambda }\\ {\displaystyle +(\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}-{\displaystyle \frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}}{k}_1)(1+\mathrm{hr\delta }(r-1))}\\ {\displaystyle +\mathrm{h\delta }(r-1{)}^2(\frac{h{r}^2{s}_0}{1+h+h{s}_0(r-1)}+{\displaystyle \frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}}{k}_2)=0.}\end{array}\end{array}$(66) ) then (67) $\begin{array}{rl}{\lambda }_1+{\lambda }_2& =\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}\\ & -\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1+1+\mathrm{hr\delta }(r-1),\end{array}$(67) and (68) $\begin{array}{rl}{\lambda }_1{\lambda }_2& =(\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}\\ & -\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1)(1+\mathrm{hr\delta }(r-1))\\ & +\mathrm{h\delta }(r-1{)}^2(\frac{h{r}^2{s}_0}{1+h+h{s}_0(r-1)}\\ & +\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_2).\end{array}$(68) Next we take ${\lambda }_1=\pm 1$ and ${\lambda }_1{\lambda }_2=1$ to get lines of marginal stability for (64(64) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-{\displaystyle \frac{h(s_0-k_1(x_t-x)-k_2(y_t-y))x_t{y}_t}{(x_t+y_t)(1+h{x}_t)}},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}).}}\end{array}\end{array}$(64) ). Also, these restriction give the fact that characteristics roots satisfying $|{\lambda }_{1,2}|<1$. Assuming that ${\lambda }_1{\lambda }_2=1$ then from (68(68) $\begin{array}{rl}{\lambda }_1{\lambda }_2& =(\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}\\ & -\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1)(1+\mathrm{hr\delta }(r-1))\\ & +\mathrm{h\delta }(r-1{)}^2(\frac{h{r}^2{s}_0}{1+h+h{s}_0(r-1)}\\ & +\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_2).\end{array}$(68) ) one has (69) $\begin{array}{rl}& L_1:h(r-1)(1-s_0+s_{0}r)(1+\mathrm{hr\delta }(r-1))k_1\\ & -h^2\delta (r-1{)}^3(1-s_0+s_{0}r)k_2\\ & +h-\mathrm{hr\delta }(r-1)+h{s}_0(r^2-1)=0.\end{array}$(69) If ${\lambda }_1=1$ then from (67(67) $\begin{array}{rl}{\lambda }_1+{\lambda }_2& =\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}\\ & -\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1+1+\mathrm{hr\delta }(r-1),\end{array}$(67) ) and (68(68) $\begin{array}{rl}{\lambda }_1{\lambda }_2& =(\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}\\ & -\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1)(1+\mathrm{hr\delta }(r-1))\\ & +\mathrm{h\delta }(r-1{)}^2(\frac{h{r}^2{s}_0}{1+h+h{s}_0(r-1)}\\ & +\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_2).\end{array}$(68) ), one gets (70) $\begin{array}{rl}& L_2:h^2\mathrm{r\delta }(r-1{)}^2(1-s_0+s_{0}r)k_1\\ & -h^2\delta (r-1{)}^3(1-s_0+s_{0}r)k_2\\ & +h^2\mathrm{r\delta }(r-1)(1-s_0+s_{0}r))=0.\end{array}$(70) Finally, if ${\lambda }_1=-1$ then from (67(67) $\begin{array}{rl}{\lambda }_1+{\lambda }_2& =\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}\\ & -\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1+1+\mathrm{hr\delta }(r-1),\end{array}$(67) ) and (68(68) $\begin{array}{rl}{\lambda }_1{\lambda }_2& =(\frac{1-\mathrm{hr}{s}_0(r-1)}{1+h+h{s}_0(r-1)}\\ & -\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_1)(1+\mathrm{hr\delta }(r-1))\\ & +\mathrm{h\delta }(r-1{)}^2(\frac{h{r}^2{s}_0}{1+h+h{s}_0(r-1)}\\ & +\frac{h(r-1)(1-s_0+s_{0}r)}{1+h+h{s}_0(r-1)}{k}_2).\end{array}$(68) ), one gets (71) $\begin{array}{rl}& L_3:h(r-1)(1-s_0+s_{0}r)(2+\mathrm{hr\delta }(r-1))k_1\\ & -h^2\delta (r-1{)}^3(1-s_0+s_{0}r)k_2-4\\ & -2h+2h{s}_0(r-1{)}^2\\ & -\mathrm{hr\delta }(-1+r)(2+h(1+s_{0}r-s_0))=0.\end{array}$(71) Then, stable eigenvalues lie within the triangular region $(k_1{k}_2)$-plane bounded by the straight lines $L_{1,2,3}$ for parametric values h, s, r and δ.

6.2. By hybrid control feedback

We use a hybrid control feedback method to control the chaos due to the emergence of flip bifurcation in the model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) by existing theory [31]. If the model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) undergoes bifurcation at ${\phi }_2$ then one can write corresponding controlled model as (72) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\alpha (\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)})+(1-\alpha )x_t,}\\ {\displaystyle y_{t+1}=\alpha (y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}))+(1-\alpha )y_t,}\end{array}\end{array}$(72) where $0<\alpha <1$ and controlled strategy (72(72) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\alpha (\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)})+(1-\alpha )x_t,}\\ {\displaystyle y_{t+1}=\alpha (y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}))+(1-\alpha )y_t,}\end{array}\end{array}$(72) ) is a combination of feedback control as well as parameter perturbation. The $V{|}_{{\phi }_2}$ is (73) $\begin{array}{r}V{|}_{{\phi }_2}=\left(\begin{array}{cc}{\displaystyle 1+\frac{\mathrm{\alpha h}(-1+s-s{r}^2)}{1+h+\mathrm{hs}(r-1)}}& {\displaystyle \frac{-\mathrm{\alpha h}{r}^{2}s}{1+h+\mathrm{hs}(r-1)}}\\ \mathrm{\alpha h\delta }(r-1{)}^2& 1+\mathrm{\alpha hr\delta }(r-1)\end{array}\right).\end{array}$(73) The characteristics equation of $V{|}_{{\phi }_2}$ evaluated at ${\phi }_2$ is (74) $\begin{array}{r}{\lambda }^2-{\mathrm{\Lambda }}_1\lambda +{\mathrm{\Lambda }}_2=0\end{array}$(74) where (75) $\begin{array}{rl}{\mathrm{\Lambda }}_1& =2+\frac{\mathrm{\alpha h}(-1+s-s{r}^2)}{1+h+\mathrm{hs}(r-1)}+\mathrm{\alpha h\delta r}(r-1),\\ {\mathrm{\Lambda }}_2& =(1+\frac{\mathrm{\alpha h}(-1+s-s{r}^2)}{1+h(1+s(r-1))})(1+\mathrm{\alpha hr\delta }(r-1))\\ & +\frac{{\alpha }^2{h}^2{r}^2\mathrm{s\delta }(r-1{)}^2}{1+h(1+\mathrm{hs}(r-1))}.\end{array}$(75) So, based on stability theory one has the result:

Lemma 1

Equilibrium ${\phi }_2$ of model (72(72) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\alpha (\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)})+(1-\alpha )x_t,}\\ {\displaystyle y_{t+1}=\alpha (y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}))+(1-\alpha )y_t,}\end{array}\end{array}$(72) ) is a sink iff $|2+\frac{\mathrm{\alpha h}(-1+s-s{r}^2)}{1+h+\mathrm{hs}(r-1)}+\mathrm{\alpha h\delta r}(r-1)|<1+(1+\frac{\mathrm{\alpha h}(-1+s-s{r}^2)}{1+h+\mathrm{hs}(r-1)})(1+\mathrm{\alpha hr\delta }(r-1))+\frac{{\alpha }^2{h}^2{r}^2\mathrm{s\delta }(r-1{)}^2}{1+h+\mathrm{hs}(r-1)}<2$.

Proof.

The necessary and sufficient condition for roots of (74(74) $\begin{array}{r}{\lambda }^2-{\mathrm{\Lambda }}_1\lambda +{\mathrm{\Lambda }}_2=0\end{array}$(74) ) satisfying $|{\lambda }_{1,2}|$ imply that ${\phi }_2$ of model (72(72) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\alpha (\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)})+(1-\alpha )x_t,}\\ {\displaystyle y_{t+1}=\alpha (y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}))+(1-\alpha )y_t,}\end{array}\end{array}$(72) ) is a sink iff $|2+\frac{\mathrm{\alpha h}(-1+s-s{r}^2)}{1+h+\mathrm{hs}(r-1)}+\mathrm{\alpha h\delta r}(r-1)|<1+(1+\frac{\mathrm{\alpha h}(-1+s-s{r}^2)}{1+h+\mathrm{hs}(r-1)})(1+\mathrm{\alpha hr\delta }(r-1))+\frac{{\alpha }^2{h}^2{r}^2\mathrm{s\delta }(r-1{)}^2}{1+h+\mathrm{hs}(r-1)}<2$.

7. Numerical simulations

Example 7.1

If h = 0.6, s = 1.37, r = 0.5, $\delta =[0.1,1.1]$ with $(x_0,y_0)=(0.03,0.04)$ then at $\delta =0.11000000000000032$ discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) undergoes the N-S bifurcation. The bifurcation diagrams and Maximum Lyapunov exponent are drawn in Figure 1. Further, at $(h,s,r,\delta )=(0.6,1.37,0.5,0.11000000000000032)$ model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) has interior equilibrium point ${\phi }_2=(0.31499999999999995,0.31499999999999995)$ and moreover from (24(24) $\begin{array}{r}V{|}_{{\phi }_2}:=\left(\begin{array}{ccc}{\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}& {\displaystyle \frac{-h{r}^{2}s}{1+h+\mathrm{hs}(r-1)}\mathrm{h\delta }(r-1{)}^2}& 1+\mathrm{hr\delta }(r-1)\end{array}\right),\end{array}$(24) ) one gets: (76) $\begin{array}{rl}V{|}_{{\phi }_2}& =(0.31499999999999995,0.31499999999999995)\\ & =(\begin{array}{c}1.0138772077375944\\ 0.016500000000000046\end{array}\\ & \begin{array}{c}-0.17283431455004208\\ 0.9834999999999999\end{array}),\end{array}$(76) with (77) $\begin{array}{r}{\lambda }^2-1.9973772077375944\lambda +1=0.\end{array}$(77) The roots of (77(77) $\begin{array}{r}{\lambda }^2-1.9973772077375944\lambda +1=0.\end{array}$(77) ) are ${\lambda }_{1,2}=0.9986886038687972\pm 0.05119641103234161\iota$ with $|{\lambda }_{1,2}|=1$, and so from (i) of Theorem 2.5 one has the Neimark–Sacker bifurcation set $(h,s,r,\delta )=(0.6,1.37,0.5,0.11000000000000032)\in \mathcal{N}{|}_{{\phi }_2=(0.31499999999999995,0.31499999999999995)}$. Furthermore, by fixing h = 0.6, s = 1.37, r = 0.5, and varying the value of bifurcation parameter $\delta <0.11000000000000032$ then one can obtain that respective interior equilibrium solution is a stable focus. To verify this if $\delta =0.10000<0.11000000000000032$ then Figure 2(a) indicates that ${\phi }_2=(0.31499999999999995,0.31499999999999995)$ is a stable focus. Similarly if one varies $\delta =0.10299$, 0.10500, 0.10880, 0.10987, 0.10999<0.11000000000000032 then Figure 2(b–f) also show that, for each variation of δ, the equilibrium solution ${\phi }_2=(0.31499999999999995,0.31499999999999995)$ is also stable focus. On the other hand, if one varies the bifurcation parameter $\delta >0.11000000000000032$ then we can concluded that respective interior equilibrium solution is an unstable focus, and as a result supercritical Neimark–Sacker bifurcation must take place. For instance, if h = 0.6, s = 1.37, r = 0.5 then from (37(37) $\begin{array}{r}\frac{d|{\lambda }_{1,2}|}{\mathrm{dϵ}}{|}_{ϵ=0}=\frac{\mathrm{hr}(r-1)}{2(1+h+\mathrm{hs}(r-1))}\ne 0.\end{array}$(37) ) the non-degenerate conditions, i.e. $\frac{d|{\lambda }_{1,2}|}{\mathrm{dϵ}}{|}_{ϵ=0}=-0.06307821698906643\ne 0$ holds, and moreover if $\delta =0.11000000000000032$ then from (A1(A1) $\begin{array}{rl}{\mathrm{\Omega }}_{11}^1& =\frac{(1+h){(x^{\ast }+y^{\ast })}^2+h^{2}s{x}^{\ast 2}{y}^{\ast }-\mathrm{hs}{y}^{\ast 2}}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{12}^1& =\frac{-\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^1& =\frac{\begin{array}{c}-h{(x^{\ast }+y^{\ast })}^3-h^2{(x^{\ast }+y^{\ast })}^3-h^{3}s{x}^{\ast 3}{y}^{\ast }\\ +\mathrm{hs}{y}^{\ast 2}+3h^{2}s{x}^{\ast }{y}^{\ast 2}+h^{2}s{y}^{\ast 3}\end{array}}{{(1+h{x}^{\ast })}^3{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^1& =\frac{h^{2}s{x}^2(x^{\ast }-y^{\ast })-2\mathrm{hs}{x}^{\ast }{y}^{\ast }}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{15}^1& =\frac{\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{11}^2& =\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},{\mathrm{\Omega }}_{12}^2=1-\mathrm{hr}{\delta }^{\ast }+\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^2& =-\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^2& =\frac{2h{\delta }^{\ast }{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},{\mathrm{\Omega }}_{15}^2=-\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3}.\end{array}$(A1) ) and (A2(A2) $\begin{array}{rl}\eta & =\frac{1}{2}(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr\delta }(r-1)),\\ \zeta & =\frac{1}{2}\sqrt{\begin{array}{c}4(\left({\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}\right)(1+\mathrm{hr\delta }(r-1))+{\displaystyle \frac{h^2{r}^2\mathrm{s\delta }{(r-1)}^2}{1+h+\mathrm{hs}(r-1)}})-{({\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}+1+\mathrm{hr\delta }(r-1))}^2\end{array}},\\ r_{11}& ={\mathrm{\Omega }}_{13}^1{\mathrm{\Omega }}_{12}^1+\eta {\mathrm{\Omega }}_{14}^1-{\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{14}^1-{\eta }^2\frac{{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}+\frac{{\left({\mathrm{\Omega }}_{11}^1\right)}^2{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}-\frac{2\eta {\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{12}& =-\zeta {\mathrm{\Omega }}_{14}^1-\frac{2\mathrm{\zeta \eta }{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}+\frac{2\zeta {\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{13}& ={\zeta }^2\frac{{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{21}& =\frac{1}{\zeta }((\eta -{\mathrm{\Omega }}_{11}^1){\mathrm{\Omega }}_{13}^1{\mathrm{\Omega }}_{12}^1+{(\eta -{\mathrm{\Omega }}_{11}^1)}^2{\mathrm{\Omega }}_{14}^1\frac{{(\eta -{\mathrm{\Omega }}_{11}^1)}^3{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}-{\mathrm{\Omega }}_{13}^2{\left({\mathrm{\Omega }}_{12}^1\right)}^2{\mathrm{\Omega }}_{14}^2{\mathrm{\Omega }}_{12}^1(\eta -{\mathrm{\Omega }}_{11}^1)-{\mathrm{\Omega }}_{15}^2{(\eta -{\mathrm{\Omega }}_{11}^1)}^2),\\ r_{22}& =-(\eta -{\mathrm{\Omega }}_{11}^1){\mathrm{\Omega }}_{14}^1-\frac{2(\eta -{\mathrm{\Omega }}_{11}^1)}{{\mathrm{\Omega }}_{12}^1}{\mathrm{\Omega }}_{15}^1+{\mathrm{\Omega }}_{14}^2{\mathrm{\Omega }}_{12}^1\\ & +2{\mathrm{\Omega }}_{15}^2(\eta -{\mathrm{\Omega }}_{11}^1),\\ r_{23}& =\zeta \frac{{\mathrm{\Omega }}_{15}^1(\eta -{\mathrm{\Omega }}_{11}^1)}{{\mathrm{\Omega }}_{12}^1}-\zeta {\mathrm{\Omega }}_{15}^2.\end{array}$(A2) ) one gets: (78) $\begin{array}{r}\{\begin{array}{l}{\mathrm{\Omega }}_{11}^1=1.0138772077375944,\\ {\mathrm{\Omega }}_{12}^1=-0.17283431455004208,\\ {\mathrm{\Omega }}_{13}^1=-0.2372883502540626,\\ {\mathrm{\Omega }}_{14}^1=-0.46146372048017426,\\ {\mathrm{\Omega }}_{15}^1=0.2743401818254637,\\ {\mathrm{\Omega }}_{11}^2=0.016500000000000046,\\ {\mathrm{\Omega }}_{12}^2=0.9834999999999999,\\ {\mathrm{\Omega }}_{13}^2=-0.02619047619047627,\\ {\mathrm{\Omega }}_{14}^2=0.05238095238095254,\\ {\mathrm{\Omega }}_{15}^2=-0.02619047619047627,\end{array}\end{array}$(78) and (79) $\begin{array}{r}\{\begin{array}{l}\eta =0.9986886038687973,\\ \zeta =\mathrm{0.05119641103234052.}\end{array}\end{array}$(79) Utilizing (78(78) $\begin{array}{r}\{\begin{array}{l}{\mathrm{\Omega }}_{11}^1=1.0138772077375944,\\ {\mathrm{\Omega }}_{12}^1=-0.17283431455004208,\\ {\mathrm{\Omega }}_{13}^1=-0.2372883502540626,\\ {\mathrm{\Omega }}_{14}^1=-0.46146372048017426,\\ {\mathrm{\Omega }}_{15}^1=0.2743401818254637,\\ {\mathrm{\Omega }}_{11}^2=0.016500000000000046,\\ {\mathrm{\Omega }}_{12}^2=0.9834999999999999,\\ {\mathrm{\Omega }}_{13}^2=-0.02619047619047627,\\ {\mathrm{\Omega }}_{14}^2=0.05238095238095254,\\ {\mathrm{\Omega }}_{15}^2=-0.02619047619047627,\end{array}\end{array}$(78) ) and (79(79) $\begin{array}{r}\{\begin{array}{l}\eta =0.9986886038687973,\\ \zeta =\mathrm{0.05119641103234052.}\end{array}\end{array}$(79) ) into (A2(A2) $\begin{array}{rl}\eta & =\frac{1}{2}(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr\delta }(r-1)),\\ \zeta & =\frac{1}{2}\sqrt{\begin{array}{c}4(\left({\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}\right)(1+\mathrm{hr\delta }(r-1))+{\displaystyle \frac{h^2{r}^2\mathrm{s\delta }{(r-1)}^2}{1+h+\mathrm{hs}(r-1)}})-{({\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}+1+\mathrm{hr\delta }(r-1))}^2\end{array}},\\ r_{11}& ={\mathrm{\Omega }}_{13}^1{\mathrm{\Omega }}_{12}^1+\eta {\mathrm{\Omega }}_{14}^1-{\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{14}^1-{\eta }^2\frac{{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}+\frac{{\left({\mathrm{\Omega }}_{11}^1\right)}^2{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}-\frac{2\eta {\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{12}& =-\zeta {\mathrm{\Omega }}_{14}^1-\frac{2\mathrm{\zeta \eta }{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}+\frac{2\zeta {\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{13}& ={\zeta }^2\frac{{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{21}& =\frac{1}{\zeta }((\eta -{\mathrm{\Omega }}_{11}^1){\mathrm{\Omega }}_{13}^1{\mathrm{\Omega }}_{12}^1+{(\eta -{\mathrm{\Omega }}_{11}^1)}^2{\mathrm{\Omega }}_{14}^1\frac{{(\eta -{\mathrm{\Omega }}_{11}^1)}^3{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}-{\mathrm{\Omega }}_{13}^2{\left({\mathrm{\Omega }}_{12}^1\right)}^2{\mathrm{\Omega }}_{14}^2{\mathrm{\Omega }}_{12}^1(\eta -{\mathrm{\Omega }}_{11}^1)-{\mathrm{\Omega }}_{15}^2{(\eta -{\mathrm{\Omega }}_{11}^1)}^2),\\ r_{22}& =-(\eta -{\mathrm{\Omega }}_{11}^1){\mathrm{\Omega }}_{14}^1-\frac{2(\eta -{\mathrm{\Omega }}_{11}^1)}{{\mathrm{\Omega }}_{12}^1}{\mathrm{\Omega }}_{15}^1+{\mathrm{\Omega }}_{14}^2{\mathrm{\Omega }}_{12}^1\\ & +2{\mathrm{\Omega }}_{15}^2(\eta -{\mathrm{\Omega }}_{11}^1),\\ r_{23}& =\zeta \frac{{\mathrm{\Omega }}_{15}^1(\eta -{\mathrm{\Omega }}_{11}^1)}{{\mathrm{\Omega }}_{12}^1}-\zeta {\mathrm{\Omega }}_{15}^2.\end{array}$(A2) ), one gets: (80) $\begin{array}{r}\{\begin{array}{l}{r}_{11}=3.2139366880763243,\\ r_{12}=0.021156708511249955,\\ r_{13}=-0.0041604325437973955,\\ r_{21}=-0.0014241879052949461,\\ r_{22}=-0.0634844121741236,\\ r_{23}=\mathrm{0.0025751472836634173.}\end{array}\end{array}$(80) Now in view of (80(80) $\begin{array}{r}\{\begin{array}{l}{r}_{11}=3.2139366880763243,\\ r_{12}=0.021156708511249955,\\ r_{13}=-0.0041604325437973955,\\ r_{21}=-0.0014241879052949461,\\ r_{22}=-0.0634844121741236,\\ r_{23}=\mathrm{0.0025751472836634173.}\end{array}\end{array}$(80) ) and (47(47) $\begin{array}{rl}{\varrho }_{02}& =\frac{1}{4}(r_{11}-r_{13}+r_{22}+\iota (r_{21}-r_{23}+r_{12})),\\ {\varrho }_{11}& =\frac{1}{2}(r_{11}+r_{13}+\iota (r_{21}+r_{23})),\\ {\varrho }_{20}& =\frac{1}{4}(r_{11}-r_{13}+r_{22}+\iota (r_{21}-r_{23}-r_{12})),\\ {\varrho }_{21}& =0.\end{array}$(47) ), one obtains (81) $\begin{array}{r}\{\begin{array}{l}{\varrho }_{02}=0.7886531771114995+0.004289343330572898\iota ,\\ {\varrho }_{11}=1.6048881277662634+0.0005754796891842356\iota ,\\ {\varrho }_{20}=0.7886531771114995-0.0062890109250520795\iota ,\\ {\varrho }_{21}=0.\end{array}\end{array}$(81) Now finally, using (81(81) $\begin{array}{r}\{\begin{array}{l}{\varrho }_{02}=0.7886531771114995+0.004289343330572898\iota ,\\ {\varrho }_{11}=1.6048881277662634+0.0005754796891842356\iota ,\\ {\varrho }_{20}=0.7886531771114995-0.0062890109250520795\iota ,\\ {\varrho }_{21}=0.\end{array}\end{array}$(81) ) along with $\lambda ,\overline{\lambda }=0.9986886038687973\pm 0.05119641103234052\iota$ into (45(45) $\begin{array}{rl}\ell & :=-\mathrm{\Re }\left(\frac{(1-2\overline{\lambda }){\overline{\lambda }}^2}{1-\lambda }{\varrho }_{11}{\varrho }_{20}\right)-\frac{1}{2}\Vert {\varrho }_{11}{\Vert }^2\\ & -\Vert {\varrho }_{02}{\Vert }^2+\mathrm{\Re }\left(\overline{\lambda }{\varrho }_{21}\right),\end{array}$(45) ) one gets: $\ell =3.9314493377802>0$ which give the fact that closed invariant curve must exists, and hence discrete model undergoes subcritical Neimark–Sacker at indicated interior fixed point ${\phi }_2=(0.31499999999999995,0.31499999999999995)$ (see Figure 3(a)). In similar manner, we can also obtain that if one varies $\delta =0.125670,0.130000,0.139999,0.141100,0.149999,0.158896,0.161231,0.167890>0.11000000000000032$ then stable invariant curves also appears which are depicted in Figure 3(b–i). Therefore, in conclusion we can say that numerical simulations in Example 7.1 agree with theoretical results obtained in Theorem 2.3, (i) of Theorem 2.5 and Theorem 4.1.

Figure 1. Bifurcation diagrams (B.Ds) with MLE of discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ). (a) B.D for $x_t$. (b) B.D for $y_t$. (c) B.D for $x_t$ and $y_t$. (d) B.D for s and $x_t$. (e) B.D for s and $y_t$. (f) MLEs.

Figure 2. Phase portrait of discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) with $(0.03,0.04)$. (a) If $\delta =0.10000$. (b) If $\delta =0.10299$. (c) If $\delta =0.10500$. (d) If $\delta =0.10880$. (e) If $\delta =0.10987$. (f) If $\delta =0.10999$.

Figure 3. Invariant closed curves of discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) with $(0.03,0.04)$. (a) If $\delta =0.110999$. (b) If $\delta =0.125670$. (c) If $\delta =0.130000$. (d) If $\delta =0.139999$. (e) If $\delta =0.141100$. (f) If $\delta =0.149999$. (g) If $\delta =0.158896$. (h) If $\delta =0.161231$. (i) If $\delta =0.167890$.

Example 7.2

If h = 1.28, s = 2.5, r = 0.798, $\delta =[10.1,16.99]$ with $(x_0,y_0)=(0.495,0.125301)$ then at $\delta =11.591710484782386$ discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) undergoes the flip bifurcation. The flip bifurcation diagrams and Maximum Lyapunov exponent are drawn in Figure 4. Further at $(h,s,r,\delta )=(1.28,2.5,0.798,11.591710484782386)$ discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) has equilibrium solution ${\phi }_2=(0.4950000000000001,0.12530075187969925)$ and moreover from (24(24) $\begin{array}{r}V{|}_{{\phi }_2}:=\left(\begin{array}{ccc}{\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}& {\displaystyle \frac{-h{r}^{2}s}{1+h+\mathrm{hs}(r-1)}\mathrm{h\delta }(r-1{)}^2}& 1+\mathrm{hr\delta }(r-1)\end{array}\right),\end{array}$(24) ) one gets: (82) $\begin{array}{rl}V{|}_{{\phi }_2}& =(\begin{array}{c}0.9279059745347696\\ 0.6054248379149572\end{array}\\ & \begin{array}{c}-1.2474123408423115\\ -1.3917278250303777\end{array}),\end{array}$(82) with ${\lambda }_1=-1$ and ${\lambda }_2=0.5361781495043931\ne 1\mathrm{or}-1$, and hence based on these simulations one can obtain that for parametric values $(h,s,r,\delta )=(1.28,2.5,0.798,11.591710484782386)\in \mathcal{F}{|}_{{\phi }_2=(0.4950000000000001,0.12530075187969925)}$. Moreover in this parametric domain, from (A3(A3) $\begin{array}{rl}{\mathrm{\Omega }}_{11}^1& =\frac{(1+h){(x^{\ast }+y^{\ast })}^2+h^{2}s{x}^{\ast 2}{y}^{\ast }-\mathrm{hs}{y}^{\ast 2}}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{12}^1& =\frac{-\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^1& =\frac{\begin{array}{c}-h{(x^{\ast }+y^{\ast })}^3-h^2{(x^{\ast }+y^{\ast })}^3-h^{3}s{x}^{\ast 3}{y}^{\ast }+\mathrm{hs}{y}^{\ast 2}\\ +3h^{2}s{x}^{\ast }{y}^{\ast 2}+h^{2}s{y}^{\ast 3}\end{array}}{{(1+h{x}^{\ast })}^3{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^1& =\frac{h^{2}s{x}^2(x^{\ast }-y^{\ast })-2\mathrm{hs}{x}^{\ast }{y}^{\ast }}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{15}^1& =\frac{\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{11}^2& =\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},{\mathrm{\Omega }}_{12}^2=1-\mathrm{hr}{\delta }^{\ast }+\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^2& =-\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^2& =\frac{2h{\delta }^{\ast }{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},{\mathrm{\Omega }}_{15}^2=-\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},k_1=\frac{h{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ k_2& =-\mathrm{hr}+\frac{h{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},k_3=\frac{-h{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ k_4& =\frac{2h{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},k_5=\frac{-h{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3}.\end{array}$(A3) ), (54(54) $\begin{array}{r}\{\begin{array}{l}{v}_1=\left\{\begin{array}{rl}& {\displaystyle \frac{1}{1-{\lambda }_2}(\frac{\begin{array}{c}-(2+h(1+\mathrm{sr}-s))\\ (-2+h(-1+s(r-1{)}^2))\end{array}}{\begin{array}{c}4+h(1+s(r-1))\\ (4+h+hs(r^2-1))\end{array}}}\\ & {\displaystyle \times ({\mathrm{\Omega }}_{13}^1+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}{\mathrm{\Omega }}_{14}^1}\\ & +{\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{\mathrm{\Omega }}_{15}^1)\\ & {\displaystyle -\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{\begin{array}{c}4+h(1+s(r-1))\\ (4+h+\mathrm{hs}(r^2-1))\end{array}}}\\ & {\displaystyle \times ({\mathrm{\Omega }}_{13}^2+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}{\mathrm{\Omega }}_{14}^2}\\ & +{\left(\frac{2-\mathrm{hs}(-1+h{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{\mathrm{\Omega }}_{15}^2))\end{array}\right\},\\ v_2=\left\{\begin{array}{r}{\displaystyle \frac{1}{1-{\lambda }_2}}\\ (\frac{-h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{\begin{array}{c}4+h(1+s(r-1))\\ (4+h+\mathrm{hs}(r^2-1))\end{array}}\\ {\displaystyle \times (k_1+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}{k}_2))}\end{array}\right\},\\ v_0=v_3=0.\end{array}\end{array}$(54) ) and (56(56) $\begin{array}{r}\{\begin{array}{l}{m}_1=\left\{\begin{array}{rl}& {\displaystyle \frac{2\mathrm{hrs}(r-1)(1+h+\mathrm{hs}(r-1))}{4+h(1+s(r-1))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle \times ({\mathrm{\Omega }}_{13}^1+\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right){\mathrm{\Omega }}_{14}^1{\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{\mathrm{\Omega }}_{15}^1)}\\ & {\displaystyle +\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle \times ({\mathrm{\Omega }}_{13}^2+\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right){\mathrm{\Omega }}_{14}^2{\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{\mathrm{\Omega }}_{15}^2)}\end{array}\right\},\\ m_2={\displaystyle \frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}}(k_1+{\displaystyle \frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}}{k}_2),\\ m_3=\left\{\begin{array}{rl}& {\displaystyle \frac{2\mathrm{hrs}(r-1)(1+h+\mathrm{hs}(r-1))}{4+h(1+s(r-1))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle \times (2{\mathrm{\Omega }}_{13}^1{v}_2+{\mathrm{\Omega }}_{14}^1{v}_2(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2})}\\ & {\displaystyle +{\mathrm{\Omega }}_{15}^1{v}_2\left(2\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)\left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right)\right))}\\ & {\displaystyle +\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}(2{\mathrm{\Omega }}_{13}^2{v}_2{\displaystyle +{\mathrm{\Omega }}_{14}^2{v}_2(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}}}\\ \\ & +\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2})\\ & {\displaystyle +{\mathrm{\Omega }}_{15}^2{v}_2(2\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}\\ & \left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right))\\ & {\displaystyle +k_1{v}_1+\left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right)k_2{v}_1}\\ & {\displaystyle +k_3+\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)k_4}\\ & +{\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)}^2{k}_5)\end{array}\right\},\\ m_4=\left\{\begin{array}{rl}& {\displaystyle \frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle (k_1{v}_2+\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}{k}_2{v}_2)}\end{array}\right\},\\ m_5=\left\{\begin{array}{rl}& {\displaystyle \frac{2\mathrm{hrs}(r-1)(1+h+\mathrm{hs}(r-1))}{4+h(1+s(r-1))(4+h+\mathrm{hs}(r^2-1))}}\\ & {\displaystyle \times (2{\mathrm{\Omega }}_{13}^1{v}_1+{\mathrm{\Omega }}_{14}^1{v}_1(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2})}\\ & {\displaystyle +{\mathrm{\Omega }}_{15}^1{v}_1\left(2\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)\left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right))\right)}\\ & {\displaystyle +\frac{h{r}^{2}s(2+h(1+\mathrm{sr}-s))}{4+h(1+s(-1+r))(4+h+\mathrm{hs}(r^2-1))}(2{\mathrm{\Omega }}_{13}^2{v}_1\phantom{\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}}}\\ & {\displaystyle +{\mathrm{\Omega }}_{14}^2{v}_1(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}+\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2})}\\ & {\displaystyle +{\mathrm{\Omega }}_{15}^2{v}_1\left(2\left(\frac{2-\mathrm{hs}(-1+r{)}^2+h}{\mathrm{hs}{r}^2}\right)\left(\frac{2(1-r)(1+h(1+\mathrm{sr}-s))}{r(2+h-\mathrm{hs}+\mathrm{hsr})}\right))\right)}\end{array}\right\}.\end{array}\end{array}$(56) ) one gets: (83) $\begin{array}{rl}& \{\begin{array}{l}{\mathrm{\Omega }}_{11}^1=0.9279059745347696,\\ {\mathrm{\Omega }}_{12}^1=-1.2474123408423115,\\ {\mathrm{\Omega }}_{13}^1=-0.5982005795778922,\\ {\mathrm{\Omega }}_{14}^1=-0.040685691442227376,\\ {\mathrm{\Omega }}_{15}^1=2.0109798949336652,\\ {\mathrm{\Omega }}_{11}^2=0.6054248379149572,\\ {\mathrm{\Omega }}_{12}^2=-1.3917278250303777,\\ {\mathrm{\Omega }}_{13}^2=-0.976018223547749,\\ {\mathrm{\Omega }}_{14}^2=7.711510320703999,\\ {\mathrm{\Omega }}_{15}^2=-15.232141673073743,\\ k_1=0.05222911999999998,\\ k_2=-0.20633088,\\ k_3=-0.0841996722424242,\\ k_4=0.6652607767272725,\\ k_5=-1.3140547025454543,\end{array}\end{array}$(83) (84) $\begin{array}{rl}& \{\begin{array}{l}{v}_0=0,\\ v_1=12.929077600120031,\\ v_2=0.10826715884189564,\\ v_3=0,\end{array}\end{array}$(84) and (85) $\begin{array}{r}\{\begin{array}{l}{m}_1=-21.71571925204092,\\ m_2=-0.2165343176837912,\\ m_3=-2.143746543978618,\\ m_4=-0.0011046972084443997,\\ m_5=-\mathrm{27.41483556845985.}\end{array}\end{array}$(85) Utilizing (85(85) $\begin{array}{r}\{\begin{array}{l}{m}_1=-21.71571925204092,\\ m_2=-0.2165343176837912,\\ m_3=-2.143746543978618,\\ m_4=-0.0011046972084443997,\\ m_5=-\mathrm{27.41483556845985.}\end{array}\end{array}$(85) ) into (57(57) $\begin{array}{rl}{\ell }_1& :=(\frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_t\mathrm{\partial }ϵ}+\frac{1}{2}\frac{\mathrm{\partial }{f}_1}{\mathrm{\partial }ϵ}\frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_t^2}){|}_{(0,0)}=m_2,\\ {\ell }_2& :=(\frac{1}{6}\frac{{\mathrm{\partial }}^3{f}_1}{\mathrm{\partial }{x}_t^3}+{\left(\frac{1}{2}\frac{{\mathrm{\partial }}^2{f}_1}{\mathrm{\partial }{x}_t^2}\right)}^2){|}_{(0,0)}=m_5+m_1^2.\end{array}$(57) ) one gets: ${\ell }_1=-0.2165343176837912\ne 0$ and ${\ell }_2=444.1576270650008>0$. Since ${\ell }_2=444.1576270650008>0$ and so it can be concluded that stable period-2 points bifurcate from ${\phi }_2=(0.4950000000000001,0.12530075187969925)$. Therefore, in conclusion we can say that numerical simulations in Example 7.2 agree with theoretical results obtained in Theorem 2.4, (ii) of Theorem 2.5 and Theorem 5.1.

Figure 4. Flip B.Ds with MLE of discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ). (a) B.D for $x_t$. (b) B.D for $y_t$. (c) B.D for $x_t$ and $y_t$. (d) B.D for h and $x_t$. (e) B.D for h and $y_t$. (f) B.D for s and $x_t$. (g) B.D for s and $y_t$. (h) MLEs.

Example 7.3

If h = 0.5, r = 0.4, $\delta =0.38333333333333336$, $s\in [0.1,1.33]$ with $(x_0,y_0)=(0.02,0.03)$ then model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) undergoes Neimark–Sacker bifurcation. The maximum Lyapunov exponents with bifurcation diagrams are drawn in Figure 5. Further if s = 1.3 then model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) has equilibrium ${\phi }_2=(0.21999999999999997,0.32999999999999996)$ where roots of characteristics equations of $V{|}_{{\phi }_2=(0.21999999999999997,0.32999999999999996)}$ are ${\lambda }_{1,2}=0.9977207207207206\pm 0.06747861471996602\iota$ with $|{\lambda }_{1,2}|=1$. In order to apply OGY control method for model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ), if $s_0=0.3$ then it has equilibrium ${\phi }_2=(0.82,1.2299999999999998)$. So, model (64(64) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-{\displaystyle \frac{h(s_0-k_1(x_t-x)-k_2(y_t-y))x_t{y}_t}{(x_t+y_t)(1+h{x}_t)}},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}).}}\end{array}\end{array}$(64) ) becomes (86) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=x_t+0.5x_t(1-x_t)-{\displaystyle \frac{\begin{array}{c}0.5(0.3-k_1(x_t-0.82)\\ -k_2(y_t-1.2299999999999998))x_t{y}_t\end{array}}{x_t+y_t}},}\\ {\displaystyle y_{t+1}=y_t+0.5(0.38333333333333336)y_t(-0.4+{\displaystyle \frac{x_t}{x_t+y_t}}),}\end{array}\end{array}$(86) where its variation matrix V −BK is (87) $\begin{array}{rl}V-\mathrm{BK}& =(\begin{array}{c}0.7347517730496455\\ +0.17446808510638295k_1\\ 0.0689994\end{array}\\ & \begin{array}{c}-0.017021276595744685\\ +0.17446808510638295k_1\\ 0.9540004\end{array}).\end{array}$(87) Now for marginal stability lines $L_1,L_2$ and $L_3$ are (88) $\begin{array}{rl}& L_1:3.340151253503999k_1+16.523026880255998k_2\\ & +13.10398752=0,\end{array}$(88) (89) $\begin{array}{rl}& L_2:11.015351253503999k_1+16.523026880255998k_2\\ & -18.358918755839998=0,\end{array}$(89) and (90) $\begin{array}{rl}& L_3:-4.335048746496k_1+16.523026880255998k_2\\ & -10.60110620416=0.\end{array}$(90) So, lines (88(88) $\begin{array}{rl}& L_1:3.340151253503999k_1+16.523026880255998k_2\\ & +13.10398752=0,\end{array}$(88) ), (89(89) $\begin{array}{rl}& L_2:11.015351253503999k_1+16.523026880255998k_2\\ & -18.358918755839998=0,\end{array}$(89) ) and (90(90) $\begin{array}{rl}& L_3:-4.335048746496k_1+16.523026880255998k_2\\ & -10.60110620416=0.\end{array}$(90) ) determine triangular region that gives $|{\lambda }_{1,2}|<1$(see Figure 6).

Figure 5. B.Ds with MLE of discrete model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ). (a) B.D for $x_t$. (b) B.D for $y_t$. (c) B.D for $x_t$ and $y_t$. (d) MLEs.

Figure 6. Region of stability where $|{\lambda }_{1,2}|<1$.

Example 7.4

Finally, if $h=1.28,s=2.5,r=0.798,\delta =11.591710484782386$ with $(x_0,y_0)=(1.3,1.12)$ then model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) undergoes flip bifurcation. For this, by applying hybrid strategy to get stable orbit at ${\phi }_2=(0.4950000000000001,0.12530075187969925)$. For this model (72(72) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\alpha (\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)})+(1-\alpha )x_t,}\\ {\displaystyle y_{t+1}=\alpha (y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}))+(1-\alpha )y_t,}\end{array}\end{array}$(72) ) takes the form (91) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\alpha (x_t+1.28x_t(1-x_t)-\frac{(1.28)(2.5)x_t{y}_t}{x_t+y_t})+(1-\alpha )x_t,}\\ {\displaystyle y_{t+1}=\alpha (y_t+(1.28)(11.591710484782386)y_t(-0.798+\frac{x_t}{x_t+y_t}))+(1-\alpha )y_t,}\end{array}\end{array}$(91) where $V_{{\phi }_2=(0.4950000000000001,0.12530075187969925)}$ is (92) $\begin{array}{rl}{V}_{{\phi }_2}& =(\begin{array}{c}1-0.07209402546523036\alpha \\ 0.6054248379149572\alpha \end{array}\\ & \begin{array}{c}-1.2474123408423115\alpha \\ 1-2.3917278250303777\alpha \end{array}),\end{array}$(92) with (93) $\begin{array}{rl}& {\lambda }^2-(2-2.463821850495606\alpha )\lambda +1\\ & -2.463821850495606\alpha \\ & +0.9276437009912135{\alpha }^2=0.\end{array}$(93) Furthermore, roots of (93(93) $\begin{array}{rl}& {\lambda }^2-(2-2.463821850495606\alpha )\lambda +1\\ & -2.463821850495606\alpha \\ & +0.9276437009912135{\alpha }^2=0.\end{array}$(93) ) satisfying $|{\lambda }_{1,2}|<1$ if $0<\alpha <1$. So, for the allowed interval of control parameter α the flip bifurcation is completely eliminated. If $\alpha =0.999$ then for controlled model (91(91) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\alpha (x_t+1.28x_t(1-x_t)-\frac{(1.28)(2.5)x_t{y}_t}{x_t+y_t})+(1-\alpha )x_t,}\\ {\displaystyle y_{t+1}=\alpha (y_t+(1.28)(11.591710484782386)y_t(-0.798+\frac{x_t}{x_t+y_t}))+(1-\alpha )y_t,}\end{array}\end{array}$(91) ), plots of t vs. $x_t$ and $y_t$ are drawn in Figure 7.

Figure 7. Plot of t vs. $x_t$ and $y_t$ for controlled system (58(58) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)}=f_1(x_t,y_t,s),{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t})=f_2(x_t,y_t,s),}}\end{array}\end{array}$(58) ).

8. Conclusion

This work explored the local dynamics at fixed points, chaos, bifurcations of a discrete prey–predator model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) in the region: ${\mathbb{R}}_{+}^2=\{(x,y):x,y\ge 0\}$. We studied local dynamics at semitrivial and interior fixed points ${\phi }_{1,2}$ of model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ), and explored that ${\phi }_1$ of model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) is a sink if $h>\frac{2}{\delta (r-1)}$, never source, saddle if $0<h<\frac{2}{\delta (r-1)}$, and non-hyperbolic if $h=\frac{2}{\delta (r-1)}$. Further it proved that ${\phi }_2$ of model (13(13) $\begin{array}{r}\{\begin{array}{l}{\displaystyle x_{t+1}=\frac{x_t(1+h)}{1+h{x}_t}-\frac{\mathrm{hs}{x}_t{y}_t}{(x_t+y_t)(1+h{x}_t)},{\displaystyle y_{t+1}=y_t+\mathrm{h\delta }{y}_t(-r+\frac{x_t}{x_t+y_t}),}}\end{array}\end{array}$(13) ) is a stable focus if $0<\delta <\frac{1-s+s{r}^2}{r^2-r}$ with $r>max\{\pm \sqrt{\frac{s-1}{s}},1\}$, an unstable focus if $\delta >\frac{1-s+s{r}^2}{r^2-r}$ with $r>max\{\pm \sqrt{\frac{s-1}{s}},1\}$, non-hyperbolic condition if $\delta =\frac{1-s+s{r}^2}{r^2-r}$, stable node if $\delta >\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$ with $s>max\{\frac{2+h}{h(r-1{)}^2},\frac{2+h}{h(1-r)}\}$, an unstable node if $0<\delta <\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$ with $s>max\{\frac{2+h}{h(r-1{)}^2},\frac{2+h}{h(1-r)}\}$, and non-hyperbolic if $\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}$. In order to studied the bifurcation analysis at ${\phi }_{1,2}$, we first identified bifurcation sets for understudied model (i) flip bifurcation set $\mathcal{F}{|}_{{\phi }_1}:=\{(h,s,\delta ,r),h=\frac{2}{\delta (r-1)}\}$ at ${\phi }_1$, (ii) Neimark–Sacker bifurcation set $\mathcal{N}{|}_{{\phi }_2}:=\{(h,s,\delta ,r),\delta =\frac{1-s+s{r}^2}{r^2-r}\}$ at ${\phi }_2$ (iii) flip bifurcation set whose mathematical expression is $\mathcal{F}{|}_{{\phi }_2}:=\{(h,s,\delta ,r),\delta =\frac{-4-2h(1-s(r-1{)}^2)}{\mathrm{hr}(r-1)(2+h+\mathrm{hs}(r-1))}\}$ at ${\phi }_2$, and then proved that at ${\phi }_1$ flip bifurcation did not take place if $(h,s,\delta ,r)\in \mathcal{F}{|}_{{\phi }_1}$ but at ${\phi }_2$ model subjected Neimark–Sacker and flip bifurcations if $(h,s,\delta ,r)\in \mathcal{N}{|}_{{\phi }_2}$ and $(h,s,\delta ,r)\in \mathcal{F}{|}_{{\phi }_2}$, respectively. Furthermore, OGY and Hybrid control strategies are utilized to control chaos in the under study discrete model due to occurrence of Neimark–Sacker and flip bifurcations. Finally numerical simulations depicted to verify theoretical results.

Acknowledgments

Authors declare that they got no funding on any part of this research.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Appendices

Appendix 1.

Expression for the coefficients quantities of system (41)

(A1) $\begin{array}{rl}{\mathrm{\Omega }}_{11}^1& =\frac{(1+h){(x^{\ast }+y^{\ast })}^2+h^{2}s{x}^{\ast 2}{y}^{\ast }-\mathrm{hs}{y}^{\ast 2}}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{12}^1& =\frac{-\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^1& =\frac{\begin{array}{c}-h{(x^{\ast }+y^{\ast })}^3-h^2{(x^{\ast }+y^{\ast })}^3-h^{3}s{x}^{\ast 3}{y}^{\ast }\\ +\mathrm{hs}{y}^{\ast 2}+3h^{2}s{x}^{\ast }{y}^{\ast 2}+h^{2}s{y}^{\ast 3}\end{array}}{{(1+h{x}^{\ast })}^3{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^1& =\frac{h^{2}s{x}^2(x^{\ast }-y^{\ast })-2\mathrm{hs}{x}^{\ast }{y}^{\ast }}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{15}^1& =\frac{\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{11}^2& =\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},{\mathrm{\Omega }}_{12}^2=1-\mathrm{hr}{\delta }^{\ast }+\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^2& =-\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^2& =\frac{2h{\delta }^{\ast }{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},{\mathrm{\Omega }}_{15}^2=-\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3}.\end{array}$(A1)

Appendix 2.

Expression for the coefficients quantities of system (44)

(A2) $\begin{array}{rl}\eta & =\frac{1}{2}(\frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}+1+\mathrm{hr\delta }(r-1)),\\ \zeta & =\frac{1}{2}\sqrt{\begin{array}{c}4(\left({\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}\right)(1+\mathrm{hr\delta }(r-1))+{\displaystyle \frac{h^2{r}^2\mathrm{s\delta }{(r-1)}^2}{1+h+\mathrm{hs}(r-1)}})-{({\displaystyle \frac{1-\mathrm{hrs}(r-1)}{1+h+\mathrm{hs}(r-1)}}+1+\mathrm{hr\delta }(r-1))}^2\end{array}},\\ r_{11}& ={\mathrm{\Omega }}_{13}^1{\mathrm{\Omega }}_{12}^1+\eta {\mathrm{\Omega }}_{14}^1-{\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{14}^1-{\eta }^2\frac{{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}+\frac{{\left({\mathrm{\Omega }}_{11}^1\right)}^2{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}-\frac{2\eta {\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{12}& =-\zeta {\mathrm{\Omega }}_{14}^1-\frac{2\mathrm{\zeta \eta }{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}+\frac{2\zeta {\mathrm{\Omega }}_{11}^1{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{13}& ={\zeta }^2\frac{{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1},\\ r_{21}& =\frac{1}{\zeta }((\eta -{\mathrm{\Omega }}_{11}^1){\mathrm{\Omega }}_{13}^1{\mathrm{\Omega }}_{12}^1+{(\eta -{\mathrm{\Omega }}_{11}^1)}^2{\mathrm{\Omega }}_{14}^1\frac{{(\eta -{\mathrm{\Omega }}_{11}^1)}^3{\mathrm{\Omega }}_{15}^1}{{\mathrm{\Omega }}_{12}^1}-{\mathrm{\Omega }}_{13}^2{\left({\mathrm{\Omega }}_{12}^1\right)}^2{\mathrm{\Omega }}_{14}^2{\mathrm{\Omega }}_{12}^1(\eta -{\mathrm{\Omega }}_{11}^1)-{\mathrm{\Omega }}_{15}^2{(\eta -{\mathrm{\Omega }}_{11}^1)}^2),\\ r_{22}& =-(\eta -{\mathrm{\Omega }}_{11}^1){\mathrm{\Omega }}_{14}^1-\frac{2(\eta -{\mathrm{\Omega }}_{11}^1)}{{\mathrm{\Omega }}_{12}^1}{\mathrm{\Omega }}_{15}^1+{\mathrm{\Omega }}_{14}^2{\mathrm{\Omega }}_{12}^1\\ & +2{\mathrm{\Omega }}_{15}^2(\eta -{\mathrm{\Omega }}_{11}^1),\\ r_{23}& =\zeta \frac{{\mathrm{\Omega }}_{15}^1(\eta -{\mathrm{\Omega }}_{11}^1)}{{\mathrm{\Omega }}_{12}^1}-\zeta {\mathrm{\Omega }}_{15}^2.\end{array}$(A2)

Appendix 3.

Expression for the coefficients quantities of system (48)

(A3) $\begin{array}{rl}{\mathrm{\Omega }}_{11}^1& =\frac{(1+h){(x^{\ast }+y^{\ast })}^2+h^{2}s{x}^{\ast 2}{y}^{\ast }-\mathrm{hs}{y}^{\ast 2}}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{12}^1& =\frac{-\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^1& =\frac{\begin{array}{c}-h{(x^{\ast }+y^{\ast })}^3-h^2{(x^{\ast }+y^{\ast })}^3-h^{3}s{x}^{\ast 3}{y}^{\ast }+\mathrm{hs}{y}^{\ast 2}\\ +3h^{2}s{x}^{\ast }{y}^{\ast 2}+h^{2}s{y}^{\ast 3}\end{array}}{{(1+h{x}^{\ast })}^3{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^1& =\frac{h^{2}s{x}^2(x^{\ast }-y^{\ast })-2\mathrm{hs}{x}^{\ast }{y}^{\ast }}{{(1+h{x}^{\ast })}^2{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{15}^1& =\frac{\mathrm{hs}{x}^{\ast 2}}{(1+h{x}^{\ast }){(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{11}^2& =\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},{\mathrm{\Omega }}_{12}^2=1-\mathrm{hr}{\delta }^{\ast }+\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ {\mathrm{\Omega }}_{13}^2& =-\frac{h{\delta }^{\ast }{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ {\mathrm{\Omega }}_{14}^2& =\frac{2h{\delta }^{\ast }{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},{\mathrm{\Omega }}_{15}^2=-\frac{h{\delta }^{\ast }{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},k_1=\frac{h{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},\\ k_2& =-\mathrm{hr}+\frac{h{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^2},k_3=\frac{-h{y}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3},\\ k_4& =\frac{2h{x}^{\ast }{y}^{\ast }}{{(x^{\ast }+y^{\ast })}^3},k_5=\frac{-h{x}^{\ast 2}}{{(x^{\ast }+y^{\ast })}^3}.\end{array}$(A3)