A Discrete Dynamics Approach to a Tumor System

Abstract

In this paper, we present a cancer system in a continuous state as well as some numerical results. We present discretization methods, e.g., the Euler method, the Taylor series expansion method, and the Runge–Kutta method, and apply them to the cancer system. We studied the stability of the fixed points in the discrete cancer system using the new version of Marotto’s theorem at a fixed point; we prove that the discrete cancer system is chaotic. Finally, we present numerical simulations, e.g., Lyapunov exponents and bifurcations diagrams.

1. Introduction

Cancer is a group of diseases involving abnormal cell growth. These cells could form a mass known as a tumor. A malignant tumor means it could invade or spread into nearby tissues. Therefore, cancer cells (linked with tumor growth) could reach distant parts of the body and form a new tumor (far from the original one). A benign tumor means the tumor can grow but not spread.

Many authors have used mathematical models. The main components of these models involve interactions among three types of cells, cancer cells with healthy host cells, and cells of the immune system. These interactions may lead to different outcomes. In the literature, there are several reviews on mathematical systems applied to tumor dynamics, e.g., [1,2,3,4].

An approach using a discrete-time system is better to describe the tumor dynamics compared to a continuous-time system (when populations have non-overlapping generations). Regarding computational and numerical simulations, discrete-time models are more efficient. For the advantages of the discrete approach compared to the continuous approach, see [4,5,6,7].

Chaos can be found in many biological systems. Chaotic systems have internal behaviors that depend (in many cases) on the initial conditions and could be suppressed as a result of the application of small perturbations.

Chaotic behaviors are complex in general, full of irregularities, and are generally undesirable in biological systems. Model complexities, in many cases, need to be reduced in mathematical modeling to generate dynamic results. Discretization methods, e.g., the Taylor series expansion or Euler and Runge–Kutta methods, are prime tools used to treat chaotic systems. These methods are of particular importance when studying differential equations applied to tumor dynamics. Selecting the best discretization approach is a problem in itself. We endeavored to compare the Euler and Runge–Kutta numerical integration approach to simulate the chaotic behavior of a multi-scroll chaotic oscillator and compare the obtained results.

In this paper, we consider the model presented by Pillis and Radunskaya; see [1]. In the second section (after the introduction), we present the cancer system in a continuous state (with some numerical results). In the third section, we present discretization methods, for example, the Euler method, the Taylor series expansion method, and the Runge–Kutta method; we applied them to the cancer system. In the fourth section, we review the stability of the fixed points in the discrete cancer system. In the fifth section, we prove that the discrete cancer system is chaotic, using the new version of Marotto’s theorem at a fixed point. In the sixth section, we present numerical simulations (e.g., Lyapunov exponents and bifurcations diagrams). Finally, we present a conclusion.

2. The Continuous Version of the Cancer System

The AIMS model is used to describe the competition and interactions among tumor cells, healthy host cells, and effector immune cells. Regarding cancer models that include interacting cells, we focused on cells near the tumor sites. Populations are based on the prey–predator models and the law of exponential growth. Although the previous models are simple, they are suitable platforms to explain important aspects of the dynamics of cancer. In this section, we will present the model to describe the biological tumor system, which is given in the form of an ordinary differential equation as follows

$$\{\begin{array}{c}\dot{N}={\rho }_{2}N(1-b_{2}N)-c_{4}NT,\\ \dot{T}={\rho }_{1}T(1-b_{1}T)-c_{2}TI-c_{3}TN,\\ \dot{I}=\left(\frac{\rho IT}{\alpha +T}\right)-c_{1}IT-d_{1}I+s,\end{array}$$

(1)

where N denotes the healthy host cells, T denotes the number of cancer cells, I denotes the effector immune cells, and ${\rho }_1,{\rho }_2,\rho ,b_1,b_2,\alpha ,c_1,c_2,c_3,c_4,d_1$, and s are positive parameters; see [1,2,4]. Here, ${\rho }_1$ represents the growth rate (in the absence of any effect) of cancer cells from other cell populations with a maximum carrying capacity of $1/b_1$; the values $c_2$ and $c_3$ refer to the ’killing rate’ of the cancer cells by the healthy host cells and effector cells, respectively; ${\rho }_2$ represents the growth rate, with a maximum carrying capacity $1/b_2$ of healthy host cells; $c_4$ represents the rate of inactivation of the healthy cells by the cancer cells. The rate (or level) of recognition of the cancer cells by the immune system depends usually on the antigenicity of the cancer cells. Due to the large complexity of this recognition process and to keep the model simpler, we assume that the stimulation of the immune system depends—in a direct way—on the number of cancer cells with positive constants, $\rho$ and $\alpha$. We consider that the effector cells are inactivated by the cancer cells at rate $c_1$ and that they die in a natural way at rate $d_1$. The value s is a constant influx of immune cells.

To simplify the study of this system (1). we reduced the number of parameters by introducing this change of variables: $x=b_{1}T,y=b_{2}N,z=\frac{I}{\alpha }\mathrm{et}\tau ={\rho }_{1}t,$ and the new parameters: $a_{12}=\frac{c_2}{b_2{\rho }_1},r_2=\frac{{\rho }_1}{{\rho }_2},a_{21}=\frac{c_4}{b_1{\rho }_1},r_3=\frac{\rho }{{\rho }_1},$ $k_3=\alpha b_1,a_{31}=\frac{c_1}{b_1{\rho }_1},$ $a_{13}=\frac{\alpha c_3}{{\rho }_1}\mathrm{and}{d}_3=\frac{d_1}{{\rho }_1}.$ then system (1) is converted to, see Figure 1, Figure 2, Figure 3 and Figure 4:

$$\{\begin{array}{c}\dot{x}=x(1-x)-a_{12}xy-a_{13}xz,\\ \dot{y}=r_{2}y(1-y)-a_{21}xy,\\ {\displaystyle \dot{z}=r_3\left(\frac{xz}{x+k_3}\right)-a_{31}xz-d_{3}z.}\end{array}$$

(2)

The Lyapunov exponents of system (2) are computed to be ${\lambda }_1=0.021909,$ ${\lambda }_2=-0.00085097$ and ${\lambda }_3=-0.54025$. The Lyapunov dimension for system (2) is $D_L=2+\frac{{\lambda }_1+{\lambda }_2}{{\lambda }_3}\simeq 2.04.$

3. Discretization Methods

3.1. Euler Discretization Method

The simplest method used for approximating the solution of (2) is called the Euler method, named after Leonhard Euler; see [7,8]. The expression of the Euler method is given in Equation (3) and the discretized model is expressed in Equation (4).

$$\{\begin{array}{l}{\displaystyle \dot{x}\left(t\right)=\frac{x(t+h)-x\left(t\right)}{h}}\\ {\displaystyle \dot{y}\left(t\right)=\frac{y(t+h)-y\left(t\right)}{h}}\\ {\displaystyle \dot{z}\left(t\right)=\frac{z(t+h)-z\left(t\right)}{h},}\end{array}$$

(3)

$$\{\begin{array}{l}{x}_{k+1}=(x_k(1-x_k)-a_{12}{x}_k{y}_k-a_{13}{x}_k{z}_k).h+x_k,\\ y_{k+1}=(r_2{y}_k(1-y_k)-a_{21}{x}_k{y}_k).h+y_k,\\ {\displaystyle \hspace{1em}{z}_{k+1}=(r_3(\frac{x_k{z}_k}{x_k+k_3})-a_{31}{x}_k{z}_k-d_3{z}_k).h+z_k.}\end{array}$$

(4)

3.2. Taylor Series Expansion Method

In this section, we give a numerical method to compute the numerical solutions of (2) by using a Taylor polynomial; see [6,8] for $x(t+h),y(t+h)$ and $z(t+h)$,

$$\{\begin{array}{l}{\displaystyle x(t+h)=x\left(t\right)+\sum _{m+1}^{\infty }\frac{1}{m!}.h^m.x^{\left(m\right)}}\\ {\displaystyle y(t+h)=y\left(t\right)+\sum _{m+1}^{\infty }\frac{1}{m!}.h^m.y^{\left(m\right)}}\\ {\displaystyle z(t+h)=z\left(t\right)+\sum _{m+1}^{\infty }\frac{1}{m!}.h^m.z^{\left(m\right)}}\end{array}$$

(5)

The Taylor series expansion method is performed for $m=2$ and h. In this setting, we obtain the equations of the discrete time state for the cancer system as follows.

$$\{\begin{array}{l}{\displaystyle x_{k+1}=x_k+h.{\dot{x}}_k+\frac{1}{2}.h^2.{\ddot{x}}_k}\\ {\displaystyle y_{k+1}=y_k+h.{\dot{y}}_k+\frac{1}{2}.h^2.{\ddot{y}}_k}\\ {\displaystyle z_{k+1}=z_k+h.{\dot{z}}_k+\frac{1}{2}.h^2.{\ddot{z}}_k,}\end{array}$$

(6)

where

$$\{\begin{array}{l}{\dot{x}}_k=x_k(1-x_k)-a_{12}{x}_k{y}_k-a_{13}{x}_k{z}_k,\\ {\dot{y}}_k=r_2{y}_k(1-y_k)-a_{21}{x}_k{y}_k,\\ {\displaystyle \hspace{1em}{\dot{z}}_k=r_3\left(\frac{x_k{z}_k}{x_k+k_3}\right)-a_{31}{x}_k{z}_k-d_3{z}_k.}\end{array}$$

(7)

3.3. Runge–Kutta Discretization Method

The Runge-Kutta fourth order method was executed for h. In this situation, the discrete cancer system (4) with the Runge-Kutta fourth order method (see [9,10]), see Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10:

$$\begin{array}{c}{\alpha }_1=h.f(x_k,y_k,z_k)=h.x_{k+1},\\ l_1=h.g(x_k,y_k,z_k)=h.y_{k+1},\\ m_1=h.p(x_k,y_k,z_k)=h.z_{k+1},\\ {\alpha }_2=h.f((x_k+\frac{1}{2}{\alpha }_1,(y_k+\frac{1}{2}{l}_1),(z_k+\frac{1}{2}{m}_1)),\\ l_2=h.g((x_k+\frac{1}{2}{\alpha }_1,(y_k+\frac{1}{2}{l}_1),(z_k+\frac{1}{2}{m}_1)),\\ m_2=h.p((x_k+\frac{1}{2}{\alpha }_1,(y_k+\frac{1}{2}{l}_1),(z_k+\frac{1}{2}{m}_1)),\\ {\alpha }_3=h.f((x_k+\frac{1}{2}{\alpha }_2,(y_k+\frac{1}{2}{l}_2),(z_k+\frac{1}{2}{m}_2)),\\ l_3=h.g((x_k+\frac{1}{2}{\alpha }_2,(y_k+\frac{1}{2}{l}_2),(z_k+\frac{1}{2}{m}_2)),\\ m_3=h.p((x_k+\frac{1}{2}{\alpha }_2,(y_k+\frac{1}{2}{l}_2),(z_k+\frac{1}{2}{m}_2)),\\ {\alpha }_4=h.f((x_k+{\alpha }_3),(y_k+l_3),(z_k+m_3)),\\ l_4=h.g((x_k+{\alpha }_3),(y_k+l_3),(z_k+m_3)),\\ m_4=h.p((x_k+{\alpha }_3),(y_k+l_3),(z_k+m_3)),\\ x_{k+1}=x_k+\frac{1}{6}({\alpha }_1+2{\alpha }_2+2{\alpha }_3+{\alpha }_4),\\ y_{k+1}=y_k+\frac{1}{6}(l_1+2l_2+2l_3+l_4),\\ z_{k+1}=z_k+\frac{1}{6}(m_1+2m_2+2m_3+m_4).\end{array}$$

(8)

4. Stability Analysis for Discrete Cancer System

To find the fixed points, the three discrete cancer equations were set to $x,y,z$ coordinates of each fixed point, determined by solving the following equations:

$$\{\begin{array}{l}x=(x(1-x)-a_{12}xy-a_{13}xz).h+x,\\ y=(r_{2}y(1-y)-a_{21}xy).h+y,\\ {\displaystyle \hspace{1em}z=(r_3(\frac{xz}{x+k_3})-a_{31}xz-d_{3}z).h+z.}\end{array}$$

To obtain the fixed points of the system (4), we set

$$\left\{\begin{array}{l}x=0,\\ x=1-a_{12}y-a_{13}z.\end{array}\right.$$

(9)

$$\left\{\begin{array}{cc}\hfill y& =0,\hfill \\ & {\displaystyle y=\frac{1}{r_2}-\frac{a_{21}}{r_2}x.}\hfill \end{array}\right.$$

(10)

$$\left\{\begin{array}{l}z=0,\\ {\displaystyle x^2+(k_3+\frac{d_3-r_3}{a_{31}})x+\frac{k_3{d}_3}{a_{31}}=0.}\end{array}\right.$$

(11)

The solutions from Equations (9)–(11) together yielded to six fixed points. We discussed their local behaviors according to their biological relevance. Now, we will look at the stabilities of these fixed points.

In this paper, we studied h in interval $[0.01,0.1]$

(1) The first fixed point is trivial and given as $v_1=(0,0,0)$, the corresponding eigenvalues are ${\lambda }_1=h+1$, ${\lambda }_2=h{r}_2+1$ and ${\lambda }_3=-h{d}_3+1$. Since h is small positive, all the parameters are positive, and $\mid {\lambda }_i\mid <1$ $(i=1,2)$; therefore,

Proposition 1.

• If $h{d}_3>2$ then $\mid {\lambda }_3\mid >1$, we have a saddle at this fixed point.
• If $h{d}_3<2$ then $\mid {\lambda }_3\mid <1$, we have a node stable at this fixed point.

(2) The second fixed point is $v_2=(0,1,0)$; the Jacobian matrix evaluated at $v_2$ is given by

$$J\left(v_2\right)=\left[\begin{array}{ccc}(-a_{12}+1)h+1& 0& 0\\ -a_{21}h& -r_{2}h+1& 0\\ 0& 0& -d_{3}h+1\end{array}\right]$$

(12)

Clearly, $J\left(v_2\right)$ has eigenvalues ${\lambda }_1=1-r_{2}h$, ${\lambda }_2=1+(-a_{12}+1)h$ and ${\lambda }_3=1-d_{3}h$. where $h\in [0.01,0.1]$. In fact, in biology, $r_2,d_3$ are smaller than $h^{-1}$. Then $\mid {\lambda }_1\mid <1$ and $\mid {\lambda }_3\mid <1$. The stability of this fixed point depends on the value of parameter $a_{12}$, if $a_{12}<1$ then ${\lambda }_2>1$, this fixed point has two stable and one unstable eigenvalue. Therefore, we have a saddle at $v_2$, and if $a_{12}>1$, then ${\lambda }_2<1$; this fixed point has three stable eigenvalues. Therefore, we have a node at this fixed point. If $a_{12}=1$, then ${\lambda }_2=1$; therefore, we cannot give any information on the stability of $v_2$. In the numerical simulations, we obtained different results depending on the values of $a_{12}$. We observed that the chaotic dynamics appeared close to $a_{12}=1$. The selection of $a_{12}<1$ provides different dynamical behaviors, such as convergence to a stable spiral. However, in this study, we focus on parameter $a_{12}$, where we have chaotic attraction.

(3) The third fixed point is $v_3=(1,0,0)$; the Jacobian matrix evaluated at $v_2$ is given by

$$J\left(v_3\right)=\left[\begin{array}{ccc}-h+1& -a_{12}h& -a_{13}h\\ 0& (r_2-a_{21})h+1& 0\\ 0& 0& ({\displaystyle \frac{r_3}{1+k_3}}-a_{31}-d_3)h+1\end{array}\right]$$

(13)

The eigenvalues of the Jacobian matrix (13) at this fixed point are obtained as ${\lambda }_1=-h+1$, ${\lambda }_2=(r_2-a_{21})h+1$ and ${\lambda }_3=({\displaystyle \frac{r_3}{1+k_3}}-a_{31}-d_3)h+1$. Then $\mid {\lambda }_1\mid <1$. Moreover, ${\lambda }_1$ is stable, and ${\lambda }_2$, ${\lambda }_3$ are stable with the selected parameters.

(4) The fourth fixed point is $v_4=(x^{*},0,z^{*})$. The Jacobian matrix evaluated at $v_4$ is given by

$$J\left(v_4\right)=\left(\begin{array}{ccc}{L}_{11}& L_{12}& L_{13}\\ 0& L_{22}& 0\\ L_{31}& 0& L_{33},\end{array}\right)$$

where

• $L_{11}=(1-a_{13}\overline{z}-2\overline{x})h+1,$
• $L_{12}=-a_{12}\overline{x}h,$
• $L_{13}=-a_{13}\overline{x}h,$
• $L_{22}=(r_2-a_{21}\overline{x})h+1,$
• $L_{31}=(\frac{r_3\overline{z}}{\overline{x}+k_3}-\frac{r_3\overline{\mathit{xz}}}{{(\overline{x}+k_3)}^2}-a_{31}\overline{z})h,$
• $L_{33}=(\frac{r_3\overline{x}}{\overline{x}+k_3}-a_{31}\overline{x}-d_3)h+1.$

The eigenvalues of the Jacobian matrix at this point are

$${\lambda }_1=L_{22}=(r_2-a_{21}\overline{x})h+1,$$

(14)

$${\lambda }_{2,3}=\frac{1}{2}\left[(L_{11}+L_{33})\mp \sqrt{{(L_{11}-L_{33})}^2+4L_{31}{L}_{13}}\right].$$

(15)

(i) 

If ${(L_{11}-L_{33})}^2+4L_{31}{L}_{13}>0$, we have three real eigenvalues.

(ii) 

If ${(L_{11}-L_{33})}^2+4L_{31}{L}_{13}<0$, we have one real and two complex eigenvalues that are stable with the selected parameter sets.

The characteristic equation of the Jacobian matrix $J\left(v_4\right)$ can be expressed in the form

$$P\left(\lambda \right)={\lambda }^3+A_2{\lambda }^2+A_1\lambda +A_0,$$

(16)

where

• $A_0=-L_{33}{L}_{22}{L}_{11}+L_{31}{L}_{13}{L}_{22},$
• $A_1=L_{11}{L}_{22}+L_{11}{L}_{33}-L_{13}{L}_{31}+L_{33}{L}_{22},$
• $A_2=-L_{33}-L_{22}-L_{11}.$

According to the Jury conditions [11], to find the asymptotically-stable region of $v_4$, it is necessary to find the region holding these conditions:

$$P\left(1\right)>0,P(-1)<0,\mid A_0\mid <A_n,\mid B_0\mid >\mid B_{n-1}\mid$$

where $B_k=\left|\begin{array}{cc}{A}_0& A_{n-k}\\ A_n& A_k\end{array}\right|.$ Then $P\left(1\right)=1+A_2+A_1+A_0,$ $P(-1)=-1+A_2-A_1+A_0,$

According to the relations $P\left(1\right)>0,P(-1)<0,\mid A_0\mid <A_n,\mid B_0\mid >\mid B_{n-1}\mid ,$ we have

$\mid A_0\mid <1$, $\mid A_0+1\mid >\mid A_1\mid$ and $\mid A_0-1\mid \mid A_0+A_1+1\mid >\mid A_0{A}_1-A_2\mid$.

(5) The fifth fixed point is ${\displaystyle v_5=(\frac{r_2(a_{12}-1)}{a_{12}{a}_{21}-r_2},{\displaystyle \frac{a_{12}-r_2}{a_{12}{a}_{21}-r_2}},0)}$, where $a_{12}{a}_{21}-r_2\ne 0$. The Jacobian matrix of system (4) at $v_5$ is given by

$$J\left(v_5\right)=\frac{1}{q}\left(\begin{array}{ccc}{M}_{11}& M_{12}& M_{13}\\ M_{21}& M_{22}& 0\\ 0& 0& M_{33}\end{array}\right),$$

where $q=a_{12}{a}_{21}-r_2$ and

• $M_{11}=-h{a_{12}}^2-r_{2}h{a}_{12}+2r_{2}h+2q$
• $M_{12}=-a_{12}{r}_2(a_{12}-1)h$
• $M_{13}=-a_{13}{r}_2(a_{12}-1)h$
• $M_{21}=-a_{21}(a_{12}-r_2)h$
• $M_{22}=-h{a}_{12}{a}_{21}{r}_2+hq{r}_2-2r_{2}h{a}_{12}+r_{2}h{a}_{21}+q$
• $M_{33}=(\frac{r_3{r}_{2}q(a_{12}-1)}{r_2(a_{12}-1)+q{k}_3}-a_{31}{r}_2(a_{12}-1)-d_{3}q)h+q$

The characteristic equation of the Jacobian matrix $J\left(v_6\right)$ can be written as

$$P^{*}\left(\lambda \right)={\lambda }^3+B_2{\lambda }^2+B_1\lambda +B_0=0.$$

where

• $B_0=-{\displaystyle \frac{M_{33}(M_{22}{M}_{11}-M_{21}{M}_{12})}{q^3}},$
• $B_1={\displaystyle \frac{M_{22}{M}_{11}+M_{33}{M}_{11}-M_{21}{M}_{12}+M_{33}{M}_{22}}{q^2}},$
• $B_2=-{\displaystyle \frac{M_{33}+M_{22}+M_{11}}{q}}.$

The eigenvalues of the Jacobian matrix at this fixed point are ${\displaystyle {\lambda }_1=\frac{M_{33}}{q},}$ and ${\displaystyle {\lambda }_{2,3}=\frac{1}{2q}(M_{22}+M_{11}\mp \sqrt{\Delta }),}$ where $\mathsf{\Delta }={M_{11}}^2-2M_{22}{M}_{11}+4M_{21}{M}_{12}+{M_{22}}^2$,

(6) The sixth fixed point is a nontrivial $v_6=(x^{*},y^{*},z^{*})$. The Jacobian matrix of system (4) at $v_6$ is given by

$$J\left(v_6\right)=\left(\begin{array}{ccc}{S}_{11}& S_{12}& S_{13}\\ S_{21}& S_{22}& 0\\ S_{31}& 0& S_{33}\end{array}\right),$$

(17)

where

• $S_{11}=(-a_{12}{y}^{*}-a_{13}{z}^{*}-2x^{*}+1)h+1,$
• $S_{12}=-a_{12}{x}^{*}h,$
• $S_{13}=-a_{13}{x}^{*}h,S_{21}=-a_{21}{y}^{*}h,$
• $S_{22}=-2y^{*}h{r}_2-x^{*}h{a}_{21}+h{r}_2+1$
• $S_{31}=z^{*}(a_{31}{\left(x^{*}\right)}^2+2a_{31}{x}^{*}{k}_3+a_{31}{k}_3^2-r_3{k}_3)h{(x^{*}+k_3)}^2,$
• $S_{33}=(\frac{r_3{x}^{*}}{x^{*}+k_3}-a_{31}{x}^{*}-d_3)h+1.$

The characteristic equation of the Jacobian matrix $J\left(v_6\right)$ can be written as

$$P^{*}\left(\lambda \right)={\lambda }^3+C_2{\lambda }^2+C_1\lambda +C_0=0.$$

According to the Jury conditions [11], to find the asymptotically-stable region of $v_6$, we need to find the region that satisfies the following conditions:

$$P^{*}\left(1\right)>0,P^{*}(-1)<0,\mid C_0\mid <C_n,\mid D_0\mid >\mid D_{n-1}\mid ,$$

where $D_k=\left|\begin{array}{cc}{C}_0& C_{n-k}\\ C_n& C_k\end{array}\right|.$

Since

$P^{*}\left(1\right)=1+C_2+C_1+C_0,$

$P^{*}(-1)=-1+C_2-C_1+C_0,$

Proposition 2.

The fixed point $v_6$ is asymptotically stable if the following conditions are satisfied:

$\mid C_0\mid <1$, $\mid C_0+1\mid >\mid C_1\mid$ and $\mid C_0-1\mid \mid C_0+C_1+1\mid >\mid C_0{C}_1-C_2\mid$.

5. Chaotic Discrete Cancer System

Marotto presented results on mathematical discrete chaos regarding n-dimensional dynamical systems. In the original Marotto theorem, there was an error corrected by Shi and Chen; see [12,13]. In this section, we will prove that system (4) exhibits chaotic dynamics with $h=0.05$ or $h=0.1$; the parameters are the following:

$$\begin{array}{c}{a}_{12}=1,a_{13}=2.5,a_{21}=1.5,a_{31}=0.2,d_3=0.5,\\ k_3=1,r_2=0.6,r_3=4.5.\end{array}$$

(18)

Theorem 1

(Marotto theorem). Consider the following n-dimensional discrete system:

$$v_{n+1}=F\left(v_n\right),n=0,1,2,\dots ,$$

(19)

where $v_n\in {\mathbb{R}}^n$ and $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$ is continuous. Let $B_r\left(v\right)$ denote the ball in ${\mathbb{R}}^n$ of radius r centered at point v and ${\overline{B}}_r\left(v\right)$, its interior. Moreover, let $\Vert v\Vert$ be the usual Euclidean norm of v in ${\mathbb{R}}^n$. Then, $\left(1\right)\Rightarrow \left(2\right)$

(1) 

All eigenvalues of the Jacobian $DF\left(v\right)$ of map (11) at the fixed point v are greater than the one with the Euclidean norm.

(2) 

There exists $s>1$ and $r>0$, such that, for all $u,v\in B_r\left(v\right)$, $\Vert F\left(u\right)-F\left(v\right)\Vert >s\Vert u-v\Vert .$

Theorem 2

(A modified version of the Marotto theorem, [13]). Consider the n-dimensional discrete dynamical system:

$$v_{n+1}=F\left(v_n\right),n=0,1,2,\dots ,$$

(20)

where $v_n\in {\mathbb{R}}^n$ and $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$, suppose that system (12) has a fixed point $v^{*}$.

Assume that

(1) 

F is continuously differentiable in the neighborhood of $v^{*}$; all the eigenvalues of $DF\left(v^{*}\right)$ have absolute values larger than 1, implying that there exists a positive constant r and a Euclidean norm, such that F expands in $B_r\left(v^{*}\right)$ in the Euclidean norm, and

(2) 

$v^{*}$ is a snap-back repeller of F with $F^m\left(v_0\right)=v^{*}$ for some $v_0\in B_r\left(v^{*}\right)$, $v_0\ne v^{*}$ and some positive integer m. Furthermore, F is continuously differentiable in some neighbourhoods of $v_0,v_1,\dots ,v_{m-1}$, respectively, and $det\left[DF\left(v_j\right)\right]\ne 0$ for $0\le j\le m-1$, where $v_j=F\left(v_{j-1}\right)$.

Then, all of the results of the Marotto Theorem hold.

A Proof of the Chaos of the Discrete Cancer System

Step 1. Let $v_2=(0,1,0)$ be the fixed point of system (4).

$F\left(v_2\right)$ given in Theorem 2 of system (4); it is continuously differentiable in $B_r\left(v_2\right)$ for some $r>0$. The Jacobian matrix evaluated at the fixed point $v_2$ is given in (12).

Moreover, (12) has eigenvalues of ${\lambda }_1=0.94$, ${\lambda }_2=1$ and ${\lambda }_3=0.95$.

Step 2. According to Definition (Theorem 2), the snap-back repeller, we need to find one point $u\in B_r\left(v_2\right),$ such that $u\ne v_2$, $F^M\left(u\right)=v_2$ and $det\left[D{F}^M\left(u\right)\right]\ne 0$ for some positive integer M.

In fact, we have

$$\left\{\begin{array}{c}(x(1-x)-a_{12}xy-a_{13}xz).h+x=x_1\\ (r_{2}y(1-y)-a_{21}xy).h+y=y_1\\ (\frac{r_{3}xz}{x+k_3}-a_{31}xz-d_{3}z).h+z=z_1\end{array}\right.$$

(21)

$$\left\{\begin{array}{c}(x_1(1-x_1)-a_{12}{x}_1{y}_1-a_{13}{x}_1{z}_1).h+x=0\\ (r_2{y}_1(1-y_1)-a_{21}{x}_1{y}_1).h+y=1\\ (\frac{r_3{x}_1{z}_1}{x_1+k_3}-a_{31}{x}_1{z}_1-d_3{z}_1).h+z=0\end{array}\right.$$

(22)

Finally, system (4) verifies the conditions of Theorem 2 with the parameters given in (18) and $h=0.1$, the fixed point $v_2$ has two stable and one unstable eigenvalue. Therefore, we have a saddle at this fixed point and there exists a point $u=(-1.1903,0.7563,2.2828)$ solution of (21) and (22), satisfying that $F^2\left(u\right)=v_2$ and $det\left(F\right(u\left)\right)=-6.6158\ne 0$ $det\left(F^2\left(u\right)\right)=27.9025\ne 0$. Thus, $v_2$ is a snap-back repeller.

6. Numerical Simulations

Lyapunov Exponents

In this subsection, we calculate the Lyapunov exponents. The Lyapunov exponents for a discrete n-dimensional systems is given in [14], with the following definition. For other techniques of mathematics agains cancer see for instance [15,16,17]:

Definition 1.

Consider the n-dimensional discrete dynamical system:

$$v_{k+1}=F\left(v_k\right),v_k\in {\mathbb{R}}^n,k=0,1,2,\dots$$

(23)

where $F:{\mathbb{R}}^n\to {\mathbb{R}}^n$ is the vector field associated with map (15), let $J\left(v\right)$ be its Jacobian evaluated at v, also define the matrix: $T_p\left(v_0\right)=J\left(v_{p-1}\right)J\left(v_{p-2}\right)\dots J\left(v_1\right)J\left(v_0\right)$.

Moreover, let $J_i(v_0,l)$ be the modulus of the $ith$ eigenvalue of the $lth$ matrix $T_p\left(v_0\right)$ where $i=1,2,\dots ,n$ and $p=0,1,2,\dots$.

Now, the Lyapunov exponents of n-dimensional discrete time models are defined by: ${\lambda }_i\left(v_0\right)=ln\left({lim}_{p\to +\infty }\left(J_i{(v_0,p)}^{\frac{1}{p}}\right)\right)$.

Therefore, the Lyapunov exponents of system (4) with parameters given in (18) and $h=0.1$ are computed to be ${\lambda }_1=0.95003,{\lambda }_2=-1.0546$ and ${\lambda }_3=-5.6174$. The Lyapunov dimension for system (4) equals the dimension of the space state; that is to say, equal to 3. Because the sum of the Lyapunov exponents is negative ${\lambda }_1+{\lambda }_2+{\lambda }_3<0$.

The Lyapunov exponents of system (4) with parameters given in (18) and $h=0.05$ are computed to be ${\lambda }_1=0.97478,{\lambda }_2=-1.0238$ and ${\lambda }_3=-5.5697$.

If at least one Lyapunov exponent is positive for some control parameter value (18), then system (4) is chaotic at that control parameter, see Figure 11, Figure 12 and Figure 13.

7. Conclusions

This paper contributes to the study of the discrete cancer system with numerical and theoretical methods. As a result of this study, it is clearly understood that the Runge–Kutta method is the best method for discretization of the cancer system. Moreover, the Taylor series expansion method has good accuracy. The Euler discretization method is less accurate but easy to perform. The simulation plots suggest that the cancer system and simulation study reveal that dynamical patterns of the cancer system are dependent on the initial parametric values of the system variables, See Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13. Hence, to obtain the system prediction of a cancer system, accurate quantification of the parametric values of different variables is important. This study could help biologists understand and appreciate the essence of measurement accuracy in different biological experiments and the power of inference through experiments.

In future work, we will study the control of chaos in this cancer system and generalize this system to the fourth-dimension.