Dynamical analysis of fractional order model of immunogenic tumors

Abstract

In this article, we examine the fractional order model of the cytotoxic T lymphocyte response to a growing tumor cell population. We investigate the long-term behavior of tumor growth and explore the conditions of tumor elimination analytically. We establish the conditions for the tumor-free equilibrium and tumor-infection equilibrium to be asymptotically stable and provide the expression of the basic reproduction number. Existence of physical significant tumor-infection equilibrium points is investigated analytically. We show that tumor growth rate, source rate of immune cells, and death rate of immune cells play vital role in tumor dynamics and system undergoes saddle-node and transcritical bifurcation based on these parameters. Furthermore, the effect of cancer treatment is discussed by varying the values of relevant parameters. Numerical simulations are presented to illustrate the analytical results.

Introduction

Mathematical biology is a fast-growing subject which opens up new and exciting branches for mathematicians and biologists. Immune system has natural capacity to detect and destroy abnormal cells which may prevent the growth of many cancers. Mathematical model provides an analytical outline to address which components of the immune system play considerable role in cancer treatment. In 1825, Gompertz used mathematical model first time to explain tumor growth by considering cell replication and death. According to his model, for large populations growth is slow whereas for small populations growth is faster. After that many mathematical models have been developed to understand interactions between the different constituent of the tumor microenvironment. For overviews of simple mathematical models on this topic, see previous works.^1–6 The tumor microenvironment comprises growth factors (containing cytokines and hormones), immune cells, the extracellular matrix, fibroblasts, signaling molecules (chemokines and cytokines), and other connective tissue cells. It is essential to understand these interactions in order to develop the influential cancer immunotherapies.

One of the main open problems in this area is to formulate the mathematical models that capture the dynamics of real data more accurately. In this direction, fractional calculus has made a significant impact on the dynamics of various diseases. ⁷ Fractional derivatives have the excellent property of capturing the memory effects which is detected in nearly all biological systems. This property cannot be covered by means of the integer order derivatives. Bolton et al. ⁷ have shown that fractional order Gompertz model fits the particular tumor growth data set with order 0.68 much better than the ordinary integer order Gompertz growth model. Rodrigues et al. ⁸ notice that dynamics of real outbreak data of dengue cannot be described by classical first-order derivative with a sufficient degree of accuracy. A fractional order model of malaria transmission among heterogeneous populations has been proposed by Pinto and Machado ⁹ to better approximate the real dynamics. HIV fractional complex-order model for drug resistance with three distinct growth rates for the CD4+ T helper cells is analyzed by Pinto and Carvalho. ¹⁰ A mathematical model of dengue fever outbreak is proposed by Diethelm ² using fractional derivative. Recently, Huo et al. ¹¹ studied a fractional order HIV model to assess the impact of vaccines in a homosexual community. They have shown that the vaccinated reproduction below unity is not threshold of HIV eradication when effectiveness and the dosage of the vaccines are low. A new critical threshold is derived in order to eradicate the HIV. Numerous fractional dynamical models have been given by Magin^12–14 for biological systems, for example, membrane charging, vestibular ocular models, one-dimensional cable model for nerve axon, viscoelastic models of cells, and tissues. Stability analysis and bifurcation techniques can help to determine qualitative dynamics of infectious disease models. Recent papers^15,16 provide a good survey of the methods available to analyze the stability of fractional differential equations. In El-Saka et al., ¹⁷ stability, persistence, and Hopf bifurcation of some continuous time fractional order differential equations are studied. Equilibrium points, existence, uniqueness, stability, and numerical solution of fractional order Lotka-Volterra predator prey and rabies system have been studied by Ahmed et al. ¹⁸ Theory of fractional differential equations can be found in Kilbas et al. ¹⁹ and Miller and Ross. ²⁰

The immune reaction to a tumor is generally adjudicated by cells. It is difficult to understand dynamics of the antitumor immune reaction in vivo. Cytotoxic T lymphocytes (CTL) and natural killer (NK) cells perform a leading role in this process. The aim of this article is to investigate fractional order mathematical model of the CTL effect to a flourishing tumor cell population. We examine the long-term behavior of tumor model by varying the values of parameters that might be changed as a result of cancer treatment. Fractional dynamics also exhibit the mechanisms of tumor dormancy. Tumor dormancy is the state when tumor cells persist for long time duration with minor or no change in population of tumor cells. ²¹ Dormant states can emerge after cancer treatment. Dormant states also observe at initial stages of tumor progression. After a long duration, small dormant tumors may escape from immune system and show uncontrolled growth by sneaking through mechanism. ²¹

This article is organized as follows. Formulation of fractional order immune tumor model is given in section “Generalized tumor-immune model.” Existence of physically significant equilibrium points and their stability are investigated in section “Equilibrium points and stability analysis.” Section “The Grünwald–Letnikov Approximation” describes the Grünwald–Letnikov scheme for tumor model. Section “Bifurcation analysis” contains the bifurcation analysis of tumor model. Numerical simulations for possible cases of bifurcation parameters are presented in section “Numerical simulations.” Discussions are given in the last section.

Generalized tumor-immune model

Immune system eliminates the very small tumors before they develop into clinically apparent by reducing their intensification on initial stages. Therefore, mathematical models are very helpful to investigate the tumor dynamics based on interactions between tumor-immune cells. Immune cells are also called effector cells. Mathematical model for the interactions among effector cells and a thriving immunogenic tumor in vivo is given by Kuznetsov et al. ²² using the following coupled system of differential equations

$$\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=s_x(x)+p_x(x,y)-d_x(x,y)-a_x(x),t\ge 0\hfill \\ \frac{\mathrm{dy}}{\mathrm{dt}}=\mathrm{yf}(y)-d_y(x,y),t\ge 0\hfill \end{array}$$

(1)

where x corresponds to the density of the effector cells and y indicates the density of the tumor cells. First equation (1) describes the dynamics of immune cells where

$$p_x(x,y)=\frac{\rho \mathrm{xy}}{\eta +y}(\mathrm{growth\; of\; effector\; cells})$$

$$d_x(x,y)=\mu \mathrm{xy}(\mathrm{death\; rate\; of\; effector\; cells\; by\; tumor\; cells})$$

both of these terms are based on connection between the effector tumor cells

$$a_x(x)=\delta x(\mathrm{natural\; death\; rate\; of\; effector\; cells})$$

$$s_x(x)=\sigma (\mathrm{constant\; treatment\; term})$$

The tumor cells include

$$f(y)=\alpha (1-\beta y)(\mathrm{growth\; rate\; of\; tumor\; cells})$$

$$d_y(x,y)=\mathrm{xy}(\mathrm{death\; rate\; of\; tumor\; cells})$$

Our aim is to examine the dynamics of the following fractional order system

$$\{\begin{array}{c}{D}^{q}x=\sigma +\frac{\rho \mathrm{xy}}{\eta +y}-\mu \mathrm{xy}-\delta x,t\ge 0\hfill \\ D^{q}y=\alpha y(1-\beta y)-\mathrm{xy},t\ge 0\hfill \end{array}$$

(2)

together with the initial conditions $x(0)=x_0\ge 0,y(0)=y_0\ge 0$ , where $D^{q}f,f:{\mathbb{R}}_{+}\to \mathbb{R}$ , represent the Caputo fractional derivative of order $0<q<1$ , defined by

$$D^{q}f(t)=\frac{1}{\Gamma (1-q)}\underset{0}{\overset{t}{\int }}{(t-s)}^{-q}f\prime (s)\mathrm{ds}$$

The first term on the right-hand side represents the source of immune-effector cells. The second term is the recruitment of tumor-specific effector cells which follow the Michaelis–Menten factor to specify the saturated effect of immune response. On the other hand, the effector cells are being eradicated by tumor cells which follow the law of mass action at a proportional rate $\mu$ . The final term of the first equation represents the natural death of effector cells at a rate $\delta$ . Second equation describes the rate of change of tumor population which follow logistic growth where $\alpha$ is the intrinsic growth rate and $1/\beta$ is the biotic capacity. The final term of second equation represents the elimination of tumor cells due to interaction with tumor-specific effector cells which follow the law of mass action. Values of these parameters are given in Table 1.

Table 1.

Parameters values.

Parameters	Description	Values
$\sigma$	Source rate of effector cells	0.1181
$\rho$	Maximum rate of effector cell proliferation	1.131
$\eta$	Half saturation constant	20.19
$\mu$	Effector cells inactivation rate by tumor cells	0.00311
$\delta$	The natural death rate of effector cells	0.3743
$\alpha$	The maximal growth rate of tumor cells	1.636
$\beta$	$\frac{1}{\beta }$ is the maximal carrying capacity	0.002

We denote

$${\mathbb{R}}_{+}^2=\{(x,y)\in {\mathbb{R}}^2;x\ge 0,y\ge 0\}$$

Lemma 2.1

Assume that the function $G:{\mathbb{R}}_{+}\times {\mathbb{R}}^2\to {\mathbb{R}}^2$ satisfies the conditions: ²³

1. $G(t,Y(t))$ is Lebesgue measurable with respect to t on ${\mathbb{R}}_{+}$ ;

2. $G(t,Y(t))$ is continuous with respect to Y on ${\mathbb{R}}^2$ ;

3. $\partial f(t,Y)/\partial Y$ is continuous with respect to Y on ${\mathbb{R}}^2$ ;

4. $\parallel G(t,X)\parallel \le \omega +\lambda \parallel Y\parallel$ , for all $t\in {\mathbb{R}}_{+},Y\in {\mathbb{R}}^2$ .

Then, the initial value problem

$$\{\begin{array}{c}{D}^{q}Y(t)=G(t,Y(t))\hfill \\ \multicolumn{1}{c}{Y(t_0)=Y_0}\end{array}$$

(3)

has a unique solution.

Theorem 2.2

There is a unique solution for the initial value problem (2) and the solution remains in ${\mathbb{R}}_{+}^2$ .

Proof

From Lemma 2.1, we obtain the unique solution on $(0,\infty )$ solving the initial value problem (2). Next, we will show the nonnegative orthant ${\mathbb{R}}_{+}^2$ is a positively invariant region. Since

$$D^{q}x{|}_{x=0}=\sigma ,D^{q}y{|}_{y=0}=0$$

on each line bounding the nonnegative octant, the vector field $(\sigma +((\rho x(t)y(t))/\eta +y(t))-\mu x(t)y(t)-\delta x(t),\alpha y(t)(1-\beta y(t))-x(t)y(t))$ points into ${\mathbb{R}}_{+}^2$ . The solution will remain in ${\mathbb{R}}_{+}^2$ . ²⁴ Hence, the solution $(x(t),y(t))$ of the fractional model (2) is nonnegative if the initial condition is nonnegative for all $t>0$ .

Equilibrium points and stability analysis

To understand the dynamics of (2), first we will find the equilibrium points of the system. Nullclins, that is, the curves along which $D^{q}x=0$ and $D^{q}y=0$ of the system (2) are described by the following:

• $D^{q}x=0\Rightarrow x=\frac{\sigma }{\delta +\mu y-\frac{\rho y}{\eta +y}}$ as long as $\delta +\mu y\ne \frac{\rho y}{\eta +y}$ .

Letting $f(y)$ be the function of tumor population given by

$$f(y)=\frac{\sigma }{\delta +\mu y-\frac{\rho y}{\eta +y}}$$

• $D^{q}y=0\Rightarrow \{\begin{array}{c}y=0\hfill \\ \multicolumn{1}{c}{\mathrm{or}}\\ \multicolumn{1}{c}{x=\alpha (1-\beta y)}\end{array}$

Letting $g(y)$ be function of tumor population given by

$$g(y)=\alpha (1-\beta y)$$

The classification of equilibrium points is as follows:

• Tumor-free equilibrium: $E_0=(\overline{x},\overline{y})=(\sigma /\delta ,0)$

• Tumor-infection equilibrium: $E=(\overline{x},\overline{y})$ , where

$$\overline{x}=g(y),\overline{y}=y$$

(4)

here y is nonnegative solution of

$$f(y)-g(y)=0$$

Depending on the parameter values, there may be from zero to three tumor-infection equilibrium points.

To discuss the local asymptotic stability, we have the following linearized system

$$\{\begin{array}{c}{D}^{\alpha }x=(\frac{\rho \overline{y}}{\eta +\overline{y}}-\mu \overline{y}-\delta )x+(\frac{\rho \overline{x}\eta }{{(\eta +\overline{y})}^2}-\mu \overline{x})y\hfill \\ D^{\alpha }y=-\overline{y}x+(\alpha -2\alpha \beta \overline{y}-\overline{x})y\hfill \end{array}$$

(5)

where $(\overline{x},\overline{y})$ denotes an equilibrium point of the system (2). The Jacobian at $E=(\overline{x},\overline{y})$ is given by

$$J(E)=\left(\begin{array}{cc}\frac{\rho \overline{y}}{\eta +\overline{y}}-\mu \overline{y}-\delta \hfill & \frac{\rho \overline{x}\eta }{{(\eta +\overline{y})}^2}-\mu \overline{x}\hfill \\ \multicolumn{1}{c}{-\overline{y}}& \alpha -2\alpha \beta \overline{y}-\overline{x}\hfill \end{array}\right)$$

The associated transcendental characteristic equation of system (2) is given by $|J(E)-\lambda I|=0$ , where I is the identity matrix. Let

$$\begin{array}{c}{P}_1=\frac{\rho \overline{y}}{\eta +\overline{y}}-\mu \overline{y}-\delta ,P_2=\alpha -2\alpha \beta \overline{y}-\overline{x},\hfill \\ P_3=\overline{y}(\frac{\rho \overline{x}\eta }{{(\eta +\overline{y})}^2}-\mu \overline{x})\hfill \\ m_1=-(P_1+P_2)\hfill \end{array}$$

(6)

and

$$m_2=P_1{P}_2+P_3$$

(7)

Then, we get

$$|J(E)-\lambda I|={\lambda }^2+m_1\lambda +m_2$$

Hence, eigenvalues of the Jacobian matrix are given by

$${\lambda }_{1,2}=\frac{1}{2}(-m_1\pm \sqrt{m_1^2-4m_2})$$

Local stability of tumor-free equilibrium point

In the absence of tumor (i.e. $y=0$ ), the system (2) has a unique tumor-free equilibrium point $E_0=(\overline{x},\overline{y})=(\sigma /\delta ,0)$ . The associated Jacobian of the model at $E_0$ is

$$J(E)=\left(\begin{array}{cc}-\delta \hfill & \frac{\rho \eta \sigma }{\delta {\eta }^2}-\frac{\mu \sigma }{\delta }\hfill \\ \multicolumn{1}{c}{0}& \alpha -\frac{\sigma }{\delta }\hfill \end{array}\right)$$

This gives

$${\lambda }_1=-\delta <0,{\lambda }_2=\alpha -\frac{\sigma }{\delta }$$

Clearly, ${\lambda }_1$ is negative and ${\lambda }_2$ is negative if and only if

$$R_0=\frac{\delta (\alpha +1)-\sigma }{\delta }<1$$

Here, $R_0$ is the basic reproduction number. According to above analysis, we can establish following theorem which gives the condition of tumor eradication.

Theorem 3.2

The tumor-free equilibrium point of the system (2) is locally asymptotically stable if $R_0<1$ and unstable if $R_0>1$ .

Local stability of tumor-infection equilibrium point

The tumor-infection equilibrium corresponds to the case where tumor infection persists $(y\ne 0)$ . These equilibrium points cannot be expressed in closed form, and we derive analytical conditions about their existence. To do this, set $f(y)$ equal to $g(y)$ , this yields a third-order polynomial for the y values of infection equilibrium points

$$A{x}^3+B{x}^2+\mathrm{Cx}+D=0$$

(8)

where

$$A=\mu \beta$$

$$B=-\mu +(\mu \eta +\delta -\rho )\beta$$

$$C=\frac{\sigma }{\alpha }+\rho -\mu \eta -\delta +\delta \eta \beta$$

$$D=\eta (\frac{\sigma }{\alpha }-\delta )$$

The infection equilibrium point of the model (2) can be obtained by substituting the solutions of (8) into (4). Let

$$\Delta =18\mathrm{ABCD}-4B^{3}D+B^2{C}^2-4A{C}^3-27A^2{D}^2$$

Notice that $A>0$ , we have following possibilities under which (8) has positive solution (Tables 2–4).

Table 2.

$\Delta >0$ .

Conditions	Condition on $R_0$	Equilibria of system (2)
$B<0$ , for all $C\in R$	$D>0$	Two distinct positive equilibria
$B\ge 0,C<0$	$D>0$	Two distinct positive equilibria
$B<0,C\ge 0$	$D=0$	Two distinct positive equilibria
$C<0$ , for all $B\in R$	$D=0$	One positive equilibrium
$B\ge 0$ , for all $C\in R$	$D<0$	One positive equilibrium
$B<0,C>0$	$D<0$	Three distinct positive equilibria
$B<0,C\le 0$	$D<0$	One positive equilibrium

Table 3.

$\Delta =0$ .

Conditions	Condition on $R_0$	Equilibria of system (2)
$B<0$ , for all $C\in R$	$D>0$	Two same positive equilibria
$B\ge 0,C<0$	$D>0$	Two same positive equilibria
$B\ge 0,C\ge 0$	$D=0$	One positive equilibrium
$B<0,C\le 0$	$D=0$	One positive equilibrium
$B\mathrm{<}0,C\mathrm{>}0$	$D=0$	Two same positive equilibria
$B\ge 0$ , for all $C\in R$	$D<0$	One positive equilibrium
$B<0,C>0$	$D<0$	Two same positive equilibria
$B<0,C\le 0$	$D<0$	One positive equilibrium

Table 4.

$\Delta <0$ .

Conditions	Condition on $R_0$	Equilibria of system (2)
For any value of $B,C\in R$	$D\ge 0$	No positive equilibrium exists
For any value of $B,C\in R$	$D<0$	One positive equilibrium

Due to $D=(\eta \delta /\alpha )(1-R_0)$ , we can state following useful theorem.

Theorem 3.4

1. If $R_0<1$ and $\Delta >0$ , then system (2) has two distinct positive tumor-infection equilibrium points if and only if (i) $B<0$ or (ii) $B\ge 0$ and $C<0$ .

2. If $R_0<1$ and $\Delta =0$ , then system (2) has two multiple positive tumor-infection equilibrium points if and only if (i) $B<0$ or (ii) $B\ge 0$ and $C<0$ .

3. If $R_0<1$ and $\Delta <0$ , then system (2) has no positive tumor-infection equilibrium point.

4. If $R_0=1$ and $\Delta >0$ , then system (2) has

5. one positive tumor-infection equilibrium points if and only if $C<0$ .

If $R_0=1$ and $\Delta =0$ , then system (2) has

one positive tumor-infection equilibrium points if and only if $B\ge 0$ and $C\ge 0$ , or

one positive tumor-infection equilibrium points if and only if $C\le 0$ , or

two multiple positive tumor-infection equilibrium points if and only if $B<0$ and $C>0$ .

If $R_0=1$ and $\Delta <0$ , then system (2) has no positive tumor-infection equilibrium point.

If $R_0>1$ and $\Delta >0$ , then system (2) has

one positive tumor-infection equilibrium points if and only if $B\ge 0$ or

three distinct positive tumor-infection equilibrium points if and only if $B<0$ and $C>0$ , or

one positive tumor-infection equilibrium points if and only if $B<0$ and $C\le 0$ .

If $R_0>1$ and $\Delta =0$ , then system (2) has

one positive tumor-infection equilibrium points if and only if $B\ge 0$ or

two multiple positive tumor-infection equilibrium points if and only if $C\le 0$ or

one positive tumor-infection equilibrium points if and only if $B<0$ and $C\le 0$ .

If $R_0>1$ and $\Delta <0$ , then system (2) has unique positive tumor-infection equilibrium point.

Using Proposition 1 given by Ahmed and Elgazzar, ²⁵ we can propose the following theorem regarding stability of tumor-infection equilibrium point.

Theorem 3.5

The tumor-infection equilibrium point of the system (2) is locally asymptotically stable if one of the following conditions is satisfied

1. $m_1>0$ and $m_2>0$ ;

2. $m_1<0$ , $4m_2>(m_1{)}^2$ , $\left|\stackrel{-1}{\tan }\left(\frac{\sqrt{4m_2-{(m_1)}^2}}{m_1}\right)\right|>\frac{\alpha \pi }{2}$

where $m_1$ and $m_2$ are defined in equations (6) and (7), respectively.

The Grünwald–Letnikov approximation

Let $0<q<1$ . Assume that the function $u:[0,T]\to R$ satisfies some smoothness conditions on every finite interval $(0,t)$ with $t\le T$ . Choosing an equidistant grid on $[0,t]$

$$0=s_0<s_1<\dots <s_{n+1}=t=(n+1)h$$

with $s_{n+1}-s_n=h$ , the Grünwald–Letnikov derivative is defined by

$$D^{q}u(t)=\underset{h\to 0}{lim}\frac{1}{h^q}{\Delta }_h^{q}u(t)$$

where

$$\frac{1}{h^q}{\Delta }_h^{q}u(t)=\frac{1}{h^q}(u(s_{n+1})-\sum _{k=1}^{n+1}{w}_k^{q}u(s_{n+1-k}))$$

and

$$w_k^q=(-1{)}^{k-1}\left(\begin{array}{c}q\hfill \\ \multicolumn{1}{c}{k}\end{array}\right)$$

These binomial coefficients are recursively connected ²⁶ by

$$w_k^q=(1-\frac{q+1}{k})w_{k-1}^q$$

where the first coefficient is $w_1^q=q$ .

Assume that there exists a unique solution $u=u(s)$ of the following Caputo fractional differential equation in the interval $[0,T]$

$$D^{q}u(t)=f(u(t)),u(s_0)=u_0$$

Let $u_k$ denote the approximation of the true solution $u(s_k)$ . Choose an equidistant grid

$$0=s_0<s_1<\dots <s_{N+1}=T,\mathrm{with}{s}_{k+1}-s_k=h$$

Then, the explicit Grünwald–Letnikov method is given by ²⁶

$$u_{n+1}=h^{q}f(u_n)+\sum _{k=1}^{n+1}{w}_k^{q}u(s_{n+1-k})+r_{n+1}^q{u}_0$$

where $r_{n+1}^q$ is the correction term defined by

$$r_k^q=\frac{t^{k+1}}{\Gamma (k+1-q)}$$

Let $x_k$ and $y_k$ denote the approximation of the true solution $x(s_k)$ and $y(s_k)$ , respectively. Then, using the Grünwald–Letnikov discretization method, we obtain the following scheme for system (2)

$$\{\begin{array}{c}{x}_{n+1}=h^q(\sigma +\frac{\rho x_n{y}_n}{\eta +y_n}-\mu x_n{y}_n-\delta x_n)+\sum _{k=1}^{n+1}{w}_k^q{x}_{n+1-k}+r_{n+1}^q{x}_0\hfill \\ y_{n+1}=h^q(\alpha y_n(1-\beta y_n)-x_n{y}_n)+\sum _{k=1}^{n+1}{w}_k^q{y}_{n+1-k}+r_{n+1}^q{y}_0\hfill \end{array}$$

Bifurcation analysis

In this section, we explore the long-term outcome of the system by using bifurcation analysis. We will investigate the effect of parameters on the tumor dynamics, that is, how the number of equilibria and their stability change by varying the different parameters.

For the parameter values given in Table 1, $R_0=2.3205>1$ , therefore tumor-free equilibrium $E_0=(0.3155,0)$ is unstable. We cannot guess whether the tumor will grow forever or it will reduce in size. We need to find the stability and position of the other equilibrium points in order to predict the long-term behavior of our tumor model. First, we will analyze the effect of variation in parameter $\alpha$ (tumor growth) and keeping all other parameters fixed as given in Table 1. We can see

• If $\alpha \in (0,0.3155]$ , then $R_0<1$ ;

• If $\alpha \in (0.3155,2]$ , then $R_0>1$ .

This implies that tumor-free equilibrium point is stable on $\alpha \in (0,0.3155]$ and it leads to unstable when $\alpha \in (0.3155,2]$ by Theorem 3.2. Therefore, change in stability of tumor-free equilibrium occurs as $\alpha$ passes through the value $0.3155$ . Furthermore, we noticed that

• If $\alpha \in (0.3155,0.8390]$ , then $\Delta <0$ and $D<0$ ;

• If $\alpha \in (0.8390,2]$ , then $\Delta >0,B<0,C>0,D<0$ .

Hence, according to Theorem 3.4 part (9), system (2) has a unique positive tumor-infection equilibrium point when $\alpha \in (0.3155,0.8390]$ and Theorem 3.4 part 7(ii) implies that system (2) has three distinct tumor-infection equilibrium points when $\alpha \in (0.8390,2]$ . Subsequently, there is saddle-node bifurcation at $\alpha \approx 0.8390$ since the number of tumor-infection equilibrium changes from one to three as $\alpha$ passes through this point. Thus, from the above analysis, we deduce that

• For $\alpha \in (0,2]$ , transcritical bifurcation occurs at the point $\alpha \approx 0.3155$ ;

• For $\alpha \in (0.3155,2]$ , there is saddle-node bifurcation at the point $\alpha \approx 0.8390$ .

From Figure 1, we can see that the tumor-free equilibrium is unstable for $\alpha =1.636$ (gray curves) and it will become stable when $\alpha$ is decreased to $0.2$ , hence the patient is ultimately cured.

Figure 1.

Effect of parameter $\alpha =1.636$ (gray curves) and $\alpha =0.2$ (red curves).

With the help of bone marrow transplant, parameter $\sigma$ (source rate of immune cells) might be varied. This therapy would raise stem cell production. Now we change $\sigma$ and keep other parameters fixed as given in Table 1 in order to investigate the outcome of variation in immune cells source rate on tumor dynamics. We get

• $\sigma \in (0,0.1931]$ implies $\Delta <0$ and $D<0$ ;

• $\sigma \in (0.1931,0.6124]$ implies $\Delta >0$ and $D<0,B<0,C>0$ .

Using same argument as in previous discussion, we conclude that there is saddle-node bifurcation at $\sigma \approx 0.1931$ since the number of tumor-infection equilibrium changes from one to three as $\sigma$ passes through this point:

• If $\sigma \in (0.1931,0.6124]$ , then $R_0>1$ ;

• If $\sigma \in (0.6124,2]$ , then $R_0<1$ .

This yields that change in stability of tumor-free equilibrium occurs as $\sigma$ passes through the value $0.6124$ . Precisely, we obtain that

• For $\sigma \in (0,0.6124]$ , transcritical bifurcation occurs at the point $\alpha \approx 0.1931$ ;

• For $\sigma \in (0.1931,2]$ , there is saddle-node bifurcation at the point $\alpha \approx 0.6124$ .

We can see with raise in $\sigma$ , tumor-free equilibrium will become stable. Physically, this represents that addition in immune cells causes eradication of tumor.

The value of $\delta$ (death rate of immune cells) can be varied with cytotoxic chemotherapy. Following the similar analysis given for $\alpha$ , we can show that

• For $\delta \in (0,0.7441]$ , transcritical bifurcation occurs at the point $\alpha \approx 0.0721$ ;

• For $\delta \in (0.0721,2]$ , there is saddle-node bifurcation at the point $\alpha \approx 0.7441$ .

This implies that with the increment in $\delta$ , other equilibrium disappeared and only stable tumor-free equilibrium left, which signifies the suppression of tumor.

Numerical simulations

In this section, we employ the Grünwald–Letnikov method to fractional model (2). In order to verify the analytical results with numerical simulations, first, we will compute the local stability of tumor-free and tumor-infection equilibrium points of fractional order tumor model (2) using Theorem 3.2 and Theorem 3.4, respectively, for selected values of bifurcation parameter $\sigma$ (see Tables 5 and 6).

Table 5.

Local stability of tumor-free equilibrium points based on $\sigma$ .

$\sigma$	$R_0$	Tumor-free equilibrium point	Local stability
0.318	2.0525	(0.5830, 0)	Unstable
0.182	2.3021	(0.3339, 0)	Unstable
0.073	2.5021	(0.1339, 0)	Unstable

Table 6.

Local stability of tumor-infection equilibrium points based on $\sigma$ .

$\sigma$	Tumor-infection equilibrium point	$m_1$	$m_2$	Local stability
0.318	(1.6028, 10.1408)	0.2316	0.3595	Stable
0.182	(1.5880, 14.6738)	0.1626	0.3708	Stable
–	(0.9074, 222.6896)	0.9293	−0.4041	Unstable
–	(0.2262, 430.8711)	2.2145	0.8422	Stable
0.073	(1.5712, 19.8117)	0.1113	0.3504	Stable
–	(1.0728, 172.1251)	0.6312	−0.4219	Unstable
–	(0.0776, 476.2981)	2.4998	1.3555	Stable

Figure 2 shows the solution of the system (2) for the parameter values given in Table 1. We obtained convergence in a shorter time for fractional order $0.5$ and $0.7$ as compared to $q=1$ . Figure 3 indicates that transients beginning with the initial condition $(0.5,50)$ and $(1.5,460)$ approach to dormant tumor steady state $E_1$ for $q=0.9$ . However, with the slight change in the initial conditions $(0.1,60)$ and $(0.6,470)$ , system (2) approaches to uncontrolled tumor steady state $E_3$ .

Figure 2.

Solution of the system (2). Upper panel represents simulation results of effector cells x with initial condition $(x,y)=(0.5,50)$ (blue curves) and $(x,y)=(0.1,60)$ (green curves). Lower panel represents simulation results of tumor cells y with initial condition $(x,y)=(0.5,50)$ (red curves) and $(x,y)=(0.1,60)$ (gray curves).

Figure 3.

Phase portraits of effector cells x and tumor cells y with $(\delta ,\sigma )=(0.3743,0.1181)$ for $q=0.9$ . Labeling: stable equilibrium point is represented by green circles (•), saddle equilibrium point with green diamonds $(♢)$ and red curve corresponds to one-dimensional stable manifold of equilibrium $E_2$ .

The approximate solutions of system (2) for fractional order $q=0.5,q=0.7,q=0.8$ and $q=1$ are displayed in Figures 4–8 showing the effect of variation in $\sigma$ . At $\sigma =0.318$ , there exists one unstable tumor-free equilibrium point $E_0=(0.583,0)$ (see Table 5) and one stable tumor-infection equilibrium point $E_1=(1.6028,10.1408)$ (see Table 6). Figure 4 shows that effector cells $x(t)\to 1.6028$ and $y(t)\to 10.1408$ as $t\to \infty$ , with different initial conditions $(0.1,370)$ and $(3,270)$ . Thus, our numerical solution approaches the stable equilibrium.

Figure 4.

Solution of the system (2) with $\sigma =0.318$ . Upper panel represents simulations of effector cells with initial condition $(0.1,370)$ (blue curves) and $(3,270)$ (green curves) and lower panel represents simulations of tumor cells with initial condition $(0.1,370)$ (red curves) and $(3,270)$ (gray curves).

Figure 5.

Solution of the system (2) with $\sigma =0.182$ . Upper panel represents simulations of effector cells and lower panel represents simulations of tumor cells with initial condition $(0.6,50)$ .

Figure 6.

Solution of the system (2) with $\sigma =0.182$ . Upper panel represents simulations of effector cells and lower panel represents simulations of tumor cells with initial condition $(0.3,110)$ .

Figure 7.

Solution of the system (2) with $\sigma =0.073$ . Upper panel represents simulations of effector cells and lower panel represents simulations of tumor cells with initial condition $(3.5,130)$ .

Figure 8.

Solution of the system (2) with $\sigma =0.073$ . Upper panel represents simulations of effector cells and lower panel represents simulations of tumor cells with initial condition $(0.1,130)$ .

Table 5 shows that when $\sigma =0.182$ , there is one tumor-free equilibrium point $E_0=(0.3339,0)$ (unstable) and Table 6 depicts that there are three tumor-infection equilibrium points $E_1=(1.5880,14.6738)$ (stable), $E_2=(0.9074,222.6896)$ (unstable), and $E_3=(0.2262,430.8711)$ (stable). Figures 5 and 6 show that effector cells x and the tumor cells BCL₁ leukemia y settle to the “dormant tumor” steady state $E_1=(1.5880,14.6738)$ with initial condition $(0.6,50)$ . However, with the increment in tumor cell, approximate solutions converge to “uncontrolled tumor state” $E_3=(0.2262,430.8711)$ with the initial condition $(0.3,110)$ . On the other hand, Figures 7 and 8 demonstrate that effector cells x and the tumor cells y approach the “dormant tumor” steady state $E_1=(1.5712,19.8117)$ with initial condition $(3.5,130)$ . But approximate solutions converge to “uncontrolled tumor state” $E_3=(0.0776,476.2918)$ by decreasing effector cells and keeping same number of tumor cells, that is, with the initial condition $(0.1,130)$ .

Discussion and conclusion

In this article, we have studied an immunogenic tumor model of fractional order. Dynamics of the model depends on the stability of the equilibrium points:

1. Tumor-free equilibrium: $E_0=(\overline{x},\overline{y})=(\sigma /\delta ,0)$ ;

2. Tumor-infection equilibrium: $E=(\overline{x},\overline{y})$ .

Medically, if the trajectories of the fractional order model will converge to tumor-free equilibrium, then tumor will eliminate in the long run (hence the patient is ultimately cured) and if trajectories will converge to the tumor-infection equilibrium point, then tumor will persist in the body. Mathematically, in order to find cancer treatment, we have to find the analytical conditions under which trajectories of the system will converge to tumor-free equilibrium. This analytical condition is given in Theorem 3.2, that is, the tumor-free equilibrium point of the system is locally asymptotically stable if

$$R_0=\frac{\delta (\alpha +1)-\sigma }{\delta }<1$$

This relation depends on three parameters:

1. $\alpha$ (tumor growth rate);

2. $\sigma$ (source rate of effector cells);

3. $\delta$ (death rate of effector cells).

Precisely, variation of these three parameters plays a key role in cancer treatment and doctors have to adjust these parameters according to above relation in order to eliminate tumor. The parameters $\alpha ,\sigma$ , and $\delta$ might be changed as a consequence of cancer treatment. For instance, effect of variation in $\alpha$ is given in Figure 1. By decreasing tumor growth $\alpha$ , all tumor-infection equilibrium points are disappeared and tumor-free equilibrium will become stable which indicates the situation of tumor eradication and the patient is ultimately cured. The value of $\sigma$ can be increased by bone marrow transplant, while the cytotoxic chemotherapy can raise the value of $\delta$ . With the help of bifurcation analysis, we have shown that increase in $\sigma$ can control the tumor growth and tumor-free equilibrium will become stable. Effect of variations in $\sigma$ is presented in Figures 4–8 for different values of fractional derivative. Another important issue from medical point of view is to investigate the factors which cause growth of tumor. Mathematically, this means, we need to investigate the analytical relations for existence and stability of tumor-infection equilibrium. This matter is treated in Theorem 3.4 which gives all the possible analytical conditions for existence of physically significant tumor-infection equilibrium point and Theorem 3.5 describes the criteria for stability of tumor-infection equilibrium point.

Numerical simulations are presented for various values of fractional order q indicating that antitumor reactivity might be due to cytotoxic effector cells. Fractional calculus may improve the understanding of biological processes because the fractional differential equations can capture past evolution of the function. The memory kernel in the solution of fractional differential equations decays as a power law whereas in ordinary differential equations the memory kernel turns into the Dirac Delta function. Mathematical models can help the physicians to choose optimal dosage. Cancer patients have different tendencies during tumor growth which is difficult to capture by integer order derivative. Fractional derivative can be varied to best fit the real data according to progression of different tumors. Thus, more reliable model can be obtained by choosing the relevant fractional index according to real data which can help clinician to suggest treatment to each individual patient.