Stability Analysis and Chaos Diagnosis of a Tumor-Host Immune Cell Interaction Model with Neimark Sacker Bifurcation

Maysoon M. Aziz, Sabaa Abduljabar Mohammed

Mathematics Department, College of Computer Science and Mathematics, University of Mosul, Mosul, Iraq.

DOI: 10.4236/oalib.1109577                

Abstract

In this paper, we take a three-dimension discrete time model of the interaction of host cells with immune cells and tumor cells fixed point analysis was performed to analyze the stability of the system. the necessary conditions have been created to control the growth of cancer cells the introduction of chemotherapy and the disordered behavior was diagnosed the system by finding exponent and dimension of Lyapunov in order, and the numerical simulation of the system was done using the iterative fixed point method, as well as study the dissipative and Neimark-Sacker bifurcation of system. Finally, the tumor cells of the system and its disorder were controlled using the adaptive control technique, and a stable and regular system was obtained.

Keywords

Discrete Tumor-Immune Model, Fixed Points, Stability, Lyapunov Exponent, Dissipative, Neimark-Sacker Bifurcation, Adaptive Control

1. Introduction

Mathematical modeling is one of the tools that can enhance the understanding of physical, chemical and biological phenomena. Physical modeling and chemical modeling are completely different disciplines compared to biological systems modeling. Biological systems are subject to the law of chemistry and physics, yet the specific functions of representing living systems are what distinguish biology from other sciences. In recent years, the use of integer order and fractional order mathematical modeling has become increasingly popular the fixed-point theory has also been used to study the properties of these systems of the study by many researchers [1] [2] [3] [4] [5].

Mathematical models representing the number of cancer cells and their interaction with Immune cells are generally a language that explains complex phenomena; and for more than a century, a lot of resources, whether economic or human, have been expended to fight cancer. The increase in the incidence of cancer among people has proven to be an important area of interdisciplinary research that requires expert medical and biological researchers and recent studies of cancer.

I realized the importance and support of mathematics and computer science in treating cancer and the immune system is the body’s first line of defense against parasites, infections, etc. Although studies of tumor immune dynamics date back to the 1890s, questions about tumor growth and its interaction with immune cells remain today. And diagnosing the tumor early is the most important stage, as the chances of detecting the tumor in the early stages are very less It plays an important role in the development of the tumor immune model and has a vital role in the analysis of parameters that affect the success or failure of the immune response against the tumor [6] [7] [8].

The study of tumor model immune reaction dynamics is complex but has gained great interest and attracted many researchers to research in this field. The dynamic system may undergo changes in the stability zone at various parameters. These changes in the system in the initial state are considered chaotic behavior. Chaotic behavior is not always unexpected and in certain cases, chaos can be predicted the chaos in the tumor can occur it often indicates that the immune reaction model leads to long-term tumor relapse and interpretation of the disorder of tumor growth is the lack of response to treatment. Thus, chaotic behavior somehow provides knowledge of some necessary controls that can be implemented for development in the fight against cancer [9].

Two-dimensional system is dissipative if $\left|\det D_f\right|\le 1$, where $D_f$ is the Jacobi matrix of f at p [10].,A dissipative map satisfies $\left|\det D_f\right|\le b<1$, so that the two Lyapunov numbers of each trajectory satisfy $L_1{L}_2\le b$ [10].

In this paper, the stability of a three-dimensional kinetic system representing the interaction of host cells with cancer cells and immune cells [11] was studied and analyzed, and system chaos was diagnosed. Thus, distributed as follows: system description in Section 2 and the fixed points existence in Section 3, stability analysis (characteristic equation roots test, Jury table, Lyapunov function test) [12] [13] [14] [15] in Section 4, dissipativity [16] in Section 5, numerical analysis and dynamic behavior of system [17], and Neimark-Sacker [18] in Section 6, Lyapunov exponent in Section 7, and adaptive control technique [15] in Section 8.

2. System Description

The nonlinear dynamical system that is three-dimensional and discrete time represent Host-Immune-Tumor Cells Interaction Model given as follows:

$\begin{array}{l}{x}_{t+1}=x_t+r{x}_t\left(1-x_t\right)-{\beta }_1\frac{x_t{z}_t}{1+z_t}\\ y_{t+1}=y_t+\rho y_t{z}_t-{\beta }_2{y}_t{z}_t-\mu y_t\\ z_{t+1}=z_t+r_1{z}_t\left(1-z_t\right)-{\beta }_3\frac{x_t{z}_t}{1+x_t}-{\beta }_4{z}_t{y}_t\end{array}$ (1)

where:

$x_t$ : Host cell population represent.

$y_t$ : Effector Immune cells represent.

$z_t$ : Tumor cells represent.

r: growth rate host cell.

${\beta }_1$ : rate of tumor cells killing host cells.

$\rho$ : rate growth effector immune cells.

${\beta }_2$ : Tumor cells inhibition rate of effector immune cells.

$\mu$ : rate of effector cells natural death.

$r_1$ : tumor cells growth rate.

${\beta }_3$ : rate host cells tumor killing.

${\beta }_4$ : rate tumor killing by effector immune cells.

When $r=0.25$, ${\beta }_1=0.85$, $\rho =0.605$, ${\beta }_2=0.02$, $\mu =0.11$, $r_1=0.2$, ${\beta }_3=0.1$, ${\beta }_4=0.2$.

3. The Fixed Points Existence

The following non-negative points are the fixed points of System (1), obtained by solving the following equations:

$\begin{array}{l}{x}_{t+1}=x_t+r{x}_t\left(1-x_t\right)-{\beta }_1\frac{x_t{z}_t}{1+z_t}=x_t\\ y_{t+1}=y_t+\rho y_t{z}_t-{\beta }_2{y}_t{z}_t-\mu y_t=y_t\\ z_{t+1}=z_t+r_1{z}_t\left(1-z_t\right)-{\beta }_3\frac{x_t{z}_t}{1+x_t}-{\beta }_4{z}_t{y}_t=z_t\end{array}$ (2)

$p_0=\left(0,0,0\right)$

$p_1=\left(1,0,0\right)$

$p_2=\left(0,0,1\right)$

$p_3=\left(0,\frac{r_1\left(\rho -\mu -{\beta }_2\right)}{{\beta }_4\left(\rho -{\beta }_2\right)},\frac{\mu }{\rho -{\beta }_2}\right)$

$p_4=\left(\frac{r\left(\mu +\rho -{\beta }_2\right)-{\beta }_1\mu }{r\left(\mu +\rho -{\beta }_2\right)},\frac{\left(\rho -{\beta }_2\right)\left({\beta }_1\mu {\beta }_3+\left(r_1-{\beta }_3\right)r\left(\mu +\rho -{\beta }_2\right)\right)-{\mu }^2{r}_{1}r\left[\frac{\mu +\rho -{\beta }_2}{\rho -{\beta }_2}\right]}{{\beta }_{4}r{\left(\mu +\rho -{\beta }_2\right)}^2},\frac{\mu }{\rho -{\beta }_2}\right)$

4. Stability Analysis

First, the Jacobi matrix of System (1) must be found before analyzing its stability

$J_{\left(x,y,z\right)}=\left[\begin{array}{ccc}{a}_{11}& 0& a_{13}\\ 0& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right]$ (3)

When:

$a_{11}=1+r\left(1-x\right)-rx-{\beta }_1\frac{z}{1+z}$

$a_{13}=-{\beta }_1\frac{x}{1+z}+{\beta }_1\frac{xz}{{\left(1+z\right)}^2}$

$a_{22}=1+\rho z-{\beta }_{2}z-\mu$

$a_{23}=\rho y-{\beta }_{1}y$

$a_{31}=\frac{{\beta }_{3}z}{1+z}$

$a_{32}=-{\beta }_{4}z$

$a_{33}=1+r_1\left(1-z\right)-r_{1}z-\frac{{\beta }_{3}x}{1+z}+\frac{{\beta }_{3}xz}{{\left(1+z\right)}^2}-{\beta }_{4}y$

a) Characteristic Equation Roots Test

Proposition (1): let

$c\left(\lambda \right)=a_3{\lambda }^3+a_2{\lambda }^2+a_1\lambda +a_0=0$ (4)

Is the characteristic equation of the system, the cases following are true:

1) Sink fixed point if all $\left|{\lambda }_i\right|<1$, $i=1,2,3$ and is asymptotically stable.

2) Saddle fixed point if at least one of $\left|{\lambda }_i\right|>1$, $i=1,2,3$ and is unstable.

3) Source if all $\left|{\lambda }_i\right|>1$ (i.e., the fixed point unstable).

Now, to check the system’s fixed points stability by characteristic equation roots substitute the fixed points into the Jacobi matrix we get characteristic polynomials are the form:

${\lambda }^3-tr\left(J\right){\lambda }^2+{\displaystyle \sum M}\left(J\right)\lambda -\det \left(J\right)=0$ (5)

Were

$tr\left(J\right)={\lambda }_1+{\lambda }_2+{\lambda }_3$

${\displaystyle \sum M}\left(J\right)={\lambda }_1{\lambda }_2+{\lambda }_1{\lambda }_3+{\lambda }_2{\lambda }_3$

$\det \left(J\right)={\lambda }_1{\lambda }_2{\lambda }_3$

To test the stability of the system, we substituted the fixed points in the Jacobi matrix.

Now we check stability of point $p_3$, substitute the point $p_3=\left(0,0.811966,0.188034\right)$ in to the Jacobi matrix weget:

$J_{\left(0,0.811966,0.188034\right)}=\left[\begin{array}{ccc}1.15468& 0& 0\\ 0& 0.99999& 0.4750\\ -0.00158& -0.37607& 0.96239\end{array}\right]$

Find (Det(λI-J) = 0) we get:

${\lambda }^3-3.1171{\lambda }^2+3.24718\lambda -1.13302=0$

From Equation (5) then

$\Rightarrow tr\left(J\right)=3.1171,{\displaystyle \sum M}\left(J\right)=3.24718,\det \left(J\right)=1.13302$

We will get the eigenvalues:

${\lambda }_1=1.1761,{\lambda }_2,{\lambda }_3=0.9705\pm 0.1466i$

Then from Proposition (1) we get

$\left|{\lambda }_1\right|=1.1761>1$ and $\left|{\lambda }_2\right|,\left|{\lambda }_3\right|=0.98150<1$

Then $p_3$ is called saddle point and is unstable.

Now we check stability of point $p_4$, substitute the point $p_4=\left(0.46187,0.61758,0.188034\right)$ in to the Jacobi matrix we get:

$J_{\left(0.46187,0.61758,0.188034\right)}=\left[\begin{array}{ccc}0.88453& 0& 0.27815\\ 0& 0.99999& 0.36128\\ -0.15827& -0.03761& 0.96855\end{array}\right]$

Find (Det(λI-J) = 0) we get:

${\lambda }^3-2.853078{\lambda }^2+2.718973\lambda -0.864326=0$

From Equation (5) then

$\Rightarrow tr\left(J\right)=2.853078,{\displaystyle \sum M}\left(J\right)=2.718973,\det \left(J\right)=0.864326$

We will get the eigenvalues:

${\lambda }_1=0.8630,{\lambda }_2,{\lambda }_3=0.9950\pm 0.1969i$

Then from Proposition (1) we get

$\left|{\lambda }_1\right|=0.6830<1$ and $\left|{\lambda }_2\right|,\left|{\lambda }_3\right|=1.00073>1$

Then $p_4$ is called saddle point and is unstable.

Similarly, we test the point $p_0,p_1,p_2$ and we get all fixed points are unstable, so System (1) is unstable.

Table 1 summarizes the stability results for all fixed points of System (1) using the roots of the characteristic equation.

b) Jury Test

Proposition (2): let

$C\left(\lambda \right)=a_3{\lambda }^3+a_2{\lambda }^2+a_1\lambda +a_0=0$

Characteristic equation from Jacobi matrix for System (1) and Table 2 is the Jury’s table defined as following:

Such that

$b_k=\left|\begin{array}{cc}{a}_0& a_{n-k}\\ a_n& a_k\end{array}\right|$, $n=3$, $k=0,1,2$

$c_k=\left|\begin{array}{cc}{b}_0& b_{n-1-k}\\ b_{n-1}& b_k\end{array}\right|$, $n=3$, $k=0,1,2$

Table 1. Stability of the system using the roots characteristic test.

Table 2. Jury table.

The fixed points are said to be stable if they satisfy the following conditions:

1) $C\left(\text{1}\right)>0$.

2) ${\left(-1\right)}^{n}C\left(-1\right)>0$.

3) $\left|b_0\right|<\left|b_{n-1}\right|$.

4) $\left|c_0\right|<\left|c_{n-2}\right|$.

Otherwise, conditions, the fixed points are unstable.

We check the stability of point $p_0$ from the equation following we got by substituting $p_0=\left(0,0,0\right)$ into the Jacobi matrix

${\lambda }^3-3.34{\lambda }^2+3.6805\lambda -1.335=0$

We form Table 3.

From Proposition (2) Condition (3) is not achieved therefore the point $p_0$ is unstable.

Table 3. Jury table.

and in the same way test all of the fixed points $p_1,p_2,p_3,p_4$ we get all points are unstable, so System (1) is unstable.

c) Lyapunov Function Test

Proposition (3): Let the Lyapunov function equation of System (1) as following: $V_{\left(x_t,y_t,z_t\right)}=x_t^2+y_t^2+z_t^2>0$ When $x_t,y_t,z_t>0$ by using ΔV we get:

$\Delta V=V{\left(x_{t+1}+y_{t+1}+z_{t+1}\right)}^2-V{\left(x_t+y_t+z_t\right)}^2$ (6)

Then called is stable fixed point if they are $\Delta V\le 0$, otherwise the point is unstable.

We check stability of System (1) by Lyapunov function method

Now test the stability of point $p_3$ and $p_4$ substitute in Equation (6) we get:

$\begin{array}{l}\Delta V\left(0,0.811966,0.188034\right)\\ =[{\left(0+0.25\left(0\right)\left(1-0\right)-0.85\frac{0\ast 0.188034}{1+0.188034}\right)}^2+(0.811966+0.605\ast 0.811966\\ \ast 0.188034-0.02\ast 0.811966\ast 0.188034-{0.11\ast 0.811966)}^2\\ +(0.188034+0.2\ast 0.188034\left(1-0.188034\right)-0.1\frac{0\ast 0.188034}{1+0.188034}\\ {\stackrel{}{}-0.2\ast 0.188034*0.811966)}^2]-0^2-{0.811966}^2-{0.188034}^2\\ =1.2389290>0\end{array}$

$\begin{array}{l}\Delta V\left(0.461870,0.617581,0.188034\right)\\ ={[\left(0.461870+0.25\ast 0.461870\left(1-0.461870\right)-0.85\frac{0.461870\ast 0.188034}{1+0.188034}\right)}^2\\ +(0.617581+0.605\ast 0.617581\ast 0.188034-0.02\ast 0.617581\ast 0.188034\\ -{0.11\ast 0.617581)}^2+(0.188034+0.2\ast 0.188034\left(1+0.188034\right)\stackrel{}{}\\ -{0.1\frac{0.461870\ast 0.188034}{1+0.188034}-0.2\ast 0.188034\ast 0.617581)}^2-{0.461870}^2\\ -{0.617581}^2-{0.188034}^2=5.350775083\ast {10}^{-3}>0\end{array}$

Then $p_3$ and $p_4$ are unstable points.

In the same way, we tested the rest of the fixed points and got an unstable system.

5. Dissipativity

The Jacobi matrix of the point $p_1$ is:

$J=\left[\begin{array}{ccc}0.75& 0& -0.85\\ 0& 0.89& 0\\ 0& 0& 1.1\end{array}\right]$

when calculating its determinant, we get:

$\left|\det J_{\left(1,0,0\right)}\right|=0.73425<1$

Then is dissipative, and similarly we calculate the rest of the fixed points and we get a dissipative system.

6. Numerical Analysis and Dynamic Behavior

a) Fixed point iteration method

In this section Fixed point iteration method was used to find the solution of System (1) with less error (0.9) written program in MABLE gives a result with (x, y, z) = (0.29315, 0.68038, 0.12214).

b) Time Behavior of System (1)

A time series of System (1) was generated for 1000 iterations with System (1) with the values we got from the fixed point iteration method we get Figure 1: [(1): x(t) and time(t), (2): y(t) and time(t), (3): z(t) and time(t)].

We get phase space of the system in Figure 2.

c) Neimark-Sacker Bifurcation

For a Neimark-Sacker bifurcation to occur the eigenvalues ( ${\lambda }_{1,2,3}$ ) must satisfy the following conditions:

1) At the bifurcation point the true eigenvalue is not equal to ±1.

2) The absolute value of the eigenvalue equals 1.

The bifurcation from the eigenvalues at point $p_4$ was obtained from the bifurcation taking into account the above from the parameters are ${\lambda }_1=0.8630$ and ${\lambda }_{2,3}=0.9950\pm 0.1969i$ ( $\left|{\lambda }_{1,2}\right|=1$ ).

Thus, the system undergoes Neimark-Sacker bifurcation at r = 0.25 in the fixed point $p_4$ where all the conditions for a Neimark-Sacker bifurcation are satisfied. Figure 3 shows us the bifurcation diagram.

7. Lyapunov Exponent

In this part of the article, we will use the Lyapunov exponential test, to study the chaos of System (1) and if at least one value of the Lyapunov exponent is greater than zero, then this means that the system is chaotic, and the values we obtained are as follows: $L_1=1.250000$, $L_2=0.890000$, $L_3=-2.55334$.

We use the following basic formula to calculate the Lyapunov dimension $D_L=2+\frac{L_1+L_2}{\left|L_3\right|}=2+\frac{1.250000+0.890000}{\left|-2.55334\right|}=2.8381179$

Figure 1. Time series of System (1).

Figure 2. Time trajectories of the system variables $x_t,y_t,z_t$.

Figure 3. Bifurcation diagram.

Thus, System (1) is highly chaotic and Figure 4 shows it.

8. Adaptive Control Technique

We will use adaptive control technology to manipulate the chaos of the chaotic System (1) and design a law for the adaptive control with the unknown parameter r, we get the following system after adding the controllers to System (1):

Figure 4. Lyapunov exponent.

$\begin{array}{l}{x}_{t+1}=x_t+r{x}_t\left(1-x_t\right)-{\beta }_1\frac{x_t{z}_t}{1+z_t}+u_1\\ y_{t+1}=y_t+\rho y_t{z}_t-{\beta }_2{y}_t{z}_t-\mu y_t+u_2\\ z_{t+1}=z_t+r_1{z}_t\left(1-z_t\right)-{\beta }_3\frac{x_t{z}_t}{1+x_t}-{\beta }_4{z}_t{y}_t+u_3\end{array}$ (7)

As they $u_1,u_2$ are the control units for air-conditioned feeding and defined as follows:

$u_1=-x_t-\hat{r}{x}_t\left(1-x_t\right)+{\beta }_1\frac{x_t{z}_t}{1+z_t}-M{x}_t$

$u_1=-y_t-\rho y_t{z}_t+{\beta }_2{y}_t{z}_t+\mu y_t-M{y}_t$ (8)

$u_1=-z_t-r_1{z}_t\left(1-z_t\right)+{\beta }_3\frac{x_t{z}_t}{1+x_t}+{\beta }_4{z}_t{y}_t-M{z}_t$

Since $M_1,M_2$ they are positive numbers and the parameter $\hat{r}$ is an estimation parameter for parameter r and with a substitution (8) in (7)

We get:

$\begin{array}{l}{x}_{t+1}=e_r{x}_t+\left(1-x_t\right)-M_1{x}_t\\ y_{t+1}=-M_2{y}_t\\ z_{t+1}=-M_3{z}_t\end{array}$ (10)

Since $e_r=\left(r-\hat{r}\right)$ is the error for the estimating parameter and $M_1=0.3$, $M_2=0.6$, $M_3=0.2$ and $\hat{r}=0.24$ is the [10]. estimated parameter of r.

Fixed point stability analysis of the System (10).

a) Characteristic equation roots

Before analyzing the stability of the System (10), we find the Jacobi matrix of the system:

$J\left(x_t,y_t,z_t\right)=\left[\begin{array}{ccc}{e}_r\left(1-x_t\right)-e_r{x}_t-M_1& 0& 0\\ 0& -M_2& 0\\ 0& 0& -M_3\end{array}\right]$ (11)

To test the point $p_0$, we substitute the fixed point $p_0=\left(0,0,0\right)$ in Equation (11) we get:

$J_{\left(0,0,0\right)}=\left[\begin{array}{ccc}-0.3& 0& 0\\ 0& -0.6& 0\\ 0& 0& -0.2\end{array}\right]$

And by finding the determinant (Det(λI-J) = 0), we get:

${\lambda }^3+1.1{\lambda }^2+0.36\lambda +0.036=0$ (12)

$tr\left(J\right)=-1.1$, ${\displaystyle \sum M}\left(J\right)=0.36$, $\det \left(J\right)=-0.036$

We will get the eigenvalues:

${\lambda }_1=-0.6$, ${\lambda }_2=-0.3$, ${\lambda }_3=-0.2$

According to Proposition (1), we obtain:

$\left|{\lambda }_1\right|=0.6$, $\left|{\lambda }_2\right|=0.3$, $\left|{\lambda }_3\right|=0.2$, Therefore, the point $p_0$ is stable.

The same method was tested for the rest of the points, and the results of the test points obtained are shown in Table 4 which turned out to be stable points, that is System (10) is stable.

b) Jury test

To test the point $p_0$, from the coefficients of the characteristic Equation (12) for the point, we get:

$a_0=0.036,a_1=0.36,a_2=1.1,a_3=1$

Accordingly, we form Table 5, which represents a Jury’s table for point $p_0$, as follows:

From Table 5, it is clear that all the conditions of the Jury test for point $p_0$ are fulfilled (Proposition (2)) $p_0$ is stable.

In the same way, the rest of the fixed points were tested, which turned out to be stable as well the System (10) is a stable system.

c) Lyapunov function test

Table 4. Results of the test fixed points using the test for the roots of the characteristic equation.

Here we will analyze the stability of the System (10) using the Lyapunov function test, now test for the fixed point $p_1$ based on Proposition (3) we get:

$\Delta V{p}_1\left(0,0,0\right)={\left(0.01\left(1-1\right)-0.3\left(1\right)\right)}^2+{\left(-0.6\left(0\right)\right)}^2+{\left(-0.2\left(0\right)\right)}^2-1^2=-0.91<0$

Then the fixed point $p_1$ is stable. In the same way, the stability of the rest of the fixed points was tested, and we obtained that all the fixed points are stable, and this indicates that the system is stable.

d) Lyapunov Exponent of system

We diagnose the chaos in the controlled System (10) using the adaptive control technique of the Lyapunov exponent of System (10) and the following values were obtained:

$L_1=-0.2$, $L_2=-0.29$, $L_3=-0.6$

Since all the obtained values are negative, this indicates that System (10) is regular, and Figure 5 shows us the behavior of the controlled system.

Table 5. Jury table of $p_0$.

Figure 5. Lyapunov exponent after adaptive control.

9. Conclusion

The separate three-dimensional immune tumor cell model that was dealt with in this research, and after finding the fixed points of the system, five fixed points were obtained, all of which are real points, so that the system is self-attracting, and the stability of the system was analyzed using (the characteristic equation test and Jury’s table, and the Lyapunov function test), and we got an unstable system for all fixed points. The chaos of the system was diagnosed and the Lyapunov exponent was found where two values were found, which two positives are 1.250000 and 0.890000, and a third negative is −2.55334, and thus the system was very chaotic. As well as finding the Neimark-sacker bifurcation of the system at the complex eigenvalues resulting from the substitution of the fixed point $p_4$ in the Jacobi matrix, the application of adaptive control technology on the system, thus controlling the system and obtaining a new, stable and regular system with Lyapunov exponent all are negative.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Ravichandran, C., Logeswari, K., Panda, S.K. and Nisar, K.S. (2020) On New Approach of Fractional Derivative by Mittag-Leffler Kernel to Neutralintegro-Differential Systems with Impulsive Conditions. Chaos, Solitons Fractals, 139, Article ID: 110012. https://doi.org/10.1016/j.chaos.2020.110012
[2]	Panda, S.K., Karapinar, E. and Atangana, A. (2020) A Numerical Schemes and Comparisons for Fixed Point Results with Applications to the Solutions of Volterra Integral Equations in Dislocated Extended B-Metric Space. Alexandria Engineering Journal, 59, 815-827.
[3]	Panda, S.K., Abdeljawad, T. and Ravichandran, C. (2020) Novel Fixed Point Approach to Atangana-Baleanu Fractional and Lp-Fredholm Integral Equations. Alexandria Engineering Journal, 59, 1959-1970. https://doi.org/10.1016/j.aej.2019.12.027
[4]	Tuan, N.H., Mohammadi, H. and Rezapour, S. (2020) A Mathematical Model for COVID-19 Transmission by Using the Caputo Fractional Derivative. Chaos, Solitons Fractals, 140, Article ID: 110107. https://doi.org/10.1016/j.chaos.2020.110107
[5]	Rezapour, S., Mohammadi, H. and Jajarmi, A. (2020) A New Mathematical Model for Zika Virus Transmission. Advances in Difference Equations, 2020, Article No. 589. https://doi.org/10.1186/s13662-020-03044-7
[6]	Wodarz, D. (2005) Mathematical Models of Immune Effector Responses to Viral Infections Virus Control versus the Development of Pathology. Journal of Computational and Applied Mathematics, 184, 301-319. https://doi.org/10.1016/j.cam.2004.08.016
[7]	Wodarz, D., May, R.M. and Nowak, M.A. (2007) The Role of Antigen Independent Persistence of Memory Cytotoxic T Lymphocytes. International Immunology, 12, 467-477. https://doi.org/10.1093/intimm/12.4.467
[8]	Ferrari, C. (2015) Hbv and the Immune Response. Liver International, 35, 121-128. https://doi.org/10.1111/liv.12749
[9]	Khajanchi, S. and Nieto, J.J. (2019) Mathematical Modeling of Tumor-Immune Competitive System, Considering the Role of Time Delay. Applied Mathematics and Computation, 340, 180-205. https://doi.org/10.1016/j.amc.2018.08.018
[10]	Alligood, K.T., Sauer, T.D. and Yorke, J.A. (1997) Chaos: An Introduction to Dynamical System. Springer-Verlag, New York. https://doi.org/10.1007/978-3-642-59281-2
[11]	Alzabut, J., Selvam, A.G.M., Dhakshinamoorthy, V., Mohammadi, H. and Rezapour, S. (2022) On Chaos of Discrete Time Fractional Order Host-Immune-Tumor Cells Interaction Model. Journal of Applied Mathematics and Computing, 68, 4795-4820. https://doi.org/10.1007/s12190-022-01715-0
[12]	Aziz, M.M. and Jihad, O.M. (2021) Stability, Chaos Tests with Adaptive and Feedback Control Methodsfor 3D Discrete-Time Dynamical System. International Journal of Electronics Communication and Computer Engineering, 12, 31-42.
[13]	Aziz, M.M. and Jihad, O.M. (2021) Stability & Chaos Tests of 2D Discrete Time Dynamical System with Hidden Attractors. Open Access Library Journal, 8, e7501. https://doi.org/10.4236/oalib.1107501
[14]	Aziz, M.M. (2018) Stability Analysis of Mathematical Model. International Journal of Science and Research (IJSR), 7, 147-148.
[15]	Aziz, M.M. and Merie, D.M. (2020) Stability and Adaptive Control with Synchronization of 3-D Dynamical System. Open Access Library Journal, 7, Article No. e6075. https://doi.org/10.4236/oalib.1106075
[16]	Aziz, M.M. and Merie, D.M. (2020) Stability and Chaos with Mathematical Control of 4-D Dynamical System. Indonesian Journal of Electrical Engineering and Computer Science, 20, Article 1242.
[17]	Liu, X. and Xiao, D. (2007) Complex Dynamic Behaviors of a Discrete-Time Predator-Prey System. Chaos, Solitons Fractals, 32, 80-94. https://doi.org/10.1016/j.chaos.2005.10.081
[18]	Xin, B., Chen, T. and Ma, J.H. (2010) Neimark-Sacker Bifurcation in a Discrete-Time Financial System. Discrete Dynamics in Nature and Society, 2010, Article ID: 405639. https://doi.org/10.1155/2010/405639