A Dynamically Consistent Nonstandard Difference Scheme for a Discrete-Time Immunogenic Tumors Model

Abstract

This manuscript deals with the qualitative study of certain properties of an immunogenic tumors model. Mainly, we obtain a dynamically consistent discrete-time immunogenic tumors model using a nonstandard difference scheme. The existence of fixed points and their stability are discussed. It is shown that a continuous system experiences Hopf bifurcation at one and only one positive fixed point, whereas its discrete-time counterpart experiences Neimark–Sacker bifurcation at one and only one positive fixed point. It is shown that there is no chance of period-doubling bifurcation in our discrete-time system. Additionally, numerical simulations are carried out in support of our theoretical discussion.

1. Introduction

A tumor is a cluster of infections produced by the irregular evolution of cells described by DNA latent to a blowout in additional body parts. In every category of a tumor, certain cells in the body develop strangely and damage the nearby tissues. The healthy physique has a strict immune arrangement to protect it in case of tumors caused by the letdown of the immune system and the supplementary mechanism inside the body [1]. Tumor cases have been increasing very rapidly in recent eras, and it has been estimated that approximately 18.1 million persons suffer from them each year, out of which 9.6 million die [2]. The incidence of new tumors has been ascending, and it is projected that the number may reach 22.2 million by 2030 [3]. Therefore, it is necessary to improve novel, progressive, and cost-effective approaches to deal with these circumstances. Presently, the most public tumor treatments are chemotherapy [4], immunotherapy ([5,6]), surgery [7], and radiotherapy [8]. Despite all of these management options, it reverts. Thus it is essential for additional and operational treatment to be understandable. The immune reaction to a tumor is frequently cell-refereed cytotoxic T lymphocytes (CTL), and natural killer (NK) cells show a leading character. Various scientific models of the interactions between the increasing tumor and immune system have been established [9,10,11,12]. Moreover, several mathematical models describe the kinetics of cells refereed cytotoxicity in vitro [13,14,15]. With these mathematical models, several occurrences understood, statistical estimates for biologically essential factors have been acquired, and guesses made. The qualitative analysis of anti-tumor immune response in vivo is complicated and not well-understood. Freely ascending cancers have little immunogenicity, and frequently spread out of control in most creatures. Sometimes, the escape of any tumor from immune reconnaissance is connected with many different applications, namely, the damage or covering of cancer antigens, cancer-influenced disorders in safe regulation, damage to MHC class-I particles, and the addition of a variety of cancer duplicates resistant to cytolytic mechanisms [16,17,18]. Although attacking cells of the immune system may kill tumor cells, protected reconnaissance of natural cancers may be operative and significant in keeping tumor frequency low. The leading efforts to improve arrangements for immunotherapy or its grouping with other treatment approaches focus on reducing cancer mass. However, the bulk of such efforts remain ineffective. The critical explanation for this is that even after “effective” and “clinically comprehensive” elimination of cancer cells, a small number of “remaining” cancer cells remain, which may develop into subordinate cancers or “latent” metastases.

Cancer dormancy is a functioning span used to define a state in which a possibly dead cancer cells persevere for a protracted time with slight or no growth in the cancer cell population. It is commonly supposed that cancer cells do not develop at a speedy frequency throughout the dormancy period, apparently due to the nonappearance of a factor required for advanced evolution into cancer [19,20,21]. Nevertheless, a substitute probability is that quickly increasing cells are exterminated at a frequency equivalent to their creation. Undeveloped conditions develop both during essential treatment for cancer, and in the initial phases of cancer development. In truth, there is a typical arrangement in that neoplastic cells discharge from significant cancer very initial in its growth in any person [22]. The providence of these evading neoplastic compartments regulates whether the enduring cancer survives or kills the tumor. The straight contribution of CTL in the provision of a latent cancer state has been revealed in certain specific investigational models. Moreover, CTL’s different kinds of protected system cells, such as NK cells and macrophages, may preserve a latent cancer state [22].

In the past two decades, many authors have presented various clarifications for the expiry of a latent cancer state, snitching over cancers, and immune inspiration properties. Frequently, these clarifications are centred on the concepts of protected range, antigenic variation, creation by cancer cells of unlike kinds of immune cell delaying factors, group of immune suppressor cells, variations in auto-governing systems in a cancer localization area, and other new complex concepts that are very challenging to verify or invalidate experimentally. We consider that these occurrences might result from nonlinear dynamic connections between cancer and the insusceptible system [23,24,25,26]. The authors of [26] have examined a simple scientific model of a compartment-refereed reaction to a developing cancer compartment population. They have explained that their model contrasts with most models in the literature because it explains the penetration of cancer by consequence or cells and the probability of consequence or cell in the beginning. Kuznetsov [26] have studied an alternative to this model. Moreover, ref. [26] discussed the qualitative conduct of the structure using methods from the bifurcation concept. They applied the model to discuss the appliances of cancer latency and snitching through. They discussed that a non-zero frequency of consequence or cell, in the beginning, is necessary to achieve snitching through. They have found that snitching over, cancer latency, and the immune motivation of cancer development, properties which have been investigated independently, might be connected, which is similar to our model. Here, we study cancers with cells that are “immune genetic”, and consequently focus on the insusceptible attack by cytotoxic effector cells, for example, CTL, or NK cells. The communication among effector cells (EC) and cancer cells (TC) in vitro can defined using the kinetic scheme below.

Where, in Figure 1, $T,E,C,T^{*},E^{*}$ are the limited meditations of cancer cells, effector cells, effector cell–cancer cell conjugates, “critically hit” TC cells and deactivated effector cells, respectively. Critically hit cancer cells are intended to die. They similarly have been named cells “encoded to die”. The addition of deactivated effector cells is a strange feature of this model. In cases with a lesser degree of CTL, culture appears to have a limited ability to constantly destroy target cells [27]. This is because particle collapse is responsible for cytotoxic influences or other controlling effects, probably due to the discharge of particles from the cancer cell after EC and TC are conjugated. Moreover, $k_1$, $k_{-1}$, $k_2$, and $k_3$ are non-negative kinetic real numbers; $k_1$ and $k_{-1}$ refer to the degrees of binding of TC from EC and objectivity of TC to EC without injuring cells, $k_2$ is the degree at which EC to TC connections conclusively program TC for lysis, and $k_3$ is the degree at which EC to TC connections deactivate EC. Hence, we have the following system of differential equations as a model for the communication among EC and increasing immunogenic cancer in vivo [27]:

$$\left\{\begin{array}{c}\frac{dE}{dt}=s+G(C,T)-d_{1}E-k_{1}ET+(k_{-1}+k_2)C,\hfill \\ \frac{dT}{dt}=aT(1-bT)-k_{1}ET+(k_{-1}+k_3)C,\hfill \\ \frac{dC}{dt}=k_{1}ET-(k_{-1}+k_2+k_3)C,\hfill \\ \frac{d{E}^{*}}{dt}=k_{3}C-d_2{E}^{*},\hfill \\ \frac{d{T}^{*}}{dt}=k_{2}C+d_3{T}^{*},\hfill \end{array}\right.$$

(1)

where

$$G(C,T)=\frac{fC}{g+T}.$$

The authors of [28] have explained that the last two equations from the system (1) are “slaves” to the first three equations in that system, as variables $T^{*}$ and $E^{*}$ have no influence on each other or the other variables in the system. Hence, in their work, they reduced the system (1) to the first three equations (see [28]) in order to dictate the dynamics of the system (1). Moreover, the development and division of cellular conjugates C follows the time measure of numerous tens of minutes to a limited number of hours. A time interval of this type is similarly detected before the disintegration of lethally hit tumor cells by separating the cell wall or membrane. However, the growth along with the inflow of effector cells into the spleen arises on a considerably slower time measure, possibly tens of hours. This inspires the application of a quasi-steady state estimate to the third equation of system (1) (that is, $\frac{dC}{dt}\approx 0$), which yields the following system of equations (see [28]):

$$\left\{\begin{array}{c}\frac{dE}{dt}=s+\frac{pET}{g+T}-mET-dE,\hfill \\ \frac{dT}{dt}=aT(1-bT)-nET,\hfill \end{array}\right.$$

(2)

where $p=fK,m=K{k}_3,n=K{k}_2,$ and $d=d_l$ are dimensional parameters. In addition, for better study of the dynamics of model (2) it is necessary to non-dimensionalize the system (2).

Non-Dimensionalization of System

The non-dimensionalized form of model (2) is obtained by selecting an order of degree application measure for E and T cell populations, $E_0$ and $T_0$, respectively, as proposed from the tests discussed in the previous section: $E_0$ = $T_0$ = ${10}^6$ cells (see [28]). Time is scaled comparative to the degree of cancer cell deactivation, that is, $\tau =n{T}_{0}t.$ Formally, the model can be re-articulated as

$$\left\{\begin{array}{c}\frac{dx}{dt}=\sigma +\frac{\rho x\left(t\right)y\left(t\right)}{\eta +y\left(t\right)}-\mu x\left(t\right)y\left(t\right)-\delta x\left(t\right),\hfill \\ \frac{dy}{dt}=\alpha y\left(t\right)(1-\beta y\left(t\right))-x\left(t\right)y\left(t\right),\hfill \end{array}\right.$$

(3)

where $x=\frac{E}{E_0},y=\frac{T}{T_0}$ and parameter estimates for system (3) are the following (see

Table 1. Parameter estimates from [28].

System Parameters	Dimensional Parameters	Estimated Values
$\sigma$	$\sigma =\frac{s}{n{E}_0{T}_0}$	$0.001181$
$\rho$	$\rho =\frac{p}{n{T}_0}$	$11.131$
$\eta$	$\eta =\frac{g}{T_0}$	$20.19$
$\mu$	$\mu =\frac{m}{n}$	$0.00311$
$\delta$	$\delta =\frac{d}{n{T}_0}$	$0.3743$
$\alpha$	$\alpha =\frac{a}{n{T}_0}$	$1.636$
$\beta$	$\beta =bn{T}_0$	$2.0\times {10}^{-3}$

):

The authors of [29] considered the post-collapsing conduct of functionally classified curved shell sections. Moreover, they investigated different shell geometries (cylindrical, elliptical, spherical, and hyperbolic) under the biaxial and uniaxial edge density. Duan et al. [30] have discussed a chemotherapy system’s stability analysis in cancer–immune responses. Moreover, their system has more than one fixed point. In order to conduct a stability analysis of these fixed points, they computed the upper Lyapunov exponents of the linearized system for these fixed points. They show that, while one fixed point is globally asymptotically stable if the noise is weak, another fixed point is constantly unstable whether the noise is weak or strong. We refer readers to [31] for further information on cancer–immune response systems. The authors of [32] discussed the discrete-time counterpart of a tumor–immune interaction model and analyzed the bifurcation in that model in fractional form. The authors of [33] studied the prime control for a tumor behaviour mathematical model using Atangana–Baleanu–Caputo fractional derivatives. In [34], the modeling and study of the dynamics of cancer virotherapy with an immune response were considered. For further study of attractive models related to tumor dynamics, we refer the reader to [35,36,37]. It is suitable to analyze any biological system’s qualitative behaviour by discrete-time systems compared to their alternatives in differential equations. Additionally, there is superior observation and investigation of chaos in all biological systems through discrete-time mathematical systems [38]. Hence, it is motivating to study the qualitative behaviour of the system (3) in its discrete form. In recent times, numerous scientific approaches have been presented to discretize any scientific model from continuous time. The traditional method is to use ordinary difference systems such as Euler’s approximations and Runge–Kutta methods to attain this objective.

Nevertheless, mathematical unpredictability is experienced using traditional finite difference approaches. To escape from this mathematical unpredictability it is possible to use a nonstandard finite difference technique, as specified by Mickens [39]. Generally, a nonstandard finite difference scheme is directed to protect the following characteristics of the corresponding continuous-time system: boundedness, the positivity of results, stability of fixed points, and bifurcations. The development of these varieties of difference methods is not straightforward, and no usual methods can be found in the literature for building them, possibly reflecting a chief disadvantage of nonstandard difference methods. Therefore, by applying a Mickens-type nonstandard scheme to model (3), we obtain the following discrete-time mathematical model (see [39,40]):

$$\left\{\begin{array}{c}\frac{x_{n+1}-x_n}{h}=\sigma +\frac{\rho x_n{y}_n}{\eta +y_n}-\mu x_n{y}_n-\delta x_n,\hfill \\ \frac{y_{n+1}-y_n}{h}=\alpha y_n(1-\beta y_n)-x_n{y}_n\hfill \end{array}\right.$$

(4)

where $h\in [0,1)$ is the step size for discretization. Furthermore, system (4) can be transformed into the following form:

$$\left\{\begin{array}{c}{x}_{n+1}=\frac{x_n+h\left(\sigma +\frac{\rho x_n{y}_n}{\eta +y_n}\right)}{1+h\left(\delta +\mu y_n\right)},\hfill \\ y_{n+1}=\frac{y_n\left(1+h\alpha \right)}{1+h\left(x_n+\alpha \beta y_n\right)}.\hfill \end{array}\right.$$

(5)

The subsequent parts of this manuscript are directed at: Andronov–Hopf bifurcation in system (3); the boundedness of solutions of system (5); the existence of fixed points and local stability analysis of system (5); the presence and direction of Neimark–Sacker bifurcation about the positive fixed point of system (5); the control of Neimark–Sacker bifurcation in system (5); and finally, several numerical simulations which support our theoretical discussion.

2. Andronov–Hopf Bifurcation

Let $(x^{*},y^{*})$ be the positive fixed point of map (3). First, we explore the behavior of continuous system (3). In order to explore the behavior of system (3), the Jacobian matrix $V\left(E\right)$ at $(x^{*},y^{*})$ is provided by

$$V\left(E\right)=\begin{array}{c}\left[\begin{array}{cc}-\frac{\sigma }{x^{*}}& x^{*}\left(\frac{\eta \rho }{{(y^{*}+\eta )}^2}-\mu \right)\\ -y^{*}& \alpha -x^{*}-2y^{*}\alpha \beta \end{array}\right].\end{array}$$

where $detV\left(E\right)=x^{*}{y}^{*}\left(\frac{\eta \rho }{{(y^{*}+\eta )}^2}-\mu \right)+\sigma +\frac{\alpha (2\beta y^{*}-1)\sigma }{x^{*}}$ and $TrV\left(E\right)=\alpha -2\alpha \beta y^{*}-\frac{\sigma }{x^{*}}-x^{*}.$ Hence, per the Routh–Hurwitz stability criterion, $(x^{*},y^{*})$ is sink if and only if $x^{*}{y}^{*}\left(\frac{\eta \rho }{{\left(y^{*}+\eta \right)}^2}-\mu \right)+\sigma >\frac{\alpha \left(1-2\beta {\mathrm{y}}^{*}\right)\sigma }{x^{*}}$ and $\alpha \left(1-2\beta {\mathrm{y}}^{*}\right)<x^{*}+\frac{\sigma }{x^{*}}$ and source if and only if $x^{*}{y}^{*}\left(\frac{\eta \rho }{{\left(y^{*}+\eta \right)}^2}-\mu \right)+\sigma >\frac{\alpha \left(1-2\beta {\mathrm{y}}^{*}\right)\sigma }{x^{*}}$ and $\alpha \left(1-2\beta {\mathrm{y}}^{*}\right)>x^{*}+\frac{\sigma }{x^{*}}$. Moreover, system (3) experiences Hopf bifurcation if the parameters lie in the following curve:

$$T_{HB}=\{(\alpha ,\beta ,\delta ,\sigma ,\rho ,\mu ,\eta )\in R_{+}^7:\sigma =x^{*}\left(\alpha -x^{*}-2y^{*}\alpha \beta \right)\}.$$

Furthermore, Figure 2 shows the topological classification of the fixed point $(x^{*},y^{*})$. To explore the periodic nature of solutions of system (3), we study the exitances of subcritical and supercritical Hopf bifurcation. For this, we assume the next planar system:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{dx}{dt}=& f(a,x,y),\hfill \\ \hfill \frac{dy}{dt}=& g(a,x,y),\hfill \end{array}\end{array}$$

(6)

where $a\in R$ is the bifurcation parameter. Let $V(x^{*},y^{*})$ be the Jacobian matrix of (6) computed at equilibrium point $(x^{*},y^{*})$. Moreover, the eigenvalues of (6) computed at any equilibrium point $(x^{*},y^{*})$ are of the following form:

$${\lambda }_{1,2}\left(a\right)=\varphi \left(a\right)\pm \iota \phi \left(a\right).$$

Furthermore, we assume that there is a particular value of the bifurcation parameter a, say $a_0$, for which the following conditions hold true (see [40]):

 (i) 

$\varphi \left(a_0\right)=0$ and $\phi \left(a_0\right)={\phi }_0\ne 0$. Then, there exists a conjugate pair of complex eigenvalues of $V(x^{*},y^{*})$ in the condition of non-hyperbolicity.

 (ii) 

$\frac{d\varphi \left(a\right)}{da}{|}_{a=a_0}=T\ne 0$, which is known as a transversality condition, that is, the eigenvalues of $V(x^{*},y^{*})$ cross the imaginary axis with non-zero speed [40].

 (iii) 

There exists a discriminatory quantity $L\left(a_0\right)\ne 0$, which is known as the first Lyapunov exponent (FLE) and is defined as follows:

$$L\left(a_0\right)=L_1+L_2,$$

where

$$\left\{\begin{array}{c}{L}_1=\frac{1}{16}\left[f_{xxx}+f_{xyy}+g_{xxy}+g_{yyy}\right]\hfill \\ L_2=\frac{1}{16{\phi }_0}\left[f_{xy}(f_{xx}+f_{yy})-g_{xy}(g_{xx}+g_{yy})-f_{xx}{g}_{xx}+f_{yy}{g}_{yy}\right]\hfill \end{array}\right.$$

with $f_{xy}=\frac{{\partial }^{2}f(a,x,y)}{\partial x\partial y}$ computed at $(x,y)=(x^{*},y^{*})$ and $a=a_0.$

Figure 2. Topological classification of the one and only fixed point of system (3) for $0<\alpha <2,0<\beta <1,\delta =0.3743,\eta =20.19,\mu =0.00311,\rho =11.131,$ and $0<\sigma <3.5$ with initial conditions $x_0=1.6197$ and $y_0=0.82317$.

Figure 2. Topological classification of the one and only fixed point of system (3) for $0<\alpha <2,0<\beta <1,\delta =0.3743,\eta =20.19,\mu =0.00311,\rho =11.131,$ and $0<\sigma <3.5$ with initial conditions $x_0=1.6197$ and $y_0=0.82317$.

Theorem 1

([40]). Assume that conditions (i), (ii), and (iii) are satisfied; then, there exists a unique curve of periodic solutions. Moreover, this curve bifurcates from the fixed point into the region with $a>a_0$ if $L\left(a_0\right)<0$ or $L\left(a_0\right)>0$ if $a<a_0.$ In addition, the fixed point is stable whenever $a>a_0$ (respectively, $a<a_0$) and unstable for $a<a_0$ (respectively, $a>a_0$) for $T<0$ and $T>0$, respectively.

It can be seen that periodic solutions are stable (respectively, unstable) if the fixed point is unstable (respectively, stable) on the side of $a=a_0$ where the periodic solutions exist. Keeping in view the above discussion for Andronov–Hopf bifurcation, we consider the system (3). For this, we assume that $4\left(\sigma +\frac{\alpha \sigma \left(2\beta y^{*}-1\right)}{x^{*}}+x^{*}{y}^{*}\left(\frac{\eta \rho }{{\left(\eta +y^{*}\right)}^2}-\mu \right)\right)-{\left(\frac{\sigma }{x^{*}}+x^{*}+\alpha \left(2\beta y^{*}-1\right)\right)}^2>0$; then, it is easy to see that the eigenvalues of the Jacobian matrix $V(x^{*},y^{*})$ are of the form

$${\lambda }_{1,2}\left(\sigma \right)=\varphi \left(\sigma \right)\pm \iota \phi \left(\sigma \right),$$

$$\varphi \left(\sigma \right)=\frac{1}{2}\left(\alpha -\frac{\sigma }{x^{*}}-x^{*}-2\alpha \beta y^{*}\right)$$

and

$$\phi \left(\sigma \right)=\frac{1}{2}\sqrt{4\left(\sigma +\frac{\alpha \sigma \left(2\beta y^{*}-1\right)}{x^{*}}+x^{*}{y}^{*}\left(\frac{\eta \rho }{{\left(\eta +y^{*}\right)}^2}-\mu \right)\right)-{\left(\frac{\sigma }{x^{*}}+x^{*}+\alpha \left(2\beta y^{*}-1\right)\right)}^2}.$$

Next, $\varphi \left(\sigma \right)=0$ provides $\sigma ={\sigma }_0=x^{*}\left(\alpha -x^{*}-2\alpha \beta y^{*}\right)$. At $\sigma ={\sigma }_0$, we have

$$\phi \left({\sigma }_0\right)={\phi }_0=\sqrt{x^{*}\left(2\alpha +y^{*}\left(\frac{\eta \rho }{{\left(\eta +y^{*}\right)}^2}-4\alpha \beta -\mu \right)\right)-{\left(x^{*}\right)}^2-{\alpha }^2{\left(1-2\beta y^{*}\right)}^2}\ne 0.$$

For the transversality condition, we can see that

$${\frac{d\varphi \left(\sigma \right)}{d\sigma }|}_{\sigma ={\sigma }_0}=-\frac{1}{2x^{*}}<0.$$

In order to shift the fixed point of system (3) to origin $(0,0)$, we consider the following translations:

$$u\left(t\right)=x\left(t\right)-x^{*},v\left(t\right)=y\left(t\right)-y^{*}.$$

Moreover, by implementing this transformation on system (3), we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{du}{dt}=& \sigma +\frac{\rho (u+x^{*})(v+y^{*})}{\eta +y}-\mu (u+x^{*})(v+y^{*})-\delta (u+x^{*}),\hfill \\ \hfill \frac{dv}{dt}=& \alpha (v+y^{*})(1-\beta (v+y^{*}))-(u+x^{*})(v+y^{*}),\hfill \end{array}\end{array}$$

(7)

Application of Taylor series expansion on $(u,v)=(0,0)$ provides the following system:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{c}\left[\begin{array}{c}\frac{du}{dt}\\ \frac{dv}{dt}\end{array}\right]=B\left[\begin{array}{c}u\\ v\end{array}\right]+\left[\begin{array}{c}\left(\frac{1}{\eta +y^{*}}\left(\rho -\frac{\rho \left(\eta y^{*}+{y^{*}}^2\right)}{{\left(\eta +y^{*}\right)}^2}\right)-\mu \right)vu-\frac{\rho x^{*}{v}^2\eta }{{\left(\eta +y^{*}\right)}^3}-\frac{\rho \eta v^{2}u}{{\left(\eta +y^{*}\right)}^3}+\frac{\rho v^3\eta x^{*}}{{\left(\eta +y^{*}\right)}^4}\\ -\alpha v^2\beta -uv\end{array}\right]\end{array}\end{array}\end{array}$$

(8)

where

$$B=\begin{array}{c}\left[\begin{array}{cc}-\frac{\sigma }{x^{*}}& x^{*}\left(\frac{\eta \rho }{{(y^{*}+\eta )}^2}-\mu \right)\\ -y^{*}& \alpha -x^{*}-2y^{*}\alpha \beta \end{array}\right].\end{array}$$

Next, we want to convert matrix B into its canonical form. For this, the following similarity transformation is considered:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{c}\left[\begin{array}{c}u\left(t\right)\\ v\left(t\right)\end{array}\right]=\left[\begin{array}{cc}{x}^{*}\left(\frac{\eta \rho }{{(y^{*}+\eta )}^2}-\mu \right)& 0\\ \frac{\sigma }{x^{*}}& \sqrt{x^{*}{y}^{*}\left(\frac{\eta \rho }{{(y^{*}+\eta )}^2}-\mu \right)+\sigma +\frac{\alpha (2y^{*}\beta -1)\sigma }{x^{*}}}\end{array}\right]\left[\begin{array}{c}w\left(t\right)\\ z\left(t\right)\end{array}\right]\end{array}\end{array}\end{array}$$

(9)

From (8) and (9), it follows that

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{c}\left[\begin{array}{c}\frac{dw}{dt}\\ \frac{dz}{dt}\end{array}\right]=\left[\begin{array}{cc}0& {\phi }_0\\ -{\phi }_0& 0\end{array}\right]\left[\begin{array}{c}w\left(t\right)\\ z\left(t\right)\end{array}\right]\end{array}\end{array}\end{array}$$

(10)

where

$$f(w,z)=\frac{uv{(y^{*}+\eta )}^4\mu -v(v-y^{*}-\eta )\eta (v{x}^{*}-u(y^{*}+\eta ))\rho }{x^{*}{(y^{*}+\eta )}^4\mu -x^{*}\eta {(y^{*}+\eta )}^2\rho },$$

and

$$g(w,z)=\frac{v\left(\frac{(x^{*}+\alpha (2y^{*}\beta -1))\left(u{(y^{*}+\eta )}^4\mu -\eta (-v+y^{*}+\eta )(-v{x}^{*}+u(y^{*}+\eta ))\rho \right)}{x^{*}{(y^{*}+\eta )}^4\mu -x^{*}\eta {(y^{*}+\eta )}^2\rho }-u-v\alpha \beta \right)}{\sqrt{x^{*}\left(\alpha (2-4y^{*}\beta )-y^{*}\mu +\frac{y^{*}\eta \rho }{{(y^{*}+\eta )}^2}\right)-{x^{*}}^2-{\alpha }^2{(1-2y^{*}\beta )}^2}},$$

$$u=x^{*}\left(\frac{\eta \rho }{{(y^{*}+\eta )}^2}-\mu \right)w,v=\left(\frac{\sigma }{x^{*}}\right)w+\left(\sqrt{x^{*}{y}^{*}\left(\frac{\eta \rho }{{(y^{*}+\eta )}^2}-\mu \right)+\sigma +\frac{\alpha (2y^{*}\beta -1)\sigma }{x^{*}}}\right)z.$$

Then, the first Lyapunov exponent for system (3) is computed as follows:

$$L\left({\sigma }_0\right)=L_1+L_2,$$

where

$$\left\{\begin{array}{c}{L}_1=\frac{1}{16}\left(\frac{2\rho y^{*}}{{\left(\eta +y^{*}\right)}^3}-\frac{2\rho }{{\left(\eta +y^{*}\right)}^2}\right)\hfill \\ L_2=\frac{-2\alpha \beta -2\alpha \beta \left(\frac{2\rho x^{*}{y}^{*}}{{\left(\eta +y^{*}\right)}^3}-\frac{2\rho x^{*}}{{\left(\eta +y^{*}\right)}^2}\right)+\left(\frac{2\rho x^{*}{y}^{*}}{{\left(\eta +y^{*}\right)}^3}-\frac{2\rho x^{*}}{{\left(\eta +y^{*}\right)}^2}\right)\left(-\mu -\frac{\rho y^{*}}{{\left(\eta +y^{*}\right)}^2}+\frac{\rho }{\eta +y^{*}}\right)}{16\sqrt{-{\left(x^{*}\right)}^2-{\alpha }^2{\left(1-2\beta y^{*}\right)}^2+x^{*}\left(2\alpha +y^{*}\left(-4\alpha \beta -\mu +\frac{\eta \rho }{{\left(\eta +y^{*}\right)}^2}\right)\right)}}.\hfill \end{array}\right.$$

Finally, we have the following theorem:

Theorem 2.

Assume that conditions (i), (ii), and (iii) are satisfied; then, there exists a unique curve of periodic solutions. Moreover, this curve bifurcates from the fixed point into the region with $\sigma >{\sigma }_0$ if $(\rho y^{*}{\left(\eta +y^{*}\right)}^2-\rho {\left(\eta +y^{*}\right)}^3)\sqrt{-{\left(x^{*}\right)}^2-{\alpha }^2{\left(1-2\beta y^{*}\right)}^2+x^{*}\left(2\alpha +y^{*}\left(-4\alpha \beta -\mu +\frac{\eta \rho }{{\left(\eta +y^{*}\right)}^2}\right)\right)}-$ $\alpha \beta {\left(\eta +y^{*}\right)}^5+\eta \rho x^{*}\left(\eta (2\alpha \beta \eta +\eta \mu -\rho )+(2\alpha \beta +\mu )y^{*}\left(2\eta +y^{*}\right)\right)<0$ or $(\rho y^{*}{\left(\eta +y^{*}\right)}^2-\rho {\left(\eta +y^{*}\right)}^3)\sqrt{-{\left(x^{*}\right)}^2-{\alpha }^2{\left(1-2\beta y^{*}\right)}^2+x^{*}\left(2\alpha +y^{*}\left(-4\alpha \beta -\mu +\frac{\eta \rho }{{\left(\eta +y^{*}\right)}^2}\right)\right)}-$ $\alpha \beta {\left(\eta +y^{*}\right)}^5+\eta \rho x^{*}\left(\eta (2\alpha \beta \eta +\eta \mu -\rho )+(2\alpha \beta +\mu )y^{*}\left(2\eta +y^{*}\right)\right)>0$ if $\sigma <{\sigma }_0.$ In addition, the fixed point is stable whenever $\sigma >{\sigma }_0$ (respectively, $\sigma <{\sigma }_0$) and unstable for $\sigma <{\sigma }_0$ (respectively, $\sigma >{\sigma }_0$) for $T<0$ and $T>0$, respectively.

The plot of FLE is depicted in Figure 3.

3. Boundedness of Solutions

From the second equation of system (5) it follows that

$$y_{n+1}\le \frac{y_n\left(1+h\alpha \right)}{1+h\alpha \beta y_n}.$$

(11)

Consequently, by solving (11) and then applying the limit, we obtain

$$\underset{n⟶\infty }{lim}sup{y}_n\le \frac{1}{\beta },$$

(12)

for all $n\ge 0.$ In the same way, from the first equation of system (5) we obtain

$$x_{n+1}=\frac{x_n+h\left(\sigma +\frac{\rho x_n{y}_n}{\eta +y_n}\right)}{1+h\left(\delta +\mu y_n\right)}$$

$$\le \frac{x_n+h\left(\sigma +\frac{\rho x_n}{1+\eta \beta }\right)}{1+h\delta },$$

(13)

because

$$x_n+h\left(\sigma +\frac{\rho x_n{y}_n}{\eta +y_n}\right)\le x_n+h\left(\sigma +\frac{\rho x_n}{1+\eta \beta }\right)$$

yields

$$y_n\le \frac{1}{\beta }$$

for all $n\ge 0.$ Then, from (13) we have

$$x_{n+1}\le \frac{x_n(1+\frac{h\rho }{1+\eta \beta })+h\sigma }{1+h\delta }.$$

Hence, we can obtain the upper bound for $x_n$ as

$$\underset{n\to \infty }{lim}sup{x}_n\le \frac{\sigma +\beta \eta \sigma }{\delta +\beta \delta \eta -\rho },$$

(14)

for all $n\ge 0.$ Finally, we have the following theorem concerning the boundedness of all solutions of (5).

Theorem 3.

Assume that $0<y_0\le \frac{1}{\beta }$ and $0<x_0\le \frac{\sigma +\beta \eta \sigma }{\delta +\beta \delta \eta -\rho }$; then, for all $n\ge 0,$ every positive solution $(x_n,y_n)$ of system (5) is bounded and contained in the set $\left[0,\frac{\sigma +\beta \eta \sigma }{\delta +\beta \delta \eta -\rho }\right]\times \left[0,\frac{1}{\beta }\right]$ whenever $\rho <\delta (1+\beta \eta ).$

4. Existence of Fixed Points and Local Stability Analysis

It is easy to see that systems (3) and (5) have more than one fixed point by performing simple algebraic manipulation. Moreover, most of them are complex and have no biological significance. Here, the fixed points with biological significance are boundary fixed points and the unique positive fixed point. In order to obtain those significant fixed points of system (5), we consider the following system of equations:

$$\left\{\begin{array}{c}{x}^{*}=\frac{x^{*}+h\left(\sigma +\frac{\rho x^{*}{y}^{*}}{\eta +y^{*}}\right)}{1+h\left(\delta +\mu y^{*}\right)},\hfill \\ y^{*}=\frac{y^{*}\left(1+h\alpha \right)}{1+h\left(x^{*}+\alpha \beta y^{*}\right)}.\hfill \end{array}\right.$$

(15)

Furthermore, from (15) we obtain the following pair:

$$\delta x^{*}+\mu y^{*}=\sigma +\frac{\rho x^{*}{y}^{*}}{\eta +y^{*}},y^{*}=\frac{\alpha -x^{*}}{\alpha \beta }.$$

Moreover, we have two fixed points from (15), namely, $(\frac{\sigma }{\delta },0)$ and $(x^{*},y^{*}),$ where $x^{*}$ is solution of the following equation:

$$a_{11}{x}^2+a_{12}x+a_{13}=0,$$

(16)

with

$$\left\{\begin{array}{c}{a}_{11}=\mu +\alpha \beta (\rho -\delta ),,a_{12}=\alpha \beta (\alpha \beta \eta \delta -\eta \mu -\alpha +\sigma +\alpha \rho )-2\alpha \mu \hfill \\ and{a}_{13}=\eta ({\alpha }^2\beta \mu -\sigma )-{\alpha }^2(\beta \sigma -\mu ).\hfill \end{array}\right.$$

(17)

In addition, $a_{11},a_{12}>0$ and $a_{13}<0$ if we have

$$\left\{\begin{array}{c}\rho >\delta ,,\sigma +\alpha (\beta \eta \delta +\rho )<\eta \mu +\alpha (1+2\mu )\hfill \\ and{\alpha }^2\mu (\beta \eta +1)<\sigma (\eta +{\alpha }^2\beta ).\hfill \end{array}\right.$$

(18)

Hence, using Descartes’ rule of signs (see [41]), we have the following result:

Theorem 4.

Assume that $0<y_0\le \frac{1}{\beta }$ and $0<x_0\le \frac{\sigma +\beta \eta \sigma }{\delta +\beta \delta \eta -\rho }$; then, for

$$\left\{\begin{array}{c}\rho >\delta ,,\sigma +\alpha (\beta \eta \delta +\rho )<\eta \mu +\alpha (1+2\mu )\hfill \\ and{\alpha }^2\mu (\beta \eta +1)<\sigma (\eta +{\alpha }^2\beta )\hfill \end{array}\right.$$

(19)

there exists a unique positive constant solution $(x^{*},y^{*})$ of system (5) in $\left[0,\frac{\sigma +\beta \eta \sigma }{\delta +\beta \delta \eta -\rho }\right]\times \left[0,\frac{1}{\beta }\right]$ if and only if $x^{*}<\alpha$ for each $y^{*}\in \left[0,\frac{1}{\beta }\right]$.

In order to discuss the stability of system (5) about these fixed points, we compute the Jacobian matrix $J_{(x,y)}$ of system (5) about each of its fixed points $(x,y)$. Moreover, $J_{(x,y)}$ is

$$J_{(x,y)}=\begin{array}{c}\left[\begin{array}{cc}{m}_{11}& m_{12}\\ m_{21}& m_{22}\end{array}\right]\end{array}.$$

In addition, the characteristic polynomial $\mathbb{C}\left(\varphi \right)$ of $J_{(x,y)}$ is

$$\mathbb{C}\left(\varphi \right)={\varphi }^2-T\varphi +D,$$

(20)

where T and D represent the trace and determinants of $J_{(x,y)}$, respectively. The following lemma describes the condition equivalent to the Jury conditions for the local stability of fixed points (see [42]).

Lemma 1

([42]). Let $\mathbb{H}\left(\varphi \right)={\varphi }^2-T\varphi +D$ be the characteristic equation obtained from a $2\times 2$ Jacobian matrix $J_{(x,y)}.$ Moreover, let $J_{(x,y)}$ be any Jacobian matrix of system (5) about each of its equilibrium points. Additionally, assume that $\mathbb{H}\left(1\right)>0$. Then:

(i) $|{\varphi }_1|<1$ and $|{\varphi }_2|<1$ if and only if $\mathbb{H}(-1)>0$ and $D<1$

(ii) $|{\varphi }_1|>1$ and $|{\varphi }_2|>1$ if and only if $\mathbb{H}(-1)>0$ and $D>1$

(iii) $|{\varphi }_1|<1$ and $|{\varphi }_2|>1$ or $\left(\right|{\varphi }_1|>1$ and $|{\varphi }_2<\left|1\right)$ if and only if $\mathbb{H}(-1)<0$

(iv) ${\varphi }_1$ and ${\varphi }_2$ represent complex conjugates with $|{\varphi }_1|=1=|{\varphi }_2|$ if and only if $T^2-4D<0$ and $D=1$.

As ${\varphi }_1$ and ${\varphi }_2$ are characteristic values of (20), the point $(x,y)$ is sink if $|{\varphi }_1|<1$ and $|{\varphi }_2|<1$. Furthermore, it is locally asymptotically stable. The point $(x,y)$ is known as source(repeller) if $|{\varphi }_1|>1$ and $|{\varphi }_2|>1$. The point $(x,y)$ is a saddle point if $|{\varphi }_1|<1$ and $|{\varphi }_2|>1$ or $\left(|{\varphi }_1|>1and|{\varphi }_2|<1\right)$. Finally, $(x,y)$ is non-hyperbolic if condition $\left(iv\right)$ is satisfied.

First, we study the stability conditions for (5) about the fixed point $(\frac{\sigma }{\delta },0).$ The matrix $J_{(x,y)}$ evaluated at $(\frac{\sigma }{\delta },0)$ is provided by

$$J_{(\frac{\sigma }{\delta },0)}=\left[\begin{array}{cc}\frac{1}{1+h\delta }& \frac{h(\rho -\eta \mu )\sigma }{\delta (1+h\delta )\eta }\\ 0& \frac{\delta +h\alpha \delta }{\delta +h\sigma }\end{array}\right].$$

Furthermore, the eigenvalues of $J_{(\frac{\sigma }{\delta },0)}$ are ${\xi }_1=\frac{1}{1+h\delta }$ and ${\xi }_2=\frac{\delta +h\alpha \delta }{\delta +h\sigma }$ with $|{\xi }_1|<1$ for all parametric values. Hence, we have the following result related to the dynamics of (5) about $(\frac{\sigma }{\delta },0)$:

Proposition 1.

The boundary equilibrium $(\frac{\sigma }{\delta },0)$ of system (5) is source and saddle $iff$ conditions $\alpha \delta <\sigma$ and $\alpha \delta >\sigma$ are satisfied respectively (see Figure 4).

Next, our task is to explore of the local stability of system (5) about the point $(x^{*},y^{*})$. Let $J_{(x^{*},y^{*})}$ be the Jacobian matrix of system (5) about the fixed point $(x^{*},y^{*})$; then, $J_{(x^{*},y^{*})}$ has the following mathematical form:

$$J_{(x^{*},y^{*})}=\left[\begin{array}{cc}\frac{1+\frac{h\rho y^{*}}{\eta +y^{*}}}{1+h\delta +h\mu y^{*}}& -\frac{h\left(h\mu \sigma {\left(\eta +y^{*}\right)}^2+x^{*}\left(\eta (\eta \mu -(1+h\delta )\rho )+\mu y^{*}\left(2\eta +y^{*}+h\rho y^{*}\right)\right)\right)}{{\left(\eta +y^{*}\right)}^2{\left(1+h\delta +h\mu y^{*}\right)}^2}\\ -\frac{h{y}^{*}}{1+h\alpha }& 1-\frac{h\alpha \beta y^{*}}{1+h\alpha }\end{array}\right].$$

Moreover, from $J_{(x^{*},y^{*})}$ we can obtain the following characteristic polynomial:

$$\begin{array}{ccc}\hfill \mathbb{C}\left(\varphi \right)& =& {\varphi }^2-\varphi \left(1-\frac{h\alpha \beta y^{*}}{1+h\alpha }+\frac{1+\frac{h\rho y^{*}}{\eta +y^{*}}}{1+h\delta +h\mu y^{*}}\right)+D,\hfill \end{array}$$

(21)

where

$$\begin{array}{ccc}\hfill D& =& \frac{1}{1+h\delta +h\mu y^{*}}+\frac{h{y}^{*}}{1+h\alpha }\left(\frac{S{(1+h\alpha )}^2}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}-\frac{\alpha \beta }{1+h\delta +h\mu y^{*}}+\frac{\rho +h\alpha \rho -h\alpha \beta \rho y^{*}}{\left(\eta +y^{*}\right)\left(1+h\delta +h\mu y^{*}\right)}\right)\hfill \end{array}$$

and

$$\begin{array}{c}\hfill S=-\frac{h\left(h\mu \sigma {\left(\eta +y^{*}\right)}^2+x^{*}\left(\eta (\eta \mu -(1+h\delta )\rho )+\mu y^{*}\left(2\eta +y^{*}+h\rho y^{*}\right)\right)\right)}{{\left(\eta +y^{*}\right)}^2{\left(1+h\delta +h\mu y^{*}\right)}^2}.\end{array}$$

By considering (19) and taking

$$\left\{\begin{array}{c}{s}_1=h\left(h\mu \sigma {\left(\eta +y^{*}\right)}^2+x^{*}\left(\eta (\eta \mu -(1+h\delta )\rho )+\mu y^{*}\left(2\eta +y^{*}+h\rho y^{*}\right)\right)\right),\hfill \\ s_2={\left(\eta +y^{*}\right)}^2{\left(1+h\delta +h\mu y^{*}\right)}^2,\hfill \end{array}\right.$$

(22)

it follows that

$$\mathbb{C}\left(1\right)=\frac{h{y}^{*}\left((S+h(S+\alpha \beta )\delta )\eta +y^{*}\left(S(1+h\delta +h\eta \mu )+h\alpha \beta (\delta +\eta \mu -\rho )+h(S+\alpha \beta )\mu y^{*}\right)\right)}{(1+h\alpha )\left(\eta +y^{*}\right)\left(1+h\delta +h\mu y^{*}\right)}>0,$$

for

$$h{s}_2\alpha \beta \left(\delta (1+\eta )+\mu (\eta +y^{*})\right)>s_1\left(\eta (1+h\delta )+y^{*}(1+h(\delta +2\mu \eta ))\right)+h{s}_2\alpha \beta \rho .$$

Moreover, we have

$$\begin{array}{c}\hfill h=\frac{\delta \eta +y^{*}\left(\delta +\eta (\alpha \beta +\mu -S)-\rho +(\alpha \beta +\mu -S)y^{*}\right)}{y^{*}\left(S\delta \eta -\alpha (\delta +\eta \mu -\rho )+y^{*}\left(S(\delta +\eta \mu )-\alpha (\mu +\beta \rho )+S\mu y^{*}\right)\right)-\alpha \delta \eta }\end{array}$$

(23)

and

$$\begin{array}{ccc}\hfill \mathbb{C}(-1)& =& 1+\frac{(1+h\alpha )\left(1+h{x}^{*}\right)}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}+\frac{1+\frac{h\rho y^{*}}{\eta +y^{*}}}{1+h\delta +h\mu y^{*}}+\frac{1}{1+h\delta +h\mu y^{*}}\hfill \\ & +& \frac{h{y}^{*}}{1+h\alpha }\left(\frac{S{(1+h\alpha )}^2}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}-\frac{\alpha \beta }{1+h\delta +h\mu y^{*}}+\frac{\rho +h\alpha \rho -h\alpha \beta \rho y^{*}}{\left(\eta +y^{*}\right)\left(1+h\delta +h\mu y^{*}\right)}\right).\hfill \end{array}$$

Hence, for the study of the linearized stability of system (5) about $(x^{*},y^{*}),$ we have the following proposition (see Lemma 1).

Proposition 2.

Let (18) remain true; then, $(x^{*},y^{*})$ is positive equilibrium of (5). Moreover, let us assume that

$$\left\{\begin{array}{c}{\Omega }_{11}=1+h{x}^{*}+h\alpha \beta y^{*},{\Omega }_{12}=1+h\delta +h\mu y^{*},{\Omega }_{13}=\eta +y^{*},\hfill \\ {\Omega }_{14}=1+h\alpha ,\hfill \end{array}\right.$$

(24)

then,

• The fixed point $(x^{*},y^{*})$ is stable if and only if for

  $$\begin{array}{c}\eta x^{*}\left(\rho +2\mu y^{*}\right)<h\mu \sigma x^{*}{\left(y^{*}\right)}^2+\mu \left(h\sigma +x^{*}\right){\left(\eta +y^{*}\right)}^2\hfill \end{array}$$

  we have

  $$\begin{array}{c}{\Omega }_{11}^2{\Omega }_{12}{\Omega }_{13}>{\Omega }_{14}\left({\Omega }_{11}\left({\Omega }_{13}+h\rho y^{*}\right)+h{y}^{*}\left(\delta \eta +(hS\delta -\alpha \beta ){\Omega }_{13}+hS\eta \mu y^{*}-\alpha \beta \rho y^{*}\right)\right)\hfill \end{array}$$
• The fixed point $(x^{*},y^{*})$ is non-hyperbolic if and only if

  $$h=\frac{\delta \eta +y^{*}\left(\delta +\eta (\alpha \beta +\mu -S)-\rho +(\alpha \beta +\mu -S)y^{*}\right)}{y^{*}\left(S\delta \eta -\alpha (\delta +\eta \mu -\rho )+y^{*}\left(S(\delta +\eta \mu )-\alpha (\mu +\beta \rho )+S\mu y^{*}\right)\right)-\alpha \delta \eta }$$

  and

  $$\begin{array}{c}1+\frac{{\Omega }_{13}+h\rho y^{*}}{{\Omega }_{13}{\Omega }_{12}}<\frac{h\alpha \beta y^{*}}{{\Omega }_{14}}+4\left(\frac{1}{{\Omega }_{12}}+\frac{h{y}^{*}}{{\Omega }_{14}}\left(\frac{S{\Omega }_{14}^2}{{\left({\Omega }_{11}\right)}^2}-\frac{\alpha \beta }{{\Omega }_{12}}+\frac{\rho +h\alpha \rho -h\alpha \beta \rho y^{*}}{\left(\eta +y^{*}\right)\left({\Omega }_{12}\right)}\right)\right).\hfill \end{array}$$

  (25)

Remark 1.

Let (18) remain true; then, there is no chance of system (5) undergoing period-doubling bifurcation, as $\mathbb{C}(-1)>0$ for every $\sigma ,\rho ,\eta ,\delta ,\mu ,\alpha ,\beta >0$ and $\rho >\delta .$

We now have the following theorem for the possible validation of Remark 1.

Theorem 5.

Let (18) remain true and let $\rho >\delta .$ Moreover, let $(x^{*},y^{*})$ be the positive fixed point of system (5) and

$$\begin{array}{ccc}\hfill \mathbb{C}(-1)& =& 1+\frac{(1+h\alpha )\left(1+h{x}^{*}\right)}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}+\frac{1+\frac{h\rho y^{*}}{\eta +y^{*}}}{1+h\delta +h\mu y^{*}}+\frac{1}{1+h\delta +h\mu y^{*}}\hfill \\ & +& \frac{h{y}^{*}}{1+h\alpha }\left(\frac{S{(1+h\alpha )}^2}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}-\frac{\alpha \beta }{1+h\delta +h\mu y^{*}}+\frac{\rho +h\alpha \rho -h\alpha \beta \rho y^{*}}{\left(\eta +y^{*}\right)\left(1+h\delta +h\mu y^{*}\right)}\right).\hfill \end{array}$$

Then, $\mathbb{C}(-1)>0$ for every $\sigma ,\rho ,\eta ,\delta ,\mu ,\alpha ,\beta >0.$

Proof. 

Assume (18) and $\rho >\delta .$ Additionally, let $(x^{*},y^{*})$ be the positive fixed point of system (5) and

$$\begin{array}{ccc}\hfill \mathbb{C}(-1)& =& 1+\frac{(1+h\alpha )\left(1+h{x}^{*}\right)}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}+\frac{1+\frac{h\rho y^{*}}{\eta +y^{*}}}{1+h\delta +h\mu y^{*}}+\frac{1}{1+h\delta +h\mu y^{*}}\hfill \\ & +& \frac{h{y}^{*}}{1+h\alpha }\left(\frac{S{(1+h\alpha )}^2}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}-\frac{\alpha \beta }{1+h\delta +h\mu y^{*}}+\frac{\rho +h\alpha \rho -h\alpha \beta \rho y^{*}}{\left(\eta +y^{*}\right)\left(1+h\delta +h\mu y^{*}\right)}\right).\hfill \end{array}$$

Then, from

$$\begin{array}{c}\hfill S=-\frac{h\left(h\mu \sigma {\left(\eta +y^{*}\right)}^2+x^{*}\left(\eta (\eta \mu -(1+h\delta )\rho )+\mu y^{*}\left(2\eta +y^{*}+h\rho y^{*}\right)\right)\right)}{{\left(\eta +y^{*}\right)}^2{\left(1+h\delta +h\mu y^{*}\right)}^2},\end{array}$$

we have

$$\begin{array}{c}\hfill h\mu \sigma {\left(\eta +y^{*}\right)}^2+x^{*}\left(\eta (\eta \mu -(1+h\delta )\rho )+\mu y^{*}\left(2\eta +y^{*}+h\rho y^{*}\right)\right)>0\end{array}$$

$$\begin{array}{c}\hfill \iff h\mu \sigma y^{*}\left(\eta +x^{*}{y}^{*}\right)+\mu \left(h\sigma +x^{*}\right)\left({\eta }^2+\eta y^{*}+{\left(y^{*}\right)}^2\right)>\eta x^{*}\left(\rho +\mu y^{*}\right)\end{array}$$

$$\begin{array}{c}\hfill \iff \eta x^{*}\left(\rho +2\mu y^{*}\right)<h\mu \sigma x^{*}{\left(y^{*}\right)}^2+\mu \left(h\sigma +x^{*}\right){\left(\eta +y^{*}\right)}^2.\end{array}$$

(26)

Finally, under condition (26), we have $\mathbb{C}(-1)>0$ if and only if

$$\begin{array}{c}\hfill \left(\eta +y^{*}\right){\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2\left(1+h\delta +h\mu y^{*}\right)>\end{array}$$

$$\begin{array}{c}\hfill (1+h\alpha )\left(\eta +y^{*}+h\left(x^{*}\left(\eta +y^{*}+h\rho y^{*}\right)+y^{*}\left(S(\eta +h\delta \eta )+\rho +S{y}^{*}\left(1+h\delta +h\eta \mu +h\mu y^{*}\right)\right)\right)\right)\end{array}$$

$$\begin{array}{c}\hfill \iff {\Omega }_{11}^2{\Omega }_{12}{\Omega }_{13}>{\Omega }_{14}\left(h\delta \eta x^{*}+{\Omega }_{13}\left(1+h{x}^{*}\right)+h{y}^{*}\left(hS\eta (\delta +\mu )+\rho +h\rho x^{*}+S{\Omega }_{12}{y}^{*}\right)\right)\end{array}$$

$$\begin{array}{c}\hfill \iff {\Omega }_{11}^2{\Omega }_{12}{\Omega }_{13}>{\Omega }_{14}\left({\Omega }_{11}\left({\Omega }_{13}+h\rho y^{*}\right)+h{y}^{*}\left(\delta \eta +(hS\delta -\alpha \beta ){\Omega }_{13}+hS\eta \mu y^{*}-\alpha \beta \rho y^{*}\right)\right).\end{array}$$

(27)

where, ${\Omega }_{11},{\Omega }_{12},{\Omega }_{13},{\Omega }_{14}$ are defined in (24). Hence, under condition (26), the above inequality (27) is true for every choice of $\sigma ,\rho ,\eta ,\delta ,\mu ,\alpha ,\beta >0$ and $\rho >\delta .$ This completes the proof. □

5. Neimark–Sacker Bifurcation

This section is related to the bifurcation analysis of system (5) about $(x_{*},y_{*}).$ Moreover, all conditions for the existence and positivity of $(x_{*},y_{*})$ are provided in Section 2. The Neimark–Sacker bifurcation in discrete-time mathematical systems corresponds to the Hopf bifurcation in continuous-time systems. For example, when Neimark–Sacker bifurcation is supercritical, a stable centre loses its stability. A parameter, namely, the bifurcation parameter, is varied with the resulting birth of an established quasi-cycle or cycle. Moreover, we mention any of these as an invariant closed curve. Additionally, for subcritical Neimark–Sacker bifurcation, a stable centre bounded by an unstable closed arc loses its stability through a resulting vanishing of the invariant closed curve as a bifurcation parameter is varied. Here, we discuss the Neimark–Sacker bifurcation experienced by system (5) about $(x_{*},y_{*})$ under certain mathematical conditions. For further study of bifurcation theory and to better understand this surprising behaviuor of discrete-time mathematical systems, we refer readers to [42,43,44,45,46,47,48]. Here, we use the standard theory of bifurcation for study of the Neimark–Sacker bifurcation of system (5) at $(x_{*},y_{*})$. Assume that

$$\begin{array}{c}\hfill S=-\frac{h\left(h\mu \sigma {\left(\eta +y^{*}\right)}^2+x^{*}\left(\eta (\eta \mu -(1+h\delta )\rho )+\mu y^{*}\left(2\eta +y^{*}+h\rho y^{*}\right)\right)\right)}{{\left(\eta +y^{*}\right)}^2{\left(1+h\delta +h\mu y^{*}\right)}^2},\end{array}$$

then, one can see from Proposition 2 that roots ${\varphi }_1,{\varphi }_2$ of (21) are complex and satisfy $|{\varphi }_1|=|{\varphi }_2|=1$ if and only if

$$\begin{array}{c}\hfill h=\frac{\delta \eta +y^{*}\left(\delta +\eta (-S+\alpha \beta +\mu )-\rho +(-S+\alpha \beta +\mu )y^{*}\right)}{-\alpha \delta \eta +y^{*}\left(S\delta \eta -\alpha (\delta +\eta \mu -\rho )+y^{*}\left(S(\delta +\eta \mu )-\alpha (\mu +\beta \rho )+S\mu y^{*}\right)\right)}\end{array}$$

and

$$\begin{array}{c}\hfill 1-\frac{h\alpha \beta y^{*}}{1+h\alpha }+\frac{1+\frac{h\rho y^{*}}{\eta +y^{*}}}{1+h\delta +h\mu y^{*}}<4D,\end{array}$$

where

$$\begin{array}{ccc}\hfill D& =& \frac{\eta +h\alpha \eta +y^{*}\left(1+h(\alpha +S\eta -\alpha \beta \eta +hS\delta \eta +\rho +h\alpha \rho )\right)}{(1+h\alpha )\left(\eta +y^{*}\right)\left(1+h\delta +h\mu y^{*}\right)}\hfill \\ & +& \frac{y^{*}\left(h{y}^{*}\left(S(1+h\delta +h\eta \mu )-\alpha \beta (1+h\rho )+hS\mu y^{*}\right)\right)}{(1+h\alpha )\left(\eta +y^{*}\right)\left(1+h\delta +h\mu y^{*}\right)}.\hfill \end{array}$$

Furthermore, under the suppositions that

$$\left\{\begin{array}{c}\rho >\delta ,,\sigma +\alpha (\beta \eta \delta +\rho )<\eta \mu +\alpha (1+2\mu )\hfill \\ and{\alpha }^2\mu (\beta \eta +1)<\sigma (\eta +{\alpha }^2\beta )\hfill \end{array}\right.$$

(28)

we study the following set:

$${\Psi }_{*}=\left\{\sigma ,\rho ,\eta ,\delta ,\mu ,\alpha ,\beta \in {\Re }^{+},h=\frac{\delta \eta +y^{*}\left(\delta +\eta (-S+\alpha \beta +\mu )-\rho +(-S+\alpha \beta +\mu )y^{*}\right)}{-\alpha \delta \eta +y^{*}\left(S\delta \eta -\alpha (\delta +\eta \mu -\rho )+y^{*}\left(S(\delta +\eta \mu )-\alpha (\mu +\beta \rho )+S\mu y^{*}\right)\right)}\right\}.$$

Then, the positive fixed point $\overline{X}$ of system (5) experiences Neimark–Sacker bifurcation such that h is taken as the bifurcation parameter and varies slightly in the neighborhood of $\widehat{h},$ which is provided as

$$\begin{array}{c}\hfill \widehat{h}=\frac{\delta \eta +y^{*}\left(\delta +\eta (-S+\alpha \beta +\mu )-\rho +(-S+\alpha \beta +\mu )y^{*}\right)}{-\alpha \delta \eta +y^{*}\left(S\delta \eta -\alpha (\delta +\eta \mu -\rho )+y^{*}\left(S(\delta +\eta \mu )-\alpha (\mu +\beta \rho )+S\mu y^{*}\right)\right)}.\end{array}$$

In addition, assume that $(\sigma ,\rho ,\eta ,\delta ,\mu ,\alpha ,\beta )\in {\Psi }_2$; then, system (5) is characterized equivalently with the following planar map:

$$\left(\begin{array}{c}x\\ y\end{array}\right)\to \left(\begin{array}{c}\frac{x+\widehat{h}\left(\sigma +\frac{\rho xy}{\eta +y}\right)}{1+\widehat{h}\left(\delta +\mu y\right)}\hfill \\ \frac{y\left(1+\widehat{h}\alpha \right)}{1+\widehat{h}\left(x+\alpha \beta y\right)}\hfill \end{array}\right).$$

(29)

In order to discuss and analyze the normal form theory for Neimark–Sacker bifurcation for fixed point $\overline{X}=\left(x^{*},y^{*}\right)$ of (29), we suppose that $h_1$ represents a small perturbation in $\widehat{h}$. Then, the perturbed mapping for (29) can be described by the next map:

$$\left(\begin{array}{c}x\\ y\end{array}\right)\to \left(\begin{array}{c}\frac{x+(\widehat{h}+h_1)\left(\sigma +\frac{\rho xy}{\eta +y}\right)}{1+(\widehat{h}+h_1)\left(\delta +\mu y\right)}\\ \frac{y\left(1+(\widehat{h}+h_1)\alpha \right)}{1+(\widehat{h}+h_1)\left(x+\alpha \beta y\right)}\end{array}\right).$$

(30)

By taking $\overline{p}=x-x^{*},$$\overline{z}=y-x^{*}$ and $h=\widehat{h}+h_1,$ from (30) we can obtain the following mapping with an equilibrium point at $(0,0)$:

$$\left(\begin{array}{c}\overline{p}\\ \overline{z}\end{array}\right)\to M\left(\begin{array}{c}\overline{p}\\ \overline{z}\end{array}\right)+\left(\begin{array}{c}{F}_1(\overline{p},\overline{z})\\ F_2(\overline{p},\overline{z})\end{array}\right),$$

(31)

where

$$\begin{array}{c}\hfill M=\left(\begin{array}{cc}\frac{1+\frac{hy\rho }{y+\eta }}{1+h(\delta +y\mu )}& \frac{h\left(-\frac{xy\rho }{{(y+\eta )}^2}+\frac{x\rho }{y+\eta }\right)}{1+h(\delta +y\mu )}-\frac{h\mu \left(x+h\left(\frac{xy\rho }{y+\eta }+\sigma \right)\right)}{{(1+h(\delta +y\mu ))}^2}\\ -\frac{hy(1+h\alpha )}{{(1+h(x+y\alpha \beta ))}^2}& \frac{(1+hx)(1+h\alpha )}{{(1+h(x+y\alpha \beta ))}^2}\end{array}\right),\end{array}$$

(32)

and

$$\begin{array}{ccc}\hfill F_1(\overline{p},\overline{z})& =& \left(\frac{-h{(y^{*}+\eta )}^2\mu +h\left(\eta +h\delta \eta -h{y^{*}}^2\mu \right)\rho }{{(y^{*}+\eta )}^2{(1+h(\delta +x^{*}\mu ))}^2}\right)\overline{p}\overline{z}\hfill \\ & +& \left(\frac{h\left(\frac{2x^{*}{y}^{*}\rho }{{(y^{*}+\eta )}^3}-\frac{2x^{*}\rho }{{(y^{*}+\eta )}^2}\right)}{2(1+h(\delta +y^{*}\mu ))}-\frac{h^2\mu \left(-\frac{x^{*}{y}^{*}\rho }{{(y^{*}+\eta )}^2}+\frac{x^{*}\rho }{y^{*}+\eta }\right)}{{(1+h(\delta +y^{*}\mu ))}^2}+\frac{h^2{\mu }^2\left(x^{*}+h\left(\frac{x^{*}{y}^{*}\rho }{x^{*}+\eta }+\sigma \right)\right)}{{(1+h(\delta +y^{*}\mu ))}^3}\right){\overline{z}}^2\hfill \\ & +& \left(\frac{\frac{2h{y}^{*}\rho }{{(y^{*}+\eta )}^3}-\frac{2h\rho }{{(y^{*}+\eta )}^2}}{2(1+h(\delta +y^{*}\mu ))}-\frac{h\mu \left(-\frac{h{y}^{*}\rho }{{(y^{*}+\eta )}^2}+\frac{h\rho }{y^{*}+\eta }\right)}{{(1+h(\delta +y^{*}\mu ))}^2}+\frac{h^2{\mu }^2\left(1+\frac{h{y}^{*}\rho }{y^{*}+\eta }\right)}{{(1+h(\delta +y^{*}\mu ))}^3}\right)\overline{p}{\overline{z}}^2\hfill \\ & +& \left(\frac{h\left(-\frac{6x^{*}{y}^{*}\rho }{{(y^{*}+\eta )}^4}+\frac{6x^{*}\rho }{{(y^{*}+\eta )}^3}\right)}{6(1+h(\delta +y^{*}\mu ))}-\frac{h^2\mu \left(\frac{2x^{*}{y}^{*}\rho }{{(y^{*}+\eta )}^3}-\frac{2x^{*}\rho }{{(y^{*}+\eta )}^2}\right)}{2{(1+h(\delta +y^{*}\mu ))}^2}+\frac{h^3{\mu }^2\left(-\frac{x^{*}{y}^{*}\rho }{{(y+\eta )}^2}+\frac{x^{*}\rho }{y^{*}+\eta }\right)}{{(1+h(\delta +y^{*}\mu ))}^3}\right){\overline{z}}^3\hfill \\ & -& \left(\frac{h^3{\mu }^3\left(x^{*}+h\left(\frac{x^{*}{y}^{*}\rho }{y^{*}+\eta }+\sigma \right)\right)}{{(1+h(\delta +y^{*}\mu ))}^4}\right){\overline{z}}^3,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill F_2(\overline{p},\overline{z})& =& \left(\frac{h^2{y}^{*}(1+h\alpha )}{{(1+h{x}^{*}+h{y}^{*}\alpha \beta )}^3}\right){\overline{p}}^2+\left(\frac{2h^2{y}^{*}\alpha (1+h\alpha )\beta }{{(1+h(x^{*}+y^{*}\alpha \beta ))}^3}-\frac{h(1+h\alpha )}{{(1+h(x^{*}+y^{*}\alpha \beta ))}^2}\right)\overline{p}\overline{z}\hfill \\ & -& \left(\frac{h(1+h{x}^{*})\alpha (1+h\alpha )\beta }{{(1+h(x^{*}+y^{*}\alpha \beta ))}^3}\right){\overline{z}}^2-\left(-\frac{h^3{y}^{*}(1+h\alpha )}{{(1+h(x^{*}+y^{*}\alpha \beta ))}^4}\right){\overline{p}}^3\hfill \\ & +& \left(\frac{h^2(1+h\alpha )(1+h(x^{*}-2y^{*}\alpha \beta ))}{{(1+h(x^{*}+y^{*}\alpha \beta ))}^4}\right){\overline{p}}^2\overline{z}+\left(\frac{h^2\alpha (1+h\alpha )\beta (2+2h{x}^{*}-h{y}^{*}\alpha \beta )}{{(1+h(x^{*}+y^{*}\alpha \beta ))}^4}\right)\overline{p}{\overline{z}}^2\hfill \\ & +& \left(\frac{h^2(1+h{x}^{*}){\alpha }^2(1+h\alpha ){\beta }^2}{{(1+h(x^{*}+y^{*}\alpha \beta ))}^4}\right){\overline{z}}^3.\hfill \end{array}$$

The characteristic equation $\mathbb{S}\left(\varphi \right)=0$ generated by the Jacobian matrix of (31) about $(0,0)$ can be written as

$${\varphi }^2-\widehat{T}\left(h_1\right)\varphi +\widehat{D}\left(h_1\right)=0,$$

(33)

where

$$\begin{array}{ccc}\hfill \widehat{D}\left(h_1\right)& =& \frac{\eta +(\widehat{h}+h_1)\alpha \eta +y^{*}\left(1+(\widehat{h}+h_1)(\alpha +S\eta -\alpha \beta \eta +(\widehat{h}+h_1)S\delta \eta +\rho +(\widehat{h}+h_1)\alpha \rho )\right)}{(1+(\widehat{h}+h_1)\alpha )\left(\eta +y^{*}\right)\left(1+(\widehat{h}+h_1)\delta +(\widehat{h}+h_1)\mu y^{*}\right)}\hfill \\ & +& \frac{y^{*}\left((\widehat{h}+h_1)y^{*}\left(S(1+(\widehat{h}+h_1)\delta +(\widehat{h}+h_1)\eta \mu )-\alpha \beta (1+(\widehat{h}+h_1)\rho )+(\widehat{h}+h_1)S\mu y^{*}\right)\right)}{(1+(\widehat{h}+h_1)\alpha )\left(\eta +y^{*}\right)\left(1+(\widehat{h}+h_1)\delta +(\widehat{h}+h_1)\mu y^{*}\right)},\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \widehat{T}\left(h_1\right)& =& 1-\frac{(\widehat{h}+h_1)\alpha \beta y^{*}}{1+(\widehat{h}+h_1)\alpha }+\frac{1+\frac{(\widehat{h}+h_1)\rho y^{*}}{\eta +y^{*}}}{1+(\widehat{h}+h_1)\delta +(\widehat{h}+h_1)\mu y^{*}}.\hfill \end{array}$$

Consider that $(\sigma ,\rho ,\eta ,\delta ,\mu ,\alpha ,\beta )\in {\Psi }_2$; at that point, the complex solutions for (33) are calculated as follows:

$$\begin{array}{c}\hfill {\varphi }_1=\frac{\widehat{T}\left(h_1\right)-i\sqrt{4\widehat{D}\left(h_1\right)-{\widehat{T}}^2\left(h_1\right)}}{2}\end{array}$$

and

$$\begin{array}{c}\hfill {\varphi }_2=\frac{\widehat{T}\left(h_1\right)+i\sqrt{4\widehat{D}\left(h_1\right)-{\widehat{T}}^2\left(h_1\right)}}{2}.\end{array}$$

Moreover, we have

$${\left(\frac{d\mid {\varphi }_1\mid }{d{h}_1}\right)}_{h_1=0}\ne 0$$

(34)

because $|\widehat{T}\left(0\right)|<2$ as $(\sigma ,\rho ,\eta ,\delta ,\mu ,\alpha ,\beta )\in {\Psi }_2.$ Moreover, a simple computation yields that

$$\begin{array}{ccc}\hfill \widehat{T}\left(0\right)& =& 1-\frac{\widehat{h}\alpha \beta y^{*}}{1+\widehat{h}\alpha }+\frac{1+\frac{\widehat{h}\rho y^{*}}{\eta +y^{*}}}{1+\widehat{h}\delta +\widehat{h}\mu y^{*}},\hfill \end{array}$$

and we suppose that $\widehat{T}\left(0\right)\ne 0$ and $\widehat{T}\left(0\right)\ne 1$, that is,

$$\begin{array}{c}\hfill 1+\frac{1+\frac{\widehat{h}\rho y^{*}}{\eta +y^{*}}}{1+\widehat{h}\delta +\widehat{h}\mu y^{*}}\ne \frac{\widehat{h}\alpha \beta y^{*}}{1+\widehat{h}\alpha },\frac{1+\frac{\widehat{h}\rho y^{*}}{\eta +y^{*}}}{1+\widehat{h}\delta +\widehat{h}\mu y^{*}}\ne \frac{\widehat{h}\alpha \beta y^{*}}{1+\widehat{h}\alpha }.\end{array}$$

(35)

Suppose that (35) holds and $(\sigma ,\rho ,\eta ,\delta ,\mu ,\alpha ,\beta )\in {\Psi }_2.$ Then, it follows that $\widehat{T}\left(0\right)\ne \pm 2,0,-1$, that is, ${\varphi }_1^m,{\varphi }_2^m\ne 1$ for every $m\in \{1,2,3,4\}$ about $h_1=0$. Consequently, both solutions of (33) do not lie inside the intersection of the unit circle with the coordinate axes when $h_1=0.$ In the same way, we assume that $\lambda =\frac{\widehat{T}\left(0\right)}{2}$, $\omega =\frac{1}{2}\sqrt{4\widehat{D}\left(0\right)-{\widehat{T}}^2\left(0\right)}.$ Formerly, to change (31) into normal form, we used the following similarity transformation:

$$\left(\begin{array}{c}\overline{p}\\ \overline{z}\end{array}\right)=\left(\begin{array}{cc}{l}_{12}& 0\\ \lambda -l_{11}& -\omega \end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right).$$

(36)

By using (36), we obtain the next typical form for (31):

$$\left(\begin{array}{c}u\\ v\end{array}\right)\to \left(\begin{array}{cc}\lambda & -\omega \\ \omega & \lambda \end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)+\left(\begin{array}{c}\tilde{F}(u,v)\\ \tilde{G}(u,v)\end{array}\right).$$

(37)

Moreover, we have

$$\begin{array}{ccc}\hfill \tilde{F}(u,v)& =& \left({\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^2{l}_{15}+\left(\left(\lambda -l_{11}\right)u-\omega v\right)l_{13}\right)u\hfill \\ & +& \frac{{\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^3{l}_{16}+{\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^2{l}_{14}}{l_{12}}+O\left(\left(\right|u|+{\left|v\right|)}^4\right),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill \tilde{G}(u,v)& =& \left(\frac{\left(\lambda -l_{11}\right)\left({\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^2{l}_{15}+\left(\left(\lambda -l_{11}\right)u-\omega v\right)l_{13}\right)}{l_{12}\omega }\right)l_{12}u-\frac{l_{26}{l_{12}}^3{u}^3}{\omega }\hfill \\ & -& \frac{\left(\left(\lambda -l_{11}\right)u-\omega v\right)l_{27}{l_{12}}^2{u}^2}{\omega }-\left(\frac{{\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^2{l}_{28}+\left(\left(\lambda -l_{11}\right)u-\omega v\right)l_{24}}{\omega }\right)l_{12}u\hfill \\ & +& \frac{\left(\lambda -l_{11}\right)\left({\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^3{l}_{16}+{\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^2{l}_{14}\right)}{l_{12}\omega }\hfill \\ & -& \frac{{\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^3{l}_{29}+l_{12}^2{u}^2{l}_{23}+{\left(\left(\lambda -l_{11}\right)u-\omega v\right)}^2{l}_{25}}{\omega }+O\left(\left(\right|u|+{\left|v\right|)}^4\right).\hfill \end{array}$$

where $l_{11},l_{12},l_{13}$ and $l_{14}$ are respective elements of $M.$ In addition, $l_{1j}$ for $j=3,4,5,6$ and $l_{2j}$ for $j=3,4,\dots ,9$ are coefficients from expressions $F_1(\overline{p},\overline{z})$ and $F_2(\overline{p},\overline{z})$, respectively. Now, we describe the next non-zero real numbers:

$$ϝ={\left(\left[-Re\left(\frac{(1-2{\varphi }_1){\varphi }_2^2}{1-{\varphi }_1}{\xi }_{20}{\xi }_{11}\right)-\frac{1}{2}|{\xi }_{11}{|}^2-{\left|{\xi }_{02}\right|}^2+Re\left({\varphi }_2{\xi }_{21}\right)\right]\right)}_{h_1=0},$$

where

$$\begin{array}{ccc}& {\xi }_{20}=& \frac{1}{8}\left(2\left(\lambda -l_{11}\right)l_{13}-\frac{2{\omega }^2{l}_{14}}{l_{12}}+\frac{2\left(\lambda -l_{11}\right){}^2{l}_{14}}{l_{12}}\right)\hfill \\ & +& \frac{1}{4}\left(-\frac{2\left(\lambda -l_{11}\right){}^2{l}_{14}}{l_{12}}+l_{12}\left(-\frac{\left(\lambda -l_{11}\right)l_{13}}{l_{12}}+l_{24}\right)+2\left(\lambda -l_{11}\right)l_{25}\right)\hfill \\ & +& \frac{i}{4}\left(-\frac{\omega \left(\lambda -l_{11}\right)l_{14}}{l_{12}}+\frac{\left(\lambda -l_{11}\right){}^3{l}_{14}}{\omega l_{12}}-\left(-\omega l_{13}-\frac{2\omega \left(\lambda -l_{11}\right)l_{14}}{l_{12}}\right)\right)\hfill \\ & +& l_{12}\frac{i}{4}\left(\left(\frac{\left(\lambda -l_{11}\right){}^2{l}_{13}}{\omega l_{12}}-\frac{\left(\lambda -l_{11}\right)l_{24}}{\omega }\right)+\omega l_{25}-\frac{l_{12}^2{l}_{23}+\left(\lambda -l_{11}\right){}^2{l}_{25}}{\omega }\right)\hfill \end{array}$$

$$\begin{array}{ccc}& {\xi }_{11}=& \frac{1}{2}\left(\left(\lambda -l_{11}\right)l_{13}+\frac{{\omega }^2{l}_{14}}{l_{12}}+\frac{\left(\lambda -l_{11}\right){}^2{l}_{14}}{l_{12}}+i\left(\frac{\omega \left(\lambda -l_{11}\right)l_{14}}{l_{12}}+\frac{\left(\lambda -l_{11}\right){}^3{l}_{14}}{\omega l_{12}}\right)\right)\hfill \\ & +& \frac{i}{2}\left(l_{12}\left(\frac{\left(\lambda -l_{11}\right){}^2{l}_{13}}{\omega l_{12}}-\frac{\left(\lambda -l_{11}\right)l_{24}}{\omega }\right)-\omega l_{25}-\frac{l_{12}^2{l}_{23}+\left(\lambda -l_{11}\right){}^2{l}_{25}}{\omega }\right),\hfill \end{array}$$

$$\begin{array}{ccc}& {\xi }_{02}=& \frac{1}{4}\left(\left(\lambda -l_{11}\right)l_{13}-\frac{{\omega }^2{l}_{14}}{l_{12}}+\frac{\left(\lambda -l_{11}\right){}^2{l}_{14}}{l_{12}}+\frac{2\left(\lambda -l_{11}\right){}^2{l}_{14}}{l_{12}}+l_{12}\left(\frac{\left(\lambda -l_{11}\right)l_{13}}{l_{12}}-l_{24}\right)\right)\hfill \\ & +& \frac{1}{4}\left(\left(2\left(\lambda -l_{11}\right)l_{25}\right)+i\left(-\frac{\omega \left(\lambda -l_{11}\right)l_{14}}{l_{12}}+\frac{\left(\lambda -l_{11}\right){}^3{l}_{14}}{\omega l_{12}}+\omega l_{25}-\omega l_{13}\right)\right)\hfill \\ & -& \frac{i}{4}\left(\frac{2\omega \left(\lambda -l_{11}\right)l_{14}}{l_{12}}-l_{12}\left(\frac{\left(\lambda -l_{11}\right){}^2{l}_{13}}{\omega l_{12}}-\frac{\left(\lambda -l_{11}\right)l_{24}}{\omega }\right)+\frac{l_{12}^2{l}_{23}+\left(\lambda -l_{11}\right){}^2{l}_{25}}{\omega }\right),\hfill \end{array}$$

and

$$\begin{array}{ccc}& {\xi }_{21}=& \frac{1}{8}\left({\omega }^2{l}_{15}+3\left(\lambda -l_{11}\right){}^2{l}_{15}+l_{12}^2{l}_{27}-l_{12}\left(\frac{2\left(\lambda -l_{11}\right){}^2{l}_{15}}{l_{12}}-2\left(\lambda -l_{11}\right)l_{28}\right)\right)\hfill \\ & +& \frac{1}{8}\left(3{\omega }^2{l}_{29}+3\left(\lambda -l_{11}\right){}^2{l}_{29}+i\left(2\omega \left(\lambda -l_{11}\right)l_{15}+\frac{3{\omega }^3{l}_{16}}{l_{12}}+\frac{6\omega \left(\lambda -l_{11}\right){}^2{l}_{16}}{l_{12}}\right)\right)\hfill \\ & +& \frac{i}{8}\left(\frac{3\left(\lambda -l_{11}\right){}^4{l}_{16}}{\omega l_{12}}-\frac{3l_{12}^3{l}_{26}}{\omega }-\frac{3\left(\lambda -l_{11}\right)l_{12}^2{l}_{27}}{\omega }+l_{12}\left(\frac{\omega \left(\lambda -l_{11}\right)l_{15}}{l_{12}}-\omega l_{28}\right)\right)\hfill \\ & +& \frac{i}{8}\left(3l_{12}\left(\frac{\left(\lambda -l_{11}\right){}^3{l}_{15}}{\omega l_{12}}-\frac{\left(\lambda -l_{11}\right){}^2{l}_{28}}{\omega }\right)-3\omega \left(\lambda -l_{11}\right)l_{29}-\frac{3\left(\lambda -l_{11}\right){}^3{l}_{29}}{\omega }\right).\hfill \end{array}$$

Hence, we have the following theorem.

Theorem 6.

Assume that (35) holds true and $ϝ\ne 0.$ Then, the positive fixed point

$$\overline{X}=\left(x^{*},y^{*}\right)$$

of system (5) experiences Neimark–Sacker bifurcation whenever h changes in the least neighborhood of

$$\begin{array}{c}\hfill \widehat{h}=\frac{\delta \eta +y^{*}\left(\delta +\eta (-S+\alpha \beta +\mu )-\rho +(-S+\alpha \beta +\mu )y^{*}\right)}{-\alpha \delta \eta +y^{*}\left(S\delta \eta -\alpha (\delta +\eta \mu -\rho )+y^{*}\left(S(\delta +\eta \mu )-\alpha (\mu +\beta \rho )+S\mu y^{*}\right)\right)}.\end{array}$$

In addition, if $ϝ<0,(ϝ>0)$, respectively, then an attracting or repelling invariant closed curve bifurcates from the equilibrium point for $h>\widehat{h}(h<\widehat{h})$, respectively.

6. Control of Neimark–Sacker Bifurcation

The study of chaos theory and bifurcation control is a multidisciplinary area of mathematics research that concentrates on basic designs and extremely complex categorical laws of primary conditions in any dynamical systems that are supposed to have entirely arbitrary statuses of disorder and inconsistency. Generally, the leading standard of disorder defines how a minor variation in any state of a nonlinear dynamical system can result in significant changes in an advanced state (the implication being that a complex dependency on primary conditions is close) [48]. Every disordered attractor encloses a countless amount of periodic and unstable orbits. Chaotic behaviour at any time is a gesture where the state system moves in the vicinity of any of these regions for a time and then drops to a closer periodic and unstable orbit, where it hangs for a degree of time, etcetera [48]. Chaos control stabilizes any of these irregular periodic orbits by the worth of small structure perturbations. Hence, we use a simple chaos control method for system (5). Furthermore, many well-known techniques have been developed in previous decades to control chaos in any discrete dynamical system. We refer readers to [48,49,50] for additional details connected to these methods. Here, we implement a generalized hybrid control technique to control the Neimark–Sacker bifurcation (seecite [51,52,53]). The generalized hybrid control method [48] is centred on parameter perturbation and a state feedback control technique. By applying generalized hybrid control methodology (with control parameter $b\in (0,1]$) to system (5), we obtain

$$\left\{\begin{array}{c}{x}_{n+1}=sinb\left(\frac{x_n+h\left(\sigma +\frac{\rho x_n{y}_n}{\eta +y_n}\right)}{1+h\left(\delta +\mu y_n\right)}\right)+(1-sinb)x_n,\hfill \\ y_{n+1}=sinb\left(\frac{y_n\left(1+h\alpha \right)}{1+h\left(x_n+\alpha \beta y_n\right)}\right)+(1-sinb)y_n.\hfill \end{array}\right.$$

(38)

Then, system (38) is controllable provided that its fixed point $(x^{*},y^{*})$ is locally asymptotically stable. Additionally, the Jacobian matrix for system (38) at its positive fixed point $(x^{*},y^{*})$ is calculated as follows:

$$J_c=\left[\begin{array}{cc}1-sinb+\frac{sinb\left(1+\frac{h\rho y^{*}}{\eta +y^{*}}\right)}{1+h\delta +h\mu y^{*}}& -\frac{hsinb\left(h\mu \sigma {\left(\eta +y^{*}\right)}^2+x^{*}\left(\eta (\eta \mu -(1+h\delta )\rho )+\mu y^{*}\left(2\eta +y^{*}+h\rho y^{*}\right)\right)\right)}{{\left(\eta +y^{*}\right)}^2{\left(1+h\delta +h\mu y^{*}\right)}^2}\\ -\frac{hsinb(1+h\alpha )y^{*}}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}& 1-sinb+\frac{sinb(1+h\alpha )\left(1+h{x}^{*}\right)}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}\end{array}\right]$$

(39)

where

$$Trace\left[J_c\right]=2-2sinb+\frac{sinb(1+h\alpha )\left(1+h{x}^{*}\right)}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}+\frac{sinb\left(1+\frac{h\rho y^{*}}{\eta +y^{*}}\right)}{1+h\delta +h\mu y^{*}}$$

(40)

and

$$\left\{\begin{array}{c}Det\left[J_c\right]=\frac{\left(1-h\tilde{A}\delta -h\tilde{A}\mu y^{*}\right)\left(1+hsinb\alpha -h\tilde{A}\left(x^{*}+\alpha \beta y^{*}\right)\left(2+h{x}^{*}+h\alpha \beta y^{*}\right)\right)}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2\left(1+h\delta +h\mu y^{*}\right)}\hfill \\ -\frac{hsinb(1+h\alpha )x^{*}\left((1-h\tilde{A}\delta )\eta +y^{*}\left(1-h\tilde{A}(\delta +\eta \mu )+hsinb\rho -h\tilde{A}\mu y^{*}\right)\right)}{\left(\eta +y^{*}\right){\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2\left(1+h\delta +h\mu y^{*}\right)}\hfill \\ +hsinb{y}^{*}\left(\frac{S(sinb+hsinb\alpha )}{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}+\frac{\rho \left(1+sinb\left(\frac{1+h\alpha }{{\left(1+h{x}^{*}+h\alpha \beta y^{*}\right)}^2}-1\right)\right)}{\left(\eta +y^{*}\right)\left(1+h\delta +h\mu y^{*}\right)}\right),\hfill \end{array}\right.$$

(41)

with $sinb-1=\tilde{A}.$ Moreover, we have the following result:

Theorem 7.

Assume the fixed point $(x^{*},y^{*})$ of system (38). Then, $(x^{*},y^{*})$ is locally asymptotically stable ⟺ we have

$$\begin{array}{c}\hfill \left|Trace\left[J_c\right]\right|<1+Det\left[J_c\right]<2,\end{array}$$

where $Trace\left[J_c\right]$ and $Det\left[J_c\right]$ are as defined in (40) and (41), respectively.

7. Numerical Simulations

First, assume that $\alpha =1.636,\beta =2\times {10}^{-3},\delta =0.3743,\eta =20.19,\mu =0.00311,$$\rho =1.131,\sigma =0.1181$, and $h\in [0,1).$ Then, from system (5) we have $(x^{*},y^{*})=(1.61954,$ 8.2317). Moreover, in this case $x_0=1.61954$ and $y_0=8.2317$ are initial conditions. The graphical behavior of both variables is shown in Figure 5. It can be seen that $x_n$ and $y_n$ undergo Neimark–Sacker bifurcation at unique positive fixed points $(x^{*},y^{*})=(1.61954,8.2317).$ In addition, in Figure 5c the maximum Lyapunov exponents are represented. To understand the consistency between bifurcation diagrams and Lyapunov exponents, in Figure 5a, one can see that bifurcation in the first variable occurs when h passes through $h=0.5$, and the Lyapunov exponent changes from negative to positive as h crosses the horizontal line at $h=0.5.$

In this case, we have $ϝ=-0.00065841<0,$ which verifies Theorem 6. Moreover, we have $\mathbb{C}(-1)=3.2136789>0,$ which validates Theorem 5. By varying the stepsize, h, in $[0,1)$ phase portraits for system (5) can be obtained, as shown in Figure 6. Hence, we can observe that system (5) experiences Neimark–Sacker bifurcation when the parameter h certainly passes through $h=0.4917267952$ (see Figure 6c). Second, to discuss the feasibility of the designed control technique, we take $\alpha =1.636,\beta =2\times {10}^{-3},\delta =0.3743,\eta =20.19,$$\mu =0.00311,\rho =1.131,\sigma =0.1181$ and $h\in [0,1).$

Then, from system (38), we have $(x^{*},y^{*})=(1.61954,8.2317).$ Moreover, by taking b as the control parameter, it can be observed from Figure 7 that the Neimark–Sacker bifurcation at unique positive fixed point $(x^{*},y^{*})$ is effectively controlled for a large range of the control parameter $b.$ Additionally, the MLE for controlled system (38) is provided in Figure 7c. From Figure 7a, it can be seen that bifurcation in the first variable occurs when b lies between $0<b<0.1$. Moreover, the controlled system is stable in $0.1<b<0.69$, and the Lyapunov exponent changes from positive to negative at the point $b=0.1\left(approx\right)$ and negative to positive as b crosses the horizontal line at $b=0.69,$ which shows the consistency of controlled plots with corresponding MLE.

Finally, to discuss the dynamics of system (3), we take $\alpha =1.636,\beta =2\times {10}^{-3},$$\delta =0.3743,\eta =20.19,\mu =0.00311,\rho =11.131,$ and $\sigma =0.1181.$ Then, from system (3) we have $(x^{*},y^{*})=(1.519,0.823).$ Consequently, we obtain a stable system (3) for these values. In addition, the qualitative behaviour of system (3) is shown in Figure 8. On the other side, by taking $\sigma$ as bifurcation parameter and taking $\sigma =0.001181,$ it can be observed from Figure 9 that the system (3) experiences Hopf bifurcation at positive fixed point $(x^{*},y^{*}).$ Additionally, the plots of both variables for system (3) are provided in Figure Figure 9a,b. In this case, the value of FLE is calculated as approximately $L\left({\sigma }_0\right)=-0.00645890013.$

8. Conclusions

Here, we have considered an immunogenic tumor model for discretization and qualitative study. The original model (3) was presented and studied by Kuznetsov et al. [28] in its continuous form. Moreover, they studied the dynamics of (3) in its continuous form as well. This work is focused on the study of the consistent counterpart of (3) and comparing the dynamics of model (3) with its discrete-time counterpart, which has not been performed previously for this model. Hence, we first converted the system (3) into its discrete form using a consistency-preserving discretization method. For this purpose, a nonstandard finite difference method was applied to obtain a discrete counterpart of the particular model (3). Our examination exposes that the continuous system (3) undergoes Hopf bifurcation about its positive fixed point $\sigma$, which is considered a bifurcation parameter, and passes over a critical value ${\sigma }_0=x^{*}\left(\alpha -x^{*}-2y^{*}\alpha \beta \right)$, such as $\rho >\delta .$ Furthermore, the first Lyapunov exponent is calculated in the closed form, specified as follows:

$$L\left({\sigma }_0\right)=L_1+L_2,$$

where

$$\left\{\begin{array}{c}{L}_1=\frac{1}{16}\left(\frac{2\rho y^{*}}{{\left(\eta +y^{*}\right)}^3}-\frac{2\rho }{{\left(\eta +y^{*}\right)}^2}\right)\hfill \\ L_2=\frac{-2\alpha \beta -2\alpha \beta \left(\frac{2\rho x^{*}{y}^{*}}{{\left(\eta +y^{*}\right)}^3}-\frac{2\rho x^{*}}{{\left(\eta +y^{*}\right)}^2}\right)+\left(\frac{2\rho x^{*}{y}^{*}}{{\left(\eta +y^{*}\right)}^3}-\frac{2\rho x^{*}}{{\left(\eta +y^{*}\right)}^2}\right)\left(-\mu -\frac{\rho y^{*}}{{\left(\eta +y^{*}\right)}^2}+\frac{\rho }{\eta +y^{*}}\right)}{16\sqrt{-{\left(x^{*}\right)}^2-{\alpha }^2{\left(1-2\beta y^{*}\right)}^2+x^{*}\left(2\alpha +y^{*}\left(-4\alpha \beta -\mu +\frac{\eta \rho }{{\left(\eta +y^{*}\right)}^2}\right)\right)}}.\hfill \end{array}\right.$$

On the other hand, when h is taken as a bifurcation parameter, the discrete-time version, which is obtained using a nonstandard finite difference scheme, experiences Neimark–Sacker bifurcation about its positive fixed point $(x^{*},y^{*})$ whenever h passes through

$$\begin{array}{c}\hfill \widehat{h}=\frac{\delta \eta +y^{*}\left(\delta +\eta (-S+\alpha \beta +\mu )-\rho +(-S+\alpha \beta +\mu )y^{*}\right)}{-\alpha \delta \eta +y^{*}\left(S\delta \eta -\alpha (\delta +\eta \mu -\rho )+y^{*}\left(S(\delta +\eta \mu )-\alpha (\mu +\beta \rho )+S\mu y^{*}\right)\right)}.\end{array}$$

The conditions for the presence of Neimark–Sacker bifurcation are specified in Theorem 6. Mathematical simulation exposes that our discretization is bifurcation-conserving and equal with lesser step size value; the first Lyapunov exponents are approximately the same in both cases, that is, $ϝ\approx L\left({\sigma }_0\right)=-0.006832.$ From Theorem 5, it can be seen that there is no chance of period-doubling bifurcation in our discrete-time system, which shows the consistency of the discretizing technique used here. Moreover, to analyze the wide-ranging and rich dynamics of another discrete-time counterpart of the immunogenic tumors model, it is possible to use the Euler method or piecewise constant arguments with system (3). Using the Euler method or piecewise constant arguments, it is possible to discuss other types of bifurcations and chaos control. We refer readers to [44,45,46,47,48] and the references therein for further consideration. We anticipate that analysis of system (3) using the Euler method and bifurcation analysis with chaos control for the obtained system will be our future tasks.