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

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.


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]: where 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, dC dt ≈ 0), which yields the following system of equations (see [28]): where p = f K, m = Kk 3 , n = Kk 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, τ = nT 0 t. Formally, the model can be re-articulated as where x = E E 0 , y = T T 0 and parameter estimates for system (3) are the following (see Table 1): Table 1. Parameter estimates from [28].

System Parameters
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 cancerimmune 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]): where h ∈ [0, 1) is the step size for discretization. Furthermore, system (4) can be transformed into the following form: 1+h(x n +αβy n ) .
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.

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 Hence, per the Routh-Hurwitz stability criterion, (x * , y * ) is sink if and only if (3) experiences Hopf bifurcation if the parameters lie in the following curve: 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: where a ∈ 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: 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) φ(a 0 ) = 0 and ϕ(a 0 ) = ϕ 0 = 0. Then, there exists a conjugate pair of complex eigenvalues of V(x * , y * ) in the condition of non-hyperbolicity.
da | a=a 0 = T = 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(a 0 ) = 0, which is known as the first Lyapunov exponent (FLE) and is defined as follows: ∂x∂y computed at (x, y) = (x * , y * ) and 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 σ + ασ(2βy * −1) then, it is easy to see that the eigenvalues of the Jacobian matrix V(x * , y * ) are of the form For the transversality condition, we can see that In order to shift the fixed point of system (3) to origin (0, 0), we consider the following translations: Moreover, by implementing this transformation on system (3), we obtain Application of Taylor series expansion on (u, v) = (0, 0) provides the following system: Next, we want to convert matrix B into its canonical form. For this, the following similarity transformation is considered: From (8) and (9), it follows that where Then, the first Lyapunov exponent for system (3) is computed as follows: Finally, we have the following theorem: In addition, the fixed point is stable whenever σ > σ 0 (respectively, σ < σ 0 ) and unstable for σ < σ 0 (respectively, σ > σ 0 ) for T < 0 and T > 0, respectively.
The plot of FLE is depicted in Figure 3.

Boundedness of Solutions
From the second equation of system (5) it follows that y n+1 ≤ y n (1 + hα) 1 + hαβy n .
Consequently, by solving (11) and then applying the limit, we obtain lim n−→∞ for all n ≥ 0. In the same way, from the first equation of system (5) we obtain because x n + h σ + ρx n y n η + y n ≤ x n + h σ + ρx n 1 + ηβ yields y n ≤ 1 β for all n ≥ 0. Then, from (13) we have Hence, we can obtain the upper bound for x n as for all n ≥ 0. Finally, we have the following theorem concerning the boundedness of all solutions of (5).

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: 1+h(x * +αβy * ) .
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 In addition, the characteristic polynomial C(φ) of J (x,y) is 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]).

Moreover, we have
and Hence, for the study of the linearized stability of system (5) about (x * , y * ), we have the following proposition (see Lemma 1).
We now have the following theorem for the possible validation of Remark 1.

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 ∈ (0, 1]) to system (5), we obtain x n +h σ+ ρxn yn η+yn 1+h(δ+µy n ) 1+h(x n +αβy n ) + (1 − sin b)y n .
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: and with sin b − 1 =Ã. 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 where Trace[J c ] and Det[J c ] are as defined in (40) and (41), respectively.

Numerical Simulations
First, assume that α = 1.636, β = 2 × 10 −3 , δ = 0.3743, η = 20.19, µ = 0.00311, ρ = 1.131, σ = 0.1181, and h ∈ [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 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 α = 1.636, β = 2 × 10 −3 , δ = 0.3743, η = 20.19, µ = 0.00311, ρ = 1.131, σ = 0.1181 and h ∈ [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(approx) and negative to positive as b crosses the horizontal line at b = 0.69, which shows the consistency of controlled plots with corresponding MLE.

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 σ, which is considered a bifurcation parameter, and passes over a critical value σ 0 = x * (α − x * − 2y * αβ), such as ρ > δ. Furthermore, the first Lyapunov exponent is calculated in the closed form, specified as follows: 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 througĥ h = δη + y * (δ + η(−S + αβ + µ) − ρ + (−S + αβ + µ)y * ) −αδη + y * (Sδη − α(δ + ηµ − ρ) + y * (S(δ + ηµ) − α(µ + βρ) + Sµy * )) .
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, ≈ L(σ 0 ) = −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.