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.

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, 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, , , , and are non-negative kinetic real numbers; and refer to the degrees of binding of TC from EC and objectivity of TC to EC without injuring cells, is the degree at which EC to TC connections conclusively program TC for lysis, and 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 and 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, ), which yields the following system of equations (see [28]): where and are dimensional parameters. In addition, for better study of the dynamics of model (2) it is necessary to non-dimensionalize the system (2).

Figure 1

Data flow diagram for our system.

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, and , respectively, as proposed from the tests discussed in the previous section: = = cells (see [28]). Time is scaled comparative to the degree of cancer cell deactivation, that is, Formally, the model can be re-articulated as where and parameter estimates for system (3) are the following (see Table 1):

Table 1

Parameter estimates from [28].

System Parameters	Dimensional Parameters	Estimated Values
σ	σ=snE_0T_0	0.001181
ρ	ρ=pnT_0	11.131
η	η=gT_0	20.19
μ	μ=mn	0.00311
δ	δ=dnT_0	0.3743
α	α=anT_0	1.636
β	β=bnT_0	2.0×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]): where is the step size for discretization. Furthermore, system (4) can be transformed into the following form: 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 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 at is provided by where and Hence, per the Routh–Hurwitz stability criterion, is sink if and only if and and source if and only if and . Moreover, system (3) experiences Hopf bifurcation if the parameters lie in the following curve: Furthermore, Figure 2 shows the topological classification of the fixed point . 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 is the bifurcation parameter. Let be the Jacobian matrix of (6) computed at equilibrium point . Moreover, the eigenvalues of (6) computed at any equilibrium point are of the following form: Furthermore, we assume that there is a particular value of the bifurcation parameter a, say , for which the following conditions hold true (see [40]): where with computed at and

Figure 2

Topological classification of the one and only fixed point of system (3) for and with initial conditions and .

and . Then, there exists a conjugate pair of complex eigenvalues of in the condition of non-hyperbolicity.

, which is known as a transversality condition, that is, the eigenvalues of cross the imaginary axis with non-zero speed [40].

There exists a discriminatory quantity , which is known as the first Lyapunov exponent (FLE) and is defined as follows:

([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

It can be seen that periodic solutions are stable (respectively, unstable) if the fixed point is unstable (respectively, stable) on the side of 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 ; then, it is easy to see that the eigenvalues of the Jacobian matrix are of the form and Next, provides . At , we have For the transversality condition, we can see that In order to shift the fixed point of system (3) to origin , we consider the following translations: Moreover, by implementing this transformation on system (3), we obtain Application of Taylor series expansion on provides the following system: where 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 and Then, the first Lyapunov exponent for system (3) is computed as follows: where Finally, we have the following theorem:

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

The plot of FLE is depicted in Figure 3.

Figure 3

FLE for system (3) for and with initial conditions and .

3. Boundedness of Solutions

From the second equation of system (5) it follows that Consequently, by solving (11) and then applying the limit, we obtain for all In the same way, from the first equation of system (5) we obtain because yields for all Then, from (13) we have Hence, we can obtain the upper bound for as for all Finally, we have the following theorem concerning the boundedness of all solutions of (5).

Assume that

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: Furthermore, from (15) we obtain the following pair: Moreover, we have two fixed points from (15), namely, and where is solution of the following equation: with In addition, and if we have Hence, using Descartes’ rule of signs (see [41]), we have the following result:

Assume that there exists a unique positive constant solution

In order to discuss the stability of system (5) about these fixed points, we compute the Jacobian matrix of system (5) about each of its fixed points . Moreover, is In addition, the characteristic polynomial of is where T and D represent the trace and determinants of , respectively. The following lemma describes the condition equivalent to the Jury conditions for the local stability of fixed points (see [42]).

([42]). Let

(i)

(ii)

(iii)

(iv)

As

First, we study the stability conditions for (5) about the fixed point The matrix evaluated at is provided by Furthermore, the eigenvalues of are and with for all parametric values. Hence, we have the following result related to the dynamics of (5) about :

The boundary equilibrium

Next, our task is to explore of the local stability of system (5) about the point . Let be the Jacobian matrix of system (5) about the fixed point ; then, has the following mathematical form:

Moreover, from we can obtain the following characteristic polynomial: where and By considering (19) and taking it follows that for Moreover, we have and Hence, for the study of the linearized stability of system (5) about we have the following proposition (see Lemma 1).

Let ( then,

The fixed point we have

The fixed point and

Let (

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

Let ( Then,

Assume (18) and Additionally, let be the positive fixed point of system (5) and Then, from we have Finally, under condition (26), we have if and only if where, are defined in (24). Hence, under condition (26), the above inequality (27) is true for every choice of and This completes the proof. □

5. Neimark–Sacker Bifurcation

This section is related to the bifurcation analysis of system (5) about Moreover, all conditions for the existence and positivity of 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 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 . Assume that then, one can see from Proposition 2 that roots of (21) are complex and satisfy if and only if and where Furthermore, under the suppositions that we study the following set: Then, the positive fixed point of system (5) experiences Neimark–Sacker bifurcation such that h is taken as the bifurcation parameter and varies slightly in the neighborhood of which is provided as In addition, assume that ; then, system (5) is characterized equivalently with the following planar map: In order to discuss and analyze the normal form theory for Neimark–Sacker bifurcation for fixed point of (29), we suppose that represents a small perturbation in . Then, the perturbed mapping for (29) can be described by the next map: By taking and from (30) we can obtain the following mapping with an equilibrium point at : where and The characteristic equation generated by the Jacobian matrix of (31) about can be written as where Consider that ; at that point, the complex solutions for (33) are calculated as follows: and Moreover, we have because as Moreover, a simple computation yields that and we suppose that and , that is, Suppose that (35) holds and Then, it follows that , that is, for every about . Consequently, both solutions of (33) do not lie inside the intersection of the unit circle with the coordinate axes when In the same way, we assume that , Formerly, to change (31) into normal form, we used the following similarity transformation: By using (36), we obtain the next typical form for (31): Moreover, we have where and are respective elements of In addition, for and for are coefficients from expressions and , respectively. Now, we describe the next non-zero real numbers: where and Hence, we have the following theorem.

Assume that ( of system ( In addition, if

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 ) to system (5), we obtain Then, system (38) is controllable provided that its fixed point is locally asymptotically stable. Additionally, the Jacobian matrix for system (38) at its positive fixed point is calculated as follows: where and with Moreover, we have the following result:

Assume the fixed point where

7. Numerical Simulations

First, assume that , and Then, from system (5) we have 8.2317). Moreover, in this case and are initial conditions. The graphical behavior of both variables is shown in Figure 5. It can be seen that and undergo Neimark–Sacker bifurcation at unique positive fixed points 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 , and the Lyapunov exponent changes from negative to positive as h crosses the horizontal line at

Figure 5

Bifurcation diagrams for system (5) for , and with initial conditions and . (a) Bifurcation diagram for ; (b) Bifurcation diagram for ; (c) MLE.

In this case, we have which verifies Theorem 6. Moreover, we have which validates Theorem 5. By varying the stepsize, h, in 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 (see Figure 6c). Second, to discuss the feasibility of the designed control technique, we take and

Figure 6

Phase portraits for system (5) for , and with initial conditions and . (a) Phase portrait for ; (b) phase portrait for ; (c) phase portrait for ; (d) phase portrait for ; (e) phase portrait for ; (f) phase portrait for ; (g) phase portrait for ; (h) phase portrait for .

Then, from system (38), we have 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 is effectively controlled for a large range of the control parameter 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 . Moreover, the controlled system is stable in , and the Lyapunov exponent changes from positive to negative at the point and negative to positive as b crosses the horizontal line at which shows the consistency of controlled plots with corresponding MLE.

Figure 7

Controlled plots of system (38) for , and with initial conditions and . (a) Bifurcation diagram for ; (b) Bifurcation diagram for ; (c) MLE.

Finally, to discuss the dynamics of system (3), we take and Then, from system (3) we have 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 as bifurcation parameter and taking it can be observed from Figure 9 that the system (3) experiences Hopf bifurcation at positive fixed point 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

Figure 8

Plots of system (3) for and with initial conditions , and . (a) Plot of for system (3); (b) plot of for system (3); (c) stable phase portrait for system (3).

Figure 9

Plots of system (3) for and with initial conditions and . (a) Plot of for system (3); (b) plot of for system (3); (c) phase portrait for system (3).

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 , which is considered a bifurcation parameter, and passes over a critical value , such as Furthermore, the first Lyapunov exponent is calculated in the closed form, specified as follows: where 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 whenever h passes through 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, 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.