Mathematical Modeling and Analysis of Tumor Chemotherapy

Abstract

Cancer diseases lead to the second-highest death rate all over the world. For treating tumors, one of the most common schemes is chemotherapy, which can decrease the tumor size and control the progression of cancer diseases. To better understand the mechanisms of chemotherapy, we developed a mathematical model of tumor growth under chemotherapy. This model includes both immune system response and drug therapy. We characterize the symmetrical properties and dynamics of this differential equation model by finding the equilibrium points and exploring the stability and symmetry properties in a range of model parameters. Sensitivity analyses suggest that the chemotherapy drug-induced tumor mortality rate and the drug decay rate contribute significantly to the determination of treatment outcomes. Numerical simulations highlight the importance of CTL activation in tumor chemotherapy.

1. Introduction

The tumor is a complex disease, which is the second leading cause of death in the world [1,2]. The death rate caused by cancer has been increased worldwide in recent years [3]. Although substantial medical research studies have been conducted, more research studies are strongly needed to understand the mechanism of cancer initiation and destruction as well as the effects of the anti-tumor therapy. Mathematical models provide a valuable theoretical framework through which immuno-oncological mechanisms can be discovered and clarified [4,5,6,7,8]. The applications of mathematical models enable us to study the dynamics and symmetrical property of the immune system as well as to discover the mechanisms controlling tumor cells, to design the therapy schemes, and to predict the outcomes of new therapy strategies [9,10,11,12,13,14,15,16,17,18,19]. In addition, theoretical analysis, modeling methodologies, and symmetrical analysis techniques have been applied to investigate various properties of mathematical models using differential equations [20,21,22,23].

Medical practices show that surgery, chemotherapy, immunotherapy, and radiation therapy so far are the widely used clinical schemes to treat cancer diseases successfully [11,12,13,14,15,16,17,18]. However, the treatment outcome depends on a number of factors, including the type and degree of the tumor, the implementation of treatment schemes, and the intensity of the patient’s autoimmune response. Mathematical modeling provides a theoretical framework to predict the long-term treatment outcomes, which is difficult to be realized by clinical studies. Currently, there are a variety of mathematical models for tumor growth and treatment. However, each model only studies a couple of important factors in cancer diseases and the medical treatment effects [14,15,16,17,18,24,25,26,27,28,29,30,31,32].

Among the treatment schemes for tumor diseases, the most common one is chemotherapy. In the work of A.G. López et al. [11], they proposed a dynamical model of tumor growth, which includes the interactions and regulations between tumor cells, healthy host cells, and immune effector cells. By analyzing the dynamical properties of the model, they demonstrated the consistency between the theoretical analysis and empirical discovery regarding the tumor progression and immune function. They introduced a chemotherapeutic treatment scheme and proposed an approach to overcome problems in complex pharmacokinetic modeling. A follow-up study proposed a differential equations model for the interactions of tumor cells and immune cells under chemotherapy [13]. This study suggested that intense chemotherapies could result in inferior tumor control. In addition, more intense therapies might result in a stronger depletion of immune cells, which allows an early re-growth of the tumor. However, these studies neither accounted for the different cell populations of innate and specific immune responses nor contained complex pharmacokinetics.

Motivated by the need to better understand the mechanisms of chemotherapy treatment and find more efficacious therapeutic methods, we propose a mathematical model to study the effects of the immune system and chemotherapy to treat cancer diseases. The purpose of this study is to explore the combined beneficial effects of the immune system to regulate tumor growth and to find the side effects of chemotherapeutic drugs on both the tumor cells and effector cells. Compared to previous models, the main innovation of this work is the combination of important functions of the immune system and the representation of dose-response dynamics in terms of mass action rather than exponential decay. Thus, the proposed model describes three types of cells and one drug concentration. We first analyze the local geometric and symmetrical properties of the equilibria of the proposed model and then use simulations to confirm the theoretical results. In addition, sensitivity analysis finds important parameters that have more impact on the variations of tumor cell numbers, which provide predictions for designing treatment schemes.

2. The Model

2.1. Mathematical Model

We first propose a mathematical model to describe the interaction of tumor cells and the anti-tumor immune system, as well as the influence of chemotherapy on both tumor cells and the immune system. The model variables are $T\left(t\right)$ for tumor cell population at time t, $N\left(t\right)$ for NK cell population, $L\left(t\right)$ for cytotoxic T cell (CTLs) population, and $u\left(t\right)$ for the amount of drug at the tumor site. Figure 1 gives the schematic diagram for the interactions of these variables. The proposed model is based on the following assumptions.

• Both immune effector cells and chemotherapy decrease the tumor population.
• The population of effector cells decreases due to the degradation process, consumption when killing tumor cells, and the effect of chemotherapy.
• Chemotherapy drugs can affect tumor cells and immune effector cells through a mass-action mechanism.
• A higher constant input of the drug dose can result in both higher tumor and immune effector cell depletion.

The proposed model describes the population dynamics of tumor cells and the immune system under chemotherapy, given by

$$\left\{\begin{array}{c}{N}^{{}^{\prime }}\left(t\right)=aN\left(t\right)(1-bN\left(t\right))-{\alpha }_{1}N\left(t\right)T\left(t\right)-k_{N}u\left(t\right)N\left(t\right),\hfill \\ L^{{}^{\prime }}\left(t\right)=rN\left(t\right)T\left(t\right)-\mu L\left(t\right)-{\beta }_{1}L\left(t\right)T\left(t\right)-k_{L}u\left(t\right)L\left(t\right),\hfill \\ T^{{}^{\prime }}\left(t\right)=cT\left(t\right)(1-dT\left(t\right))-{\alpha }_{2}N\left(t\right)T\left(t\right)-{\beta }_{2}L\left(t\right)T\left(t\right)-k_{T}u\left(t\right)T\left(t\right),\hfill \\ u^{{}^{\prime }}\left(t\right)=v-\omega u\left(t\right),\hfill \end{array}\right.$$

(1)

with initial conditions $N\left(0\right)=N_0\ge 0$, $L\left(0\right)=L_0\ge 0$, $T\left(0\right)=T_0\ge 0$, and $u\left(0\right)=u_0\ge 0$.

In the first equation of model (1), it is assumed that NK cells grow logistically through the term $aN\left(t\right)(1-bN(t\left)\right)$. However, they are inactivated through the interactions with tumor cells based on $-{\alpha }_{1}N\left(t\right)T\left(t\right)$. In the second equation, CTLs are present in the system when the tumor cells are present. They are recruited by tumor cells through a linear term $rN\left(t\right)T\left(t\right)$. In addition, the death of CTLs is a linear process $-\mu L\left(t\right)$. In addition, the interactions with tumor cells will inactivate CTLs via $-{\beta }_{1}L\left(t\right)T\left(t\right)$. In the third equation, tumor is assumed to grow according to a logistic function as $cT\left(t\right)(1-dT(t\left)\right)$ [12,13,14,18,24,25]. In addition, the tumor cells are killed by both NK cells and CTLs, which are realized by $-{\alpha }_{2}N\left(t\right)T\left(t\right)$ and $-{\beta }_{2}L\left(t\right)T\left(t\right)$, respectively. The chemotherapy drug in the last equation has a constant source v, and linearly fades out of the system by $-\omega u\left(t\right)$. In our model (1), since the chemotherapy drug affects all three cell populations through a mass-action dynamic and the drug mortality rate differs for each type of cell, we denote the three different response coefficients by $k_N$, $k_L$, and $k_T$.

Table 1. Description and values of the parameters in model (1).

Parameters	Units	Description	Value	Reference
a	day${}^{-1}$	Growth rate of NK cells	none	none
b	cell${}^{-1}$	Inverse of NK cells capacity	$3.17\times {10}^{-6}$	fitting
c	day${}^{-1}$	Growth rate of tumor	$5.14\times {10}^{-1}$	[33]
d	cell${}^{-1}$	Inverse of tumor capacity	$1.02\times {10}^{-9}$	[33]
r	cell${}^{-1}$ day${}^{-1}$	Activation rate of CTLs	$1.1\times {10}^{-7}$	[34,35]
$\mu$	day${}^{-1}$	CTL death rate	$2.0\times {10}^{-2}$	[34]
${\alpha }_1$	cell${}^{-1}$ day${}^{-1}$	NK cell death rate	$1.0\times {10}^{-7}$	[33]
${\alpha }_2$	cell${}^{-1}$ day${}^{-1}$	Tumor death rate of NK	$6.41\times {10}^{-11}$	[33,36]
${\beta }_1$	cell${}^{-1}$ day${}^{-1}$	CTL death rate	$3.42\times {10}^{-10}$	[37]
${\beta }_2$	cell${}^{-1}$ day${}^{-1}$	Rate of CTL-induced tumor death	$3.5\times {10}^{-7}$	fitting
v	dose	Influx of drug	none	none
$\omega$	day${}^{-1}$	Drug decay rate	$9\times {10}^{-1}$	[12]
$k_N,k_L$	day${}^{-1}$	Immune cell killed by drug	$6\times {10}^{-1}$	[12]
$k_T$	day${}^{-1}$	Tumor cell killed by drug	$8\times {10}^{-1}$	[12]

lists the parameters in model (1), along with their units, descriptions, estimated values, and references.

2.2. The Reduced Model

To reduce the model complexity, we use the following non-dimensionalized state variables to simplify the model, given by

$$\overline{N}=\frac{{\alpha }_2}{\mu }N,\overline{L}=\frac{{\alpha }_1{\alpha }_2}{r\mu }L,\overline{T}=\frac{{\alpha }_1}{\mu }T,\overline{u}=\frac{k_N}{\mu },dt=\frac{1}{\mu }d\tau .$$

Then the simplified model is

$$\left\{\begin{array}{c}{\overline{N}}^{{}^{\prime }}\left(\tau \right)=\frac{a}{\mu }\overline{N}\left(\tau \right)(1-\frac{b\mu }{{\alpha }_2}\overline{N}\left(\tau \right))-\overline{N}\left(\tau \right)\overline{T}\left(\tau \right)-\overline{u}\left(\tau \right)\overline{N}\left(\tau \right),\hfill \\ {\overline{L}}^{{}^{\prime }}\left(\tau \right)=\overline{N}\left(\tau \right)\overline{T}\left(\tau \right)-\overline{L}\left(\tau \right)-\frac{{\beta }_1}{{\alpha }_1}\overline{L}\left(\tau \right)\overline{T}\left(\tau \right)-\overline{u}\left(\tau \right)\overline{L}\left(\tau \right),\hfill \\ {\overline{T}}^{{}^{\prime }}\left(\tau \right)=\frac{c}{\mu }\overline{T}\left(\tau \right)(1-\frac{d\mu }{{\alpha }_1}\overline{T}\left(\tau \right))-\overline{N}\left(\tau \right)\overline{T}\left(\tau \right)-\frac{r{\beta }_2}{{\alpha }_1{\alpha }_2}\overline{L}\left(\tau \right)\overline{T}\left(\tau \right)-\frac{k_T}{k_N}\overline{u}\left(\tau \right)\overline{T}\left(\tau \right),\hfill \\ {\overline{u}}^{{}^{\prime }}\left(\tau \right)=\frac{k_{N}v}{{\mu }^2}-\frac{\omega }{\mu }\overline{u}\left(\tau \right).\hfill \end{array}\right.$$

(2)

By rewriting $\overline{N},\overline{L},\overline{T},\overline{u},\tau$ as $N,L,T,u$, and t, respectively, we have the following dimensionless model, given by

$$\left\{\begin{array}{ccc}{N}^{{}^{\prime }}\left(t\right)\hfill & =\hfill & pN\left(t\right)(1-1.8\times {10}^{-2}N\left(t\right))-N\left(t\right)T\left(t\right)-u\left(t\right)N\left(t\right),\hfill \\ L^{{}^{\prime }}\left(t\right)\hfill & =\hfill & N\left(t\right)T\left(t\right)-L\left(t\right)-3.42\times {10}^{-3}L\left(t\right)T\left(t\right)-u\left(t\right)L\left(t\right),\hfill \\ T^{{}^{\prime }}\left(t\right)\hfill & =\hfill & 25.7T\left(t\right)(1-2.04\times {10}^{-4}T\left(t\right))-N\left(t\right)T\left(t\right)\hfill \\ \hfill & \hfill & -6.02\times {10}^{3}L\left(t\right)T\left(t\right)-ku\left(t\right)T\left(t\right),\hfill \\ u^{{}^{\prime }}\left(t\right)\hfill & =\hfill & s-fu\left(t\right),\hfill \end{array}\right.$$

(3)

where $p=\frac{a}{\mu }$, $k=\frac{k_T}{k_N}$, $s=\frac{k_{N}v}{{\mu }^2}$, $f=\frac{\omega }{\mu }$. Using the values in Table 1, the values in model (3) are $\frac{b\mu }{{\alpha }_2}=1.8\times {10}^{-2}$, $\frac{{\beta }_1}{{\alpha }_1}=3.42\times {10}^{-3}$, $\frac{c}{\mu }=2.57\times 10$, $\frac{d\mu }{{\alpha }_1}=2.04\times {10}^{-4}$, and $\frac{r{\beta }_2}{{\alpha }_1{\alpha }_2}=6.02\times {10}^3$.

3. Dynamics

This section will first find the positive invariant set of model (3). Then we study the existence of equilibrium states and their corresponding stability properties. We give the first result for the positive invariant set.

Lemma 1.

The solutions of model (3) are positive if all initial values are positive and $R_{+}^4=\left\{(N\left(t\right),L\left(t\right),T\left(t\right),u\left(t\right))\right|N\left(t\right)\ge 0,L\left(t\right)\ge 0,T\left(t\right)\ge 0,u\left(t\right)\ge 0\}$ is a positive invariant set of model (3).

Proof. 

Let $X\left(t\right)=\left(N\right(t),L(t),T(t),u(t\left)\right)$ be a positive solution of model (3) based on the initial solution $X\left(0\right)=(N_0,L_0,T_0,u_0)\in R_{+}^4$. Following the dimensionless model (3), we have the following equations

$$\begin{array}{c}N\left(t\right)=N_0{e}^{{\int }_0^t(p(1-1.8\times {10}^{-2}N\left(m\right))-T\left(m\right)-u\left(m\right))dm},\hfill \\ L^{{}^{\prime }}\left(t\right)\ge -L\left(t\right)-3.42\times {10}^{-3}L\left(t\right)T\left(t\right)-u\left(t\right)L\left(t\right),\hfill \\ T\left(t\right)=T_0{e}^{{\int }_0^t(2.57\times 10(1-2.04\times {10}^{-4}T\left(m\right))-N\left(m\right)-6.02\times {10}^{3}L\left(m\right)-ku\left(m\right))dm},\hfill \\ u\left(t\right)=u_0{e}^{{\int }_0^t(\frac{s}{u\left(m\right)}-f)dm}.\hfill \end{array}$$

(4)

If $N_0\ge 0,T_0\ge 0$, and $u_0\ge 0$, we can derive that $N\left(t\right)\ge 0,T\left(t\right)\ge 0,u\left(t\right)\ge 0$. In addition, the second inequality in (4) gives

$$\begin{array}{c}\hfill {\displaystyle L\left(t\right)\ge L_0{e}^{{\int }_0^t(-1-3.42\times {10}^{-3}T\left(m\right)-u\left(m\right))dm}.}\end{array}$$

(5)

Using the comparison principle, we can obtain that, if $L_0\ge 0$, then $L\left(t\right)\ge 0$. Thus, if $N_0\ge 0,,T_0\ge 0,u_0\ge 0$, then $N\left(t\right)\ge 0,L\left(t\right)\ge 0,T\left(t\right)\ge 0,u\left(t\right)\ge 0$. That is, all solutions of model (3) are positive when their initial values are positive. Thus, $R_{+}^4$ is a positively invariant set of model (3). □

3.1. Equilibria of Dimensionless Model

We first note that the fourth equation in (3) can be decoupled from the model. The equilibrium of variable u is $u^{\ast }=\frac{s}{f}$. Thus, we need to consider the model with the remaining three equations,

$$\left\{\begin{array}{ccc}{N}^{{}^{\prime }}\left(t\right)\hfill & =\hfill & pN\left(t\right)(1-1.8\times {10}^{-2}N\left(t\right))-N\left(t\right)T\left(t\right)-u\left(t\right)N\left(t\right)\hfill \\ L^{{}^{\prime }}\left(t\right)\hfill & =\hfill & N\left(t\right)T\left(t\right)-L\left(t\right)-3.42\times {10}^{-3}L\left(t\right)T\left(t\right)-u\left(t\right)L\left(t\right)\hfill \\ T^{{}^{\prime }}\left(t\right)\hfill & =\hfill & 25.7T\left(t\right)(1-2.04\times {10}^{-4}T\left(t\right))-N\left(t\right)T\left(t\right)\hfill \\ \hfill & \hfill & -6.02\times {10}^{3}L\left(t\right)T\left(t\right)-ku\left(t\right)T\left(t\right).\hfill \end{array}\right.$$

(6)

We set each equation equals to zero. The null-surfaces of model (6) are given by

$$\begin{array}{c}{N}^{{}^{\prime }}\left(t\right)=0⟹N=0orN=\frac{p-T-u}{1.8\times {10}^{-2}p},\hfill \\ L^{{}^{\prime }}\left(t\right)=0⟹L=\frac{(p-T-u)T}{1.8\times {10}^{-2}p(1+3.42\times {10}^{-3}T+u)},\hfill \\ T^{{}^{\prime }}\left(t\right)=0⟹T=0orT=1.91\times {10}^2(25.7-N-6.02\times {10}^{3}L-ku).\hfill \end{array}$$

(7)

where $u=\frac{s}{f}$. Then the types of equilibrium states can be given as follows.

(i) The dead state, at which all three cell populations are zero, namely

$$E_0(0,0,0,\frac{s}{f})$$

(ii) The tumor-free state, at which the population of tumor cells is zero but the NK cells survive. The equilibrium point is given by

$$E_1(\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p},0,0,\frac{s}{f}),whenp>\frac{s}{f}.$$

(iii) The tumor-present state, at which the populations of NK cells and CTLs are zero but the tumor cells survive.

$$E_2(0,0,1.91\times {10}^2(25.7-\frac{ks}{f}),\frac{s}{f}),when25.7>\frac{ks}{f}.$$

(iv) The coexisting state, at which immune and tumor cells have nonzero populations. There are the following possible equilibrium states $E_3(N_i^{\ast },L_i^{\ast },T_i^{\ast },\frac{s}{f})$,$(i=1,2,3)$, where $T_i^{\ast }$ is a positive solution to the following equation

$$A{T}^2+BT+C=0,$$

(8)

where

$$\begin{array}{cc}\hfill A& =f^2(9.44\times {10}^{-3}-3.23\times {10}^{-7}p),\hfill \\ \hfill B& =f^2+1.01fs-7.95\times {10}^{-3}{f}^{2}p-9.43\times {10}^{-5}fsp-6.16\times {10}^{-5}kspf,\hfill \\ \hfill C& =fs+s^2-5.37\times {10}^{-1}{f}^{2}p-5.37\times {10}^{-1}fsp-1.8\times {10}^{-2}kspf-1.8\times {10}^{-2}k{s}^{2}p.\hfill \end{array}$$

The coexisting state can also be classified into two cases.

(1) When $A=0$, that is, $p=2.92\times {10}^4$, we have positive equilibrium $E_3(N_1^{\ast },L_1^{\ast },T_1^{\ast },\frac{s}{f})$. Here

$$\begin{array}{cc}\hfill N_1^{\ast }& =1.9\times {10}^{-3}(2.92\times {10}^4-T^{\ast }-\frac{s}{f}),\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill L_1^{\ast }& =\frac{1.9\times {10}^{-3}(2.92\times {10}^4-T^{\ast }-\frac{s}{f})T^{\ast }}{1+3.42\times {10}^{-3}{T}^{\ast }+\frac{s}{f}},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill T_1^{\ast }& =\frac{1.57\times {10}^4{f}^2+s^2-1.57\times {10}^{4}fs-5.26\times {10}^{2}ksf-5.26\times {10}^{2}k{s}^2}{f(2.31\times {10}^{2}f+1.74s+1.8ks)}.\hfill \end{array}$$

(11)

When $0<T_1^{\ast }<2.92\times {10}^4-\frac{s}{f}$, we have that $N_1^{\ast }>0$, $L_1^{\ast }>0$. Therefore, a positive equilibrium $E_3(N_1^{\ast },L_1^{\ast },T_1^{\ast },\frac{s}{f})$ exists if and only if $k>5.23\times {10}^{-3}$ and $f>f_1^{\ast }$.

(2) When $A\ne 0$, we have a unique positive equilibrium $E_3(N_2^{\ast },L_2^{\ast },T_2^{\ast },\frac{s}{f})$ or $E_3(N_3^{\ast },L_3^{\ast },$$T_3^{\ast },\frac{s}{f})$. Here

$$\begin{array}{cc}\hfill N_i^{\ast }& =\frac{p-T_i^{\ast }-\frac{s}{f}}{1.8\times {10}^{-2}p},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill L_i^{\ast }& =\frac{(p-T_i^{\ast }-\frac{s}{f})T_i^{\ast }}{1.8\times {10}^{-2}p(1+3.42\times {10}^{-3}{T}_i^{\ast }+\frac{s}{f})},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill T_i^{\ast }& =\frac{-B\pm \sqrt{B^2-4AC}}{2A},(i=2,3).\hfill \end{array}$$

(14)

It is clear that, when $p-T_i^{\ast }-\frac{s}{f}>0$, we have $N_i^{\ast }>0$, $L_i^{\ast }>0$, $i=2,3$.

Thus, the coexisting equilibrium states exist if and only if $0<T_1^{\ast }<2.92\times {10}^4-\frac{s}{f}$, or $0<T_i^{\ast }<p-\frac{s}{f}$, $i=2,3$. Based on Vieta’s theorem, the endemic equilibrium states of model (3) are given in

Table 2. Endemic equilibrium states of model (3) and the corresponding parameter ranges.

No.	p	k	f	Equilibrium Points
1	$2.92\times {10}^4$	$(5.23\times {10}^{-3},$ +$\infty )$	$(f_1^{\ast },$  +$\infty )$	$E_3(N_1^{\ast },L_1^{\ast },T_1^{\ast },\frac{s}{f})$
2	$(1.23\times {10}^4,2.92\times {10}^4)$	$(0,\frac{1}{1.8\times {10}^{-2}p})$	$(0,f_3^{\ast })$	$E_3(N_2^{\ast },L_2^{\ast },T_2^{\ast },\frac{s}{f})$
3	$(0,5.06\times {10}^3)$	$(0,k_1^{\ast })$	$(f_2^{\ast },$  +$\infty )$	$E_3(N_3^{\ast },L_3^{\ast },T_3^{\ast },\frac{s}{f})$

.

The values of parameters in Table 2 are given by

$$\begin{array}{c}{f}_1^{\ast }=0.5s+1.68\times {10}^{-2}ks+3.18\times {10}^{-5}\sqrt{s^2(2.46\times {10}^8+2.77\times {10}^5{k}^2+4.95\times {10}^{7}k)},\hfill \\ f_2^{\ast }=\frac{-B_1+\sqrt{B_1^2-4A_1{C}_1}}{2A_1},\hfill \\ f_3^{\ast }=\frac{-B_2-\sqrt{B_2^2-4A_2{C}_2}}{2A_2},\hfill \\ k_1^{\ast }=\frac{0.99-9.37\times {10}^{-5}p}{6.16\times {10}^{-5}p},\hfill \\ A_1=4.17\times {10}^{-13}{p}^3+5.55\times {10}^{-5}p-1.41\times {10}^{-8}{p}^2+2.62\times {10}^{-2},\hfill \\ B_1=1.21\times {10}^{-10}s{p}^2-4.13\times {10}^{-6}sp+1.73\times {10}^{-2}s\hfill \\ +7.96\times {10}^{-11}ks{p}^2-6.8\times {10}^{-4}ks-2.3\times {10}^{-6}ksp,\hfill \\ C_1=-s^2(1.1\times {10}^{-10}p+6.78\times {10}^{-4}k),\hfill \\ A_2=-5.37\times {10}^{-1}p,\hfill \\ B_2=s(1-5.37\times {10}^{-1}p-1.8\times {10}^{-2}kp),\hfill \\ C_2=s^2(1-1.8\times {10}^{-2}kp).\hfill \end{array}$$

Based on the above discussion,

Table 3. The existence of equilibrium states of model (3) ($s\in (0,+\infty )$).

No.	p	k	f	Equilibrium Points
1	$(0,+\infty )$	$(0,+\infty )$	$(0,+\infty )$	$E_0(0,0,0,\frac{s}{f})$
2	$(\frac{s}{f},+\infty )$	$(0,+\infty )$	$(0,+\infty )$	$E_1(\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p},0,0,\frac{s}{f})$
3	$(0,+\infty )$	$(0,+\infty )$	$(\frac{ks}{25.7},+\infty )$	$E_2(0,0,1.91\times {10}^2(25.7-\frac{ks}{f}),\frac{s}{f})$
4	$2.92\times {10}^4$	$(5.23\times {10}^{-3},+\infty )$	$(f_1^{\ast },+\infty )$	$E_3(N_1^{\ast },L_1^{\ast },T_1^{\ast },\frac{s}{f})$
5	$(1.23\times {10}^4,2.92\times {10}^4)$	$(0,\frac{1}{1.8\times {10}^{-2}p})$	$(0,f_3^{\ast })$	$E_3(N_2^{\ast },L_2^{\ast },T_2^{\ast },\frac{s}{f})$
6	$(0,5.06\times {10}^3)$	$(0,k_1^{\ast })$	$(f_2^{\ast },+\infty )$	$E_3(N_3^{\ast },L_3^{\ast },T_3^{\ast },\frac{s}{f})$

gives the existence of the equilibrium states of model (3). Based on different parameter values, model (3) may have different numbers of equilibria from these possible states. When developing a chemotherapy scheme, the goal is to let the model reach a region in which there are only harmless equilibrium states, which may be either the tumor-free equilibrium state at $E_1(\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p},0,0,\frac{s}{f})$ or the coexisting equilibrium state with a small number of tumor cells but a large number of immune cells.

3.2. Stability of Equilibrium States

To study the stability of equilibrium states, we consider the Jacobian matrix of model (3)

$$J=\left(\begin{array}{cccc}{a}_{11}& 0& -N& -N\\ T& a_{22}& N-3.42\times {10}^{-3}L& -L\\ -T& -6.02\times {10}^{3}T& a_{33}& -kT\\ 0& 0& 0& -f\end{array}\right)$$

where

$$\begin{array}{cc}\hfill a_{11}& =p-3.6\times {10}^{-2}pN-T-\frac{s}{f},\hfill \\ \hfill a_{22}& =-1-3.42\times {10}^{-3}T-\frac{s}{f},\hfill \\ \hfill a_{33}& =25.7-1.05\times {10}^{-2}T-N-6.02\times {10}^{3}L-\frac{ks}{f}.\hfill \end{array}$$

3.2.1. Dead Equilibrium State

The Jacobian matrix of model (3) at $E_0(0,0,0,\frac{s}{f})$ is given by

$$J\left(E_0\right)=\left(\begin{array}{cccc}p-\frac{s}{f}& 0& 0& 0\\ 0& -1-\frac{s}{f}& 0& 0\\ 0& 0& 25.7-\frac{ks}{f}& 0\\ 0& 0& 0& -f\end{array}\right).$$

The eigenvalues of this matrix are

$$\begin{array}{c}\hfill {\lambda }_1=p-\frac{s}{f},{\lambda }_2=-1-\frac{s}{f}<0,{\lambda }_3=25.7-\frac{ks}{f},{\lambda }_4=-f<0.\end{array}$$

(15)

Using the Routh–Hurwitz theorem [38], the dead equilibrium state $E_0(0,0,0,\frac{s}{f})$ is locally asymptotically stable if ${\lambda }_1$ and ${\lambda }_3$ are negative, namely

$$\begin{array}{c}\hfill s>pf,s>\frac{25.7f}{k}.\end{array}$$

(16)

We have the following results about the stability of the dead equilibrium state $E_0(0,0,0,\frac{s}{f})$.

(1) $E_0$ is a locally asymptotically stable node when $s>max\{pf,\frac{25.7f}{k}\}$.

(2) $E_0$ is an unstable saddle node when $s=pf$ or $ks=25.7f$.

(3) $E_0$ is a saddle point when $s<pf$ or $s<\frac{25.7f}{k}$.

3.2.2. Tumor-Free Equilibrium State

When $p>\frac{s}{f}$, the tumor-free equilibrium state $E_1(\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p},0,0,\frac{s}{f})$ exists. The Jacobian matrix of model (3) at this state becomes

$$J\left(E_1\right)=\left(\begin{array}{cccc}-(p-\frac{s}{f})& 0& -\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p}& -\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p}\\ 0& -1-\frac{s}{f}& \frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p}& 0\\ 0& 0& 25.7-\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p}-\frac{ks}{f}& 0\\ 0& 0& 0& -f\end{array}\right).$$

The eigenvalues of the above matrix are

$$\begin{array}{c}\hfill {\lambda }_5=-(p-\frac{s}{f}),{\lambda }_6=-1-\frac{s}{f}<0,{\lambda }_7=25.7-\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p}-\frac{ks}{f},{\lambda }_8=-f<0.\end{array}$$

(17)

Since $p>\frac{s}{f}$, then ${\lambda }_5=-(p-\frac{s}{f})<0$. Therefore, the tumor-free equilibrium state $E_1$ is locally asymptotically stable if and only if ${\lambda }_7<0$. After calculation, we obtain that ${\lambda }_7$ is always negative when parameters values satisfy any one of the following conditions (18) and (19):

$$\begin{array}{ccc}\hfill 0.5374pf& <& s<pf,k>\frac{s-0.5374pf}{1.8\times {10}^{-2}ps};\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill 0& <& s<0.5374pf.\hfill \end{array}$$

(19)

Using the Routh–Hurwitz theorem [38], we can obtain the stability results of the tumor-free equilibrium state $E_1$.

(4) When model parameters satisfy one of the above inequalities (18) and (19), the tumor-free equilibrium state $E_1$ is locally asymptotically stable.

(5) When $25.7=\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p}+\frac{ks}{f}$, and $p>\frac{s}{f}$, $E_1$ is an unstable saddle node.

(6) When $25.7>\frac{p-\frac{s}{f}}{1.8\times {10}^{-2}p}+\frac{ks}{f}$, and $p>\frac{s}{f}$, $E_1$ is a saddle node.

Note that if the tumor-free equilibrium state of model (3) is unstable, the tumor cells cannot be eliminated by any amount of chemotherapy drugs.

3.2.3. Tumor-Present Equilibrium State

Note that when $25.7-\frac{ks}{f}>0$, there is the tumor-present equilibrium state $E_2(0,0,$$1.91\times {10}^2(25.7-\frac{ks}{f}),\frac{s}{f})$. The Jacobian matrix at this equilibrium state has the form

$$J\left(E_2\right)=\left(\begin{array}{cccc}\overline{a_{11}}& 0& 0& 0\\ 1.91\times {10}^2(25.7-\frac{ks}{f})& \overline{a_{22}}& 0& 0\\ -1.91\times {10}^2(25.7-\frac{ks}{f})& -1.15\times {10}^6(25.7-\frac{ks}{f})& -(25.7-\frac{ks}{f})& \overline{a_{33}}\\ 0& 0& 0& -f\end{array}\right).$$

where

$$\begin{array}{cc}\hfill \overline{a_{11}}& =p-1.91\times {10}^2(25.7-\frac{ks}{f})-\frac{s}{f},\hfill \\ \hfill \overline{a_{22}}& =-1-6.53\times {10}^{-1}(25.7-\frac{ks}{f})-\frac{s}{f},\hfill \\ \hfill \overline{a_{33}}& =-1.91\times {10}^{2}k(25.7-\frac{ks}{f}).\hfill \end{array}$$

The eigenvalues of matrix $J\left(E_2\right)$ are

$$\begin{array}{cc}{\lambda }_9\hfill & =p-1.91\times {10}^2(25.7-\frac{ks}{f})-\frac{s}{f},\hfill \end{array}$$

(20)

$$\begin{array}{cc}{\lambda }_{10}\hfill & =-1-6.53\times {10}^{-1}(25.7-\frac{ks}{f})-\frac{s}{f}<0,\hfill \end{array}$$

(21)

$$\begin{array}{cc}{\lambda }_{11}\hfill & =-(25.7-\frac{ks}{f}),\hfill \end{array}$$

(22)

$$\begin{array}{cc}{\lambda }_{12}\hfill & =-f<0.\hfill \end{array}$$

(23)

Since ${\lambda }_9=p-1.91\times {10}^2(25.7-\frac{ks}{f})-\frac{s}{f}$ and $25.7>\frac{ks}{f}$, we have that

(i) when $p-\frac{s}{f}<0$, ${\lambda }_9<0$;

(ii) when $p-\frac{s}{f}>0$ and ${\lambda }_9<0$, we can derive that $s>\frac{(p-4.91\times {10}^3)f}{1-1.91\times {10}^{2}k}$;

(iii) ${\lambda }_{11}<0$.

Thus, we can obtain the stability results of the tumor-present equilibrium state $E_2$ according to Routh–Hurwitz theorem [38].

(7) If and only if the parameters satisfy the inequalities $\left(24\right)$ or $\left(25\right)$, the tumor-present equilibrium state $E_2$ is locally asymptotically stable.

$$\begin{array}{ccc}\hfill pf& <& s<\frac{25.7f}{k};\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill \frac{(p-4.91\times {10}^3)f}{1-1.91\times {10}^{2}k}& <& s<min\{pf,\frac{25.7f}{k}\}.\hfill \end{array}$$

(25)

(8) When $p=1.91\times {10}^2(25.7-\frac{ks}{f})+\frac{s}{f}$, $25.7>\frac{ks}{f}$, the tumor-present equilibrium $E_2$ is an unstable saddle node.

(9) When $p>1.91\times {10}^2(25.7-\frac{ks}{f})+\frac{s}{f}$, $25.7>\frac{ks}{f}$, the tumor-present equilibrium $E_2$ is a saddle node.

3.2.4. Coexisting Equilibrium State

The key interest of this study is to find the existence and stability of the coexisting equilibrium state with a small number of tumor cells but a large number of immune effector cells. The Jacobian matrix of model (3) at $E_3(N_i^{\ast },L_i^{\ast },T_i^{\ast },\frac{s}{f})$, $(i=1,2,3)$ is given by

$$J\left(E_3\right)=\left(\begin{array}{cccc}{a}_{11}^{\ast }& 0& -N_i^{\ast }& -N_i^{\ast }\\ T_i^{\ast }& a_{22}^{\ast }& N_i^{\ast }-3.42\times {10}^{-3}{L}_i^{\ast }& -L_i^{\ast }\\ -T_i^{\ast }& -6.02\times {10}^3{T}_i^{\ast }& a_{33}^{\ast }& -k{T}_i^{\ast }\\ 0& 0& 0& -f\end{array}\right).$$

where

$$\begin{array}{cc}\hfill a_{11}^{\ast }& =p-3.6\times {10}^{-2}p{N}_i^{\ast }-T_i^{\ast }-\frac{s}{f},\hfill \\ \hfill a_{22}^{\ast }& =-1-3.42\times {10}^{-3}{T}_i^{\ast }-\frac{s}{f},\hfill \\ \hfill a_{33}^{\ast }& =25.7-1.05\times {10}^{-2}{T}_i^{\ast }-N_i^{\ast }-6.02\times {10}^3{L}_i^{\ast }-\frac{ks}{f}.\hfill \end{array}$$

Using the same analysis method, we can obtain the stability results of the coexisting equilibrium state $E_3(N_i^{\ast },L_i^{\ast },T_i^{\ast },\frac{s}{f})$, $(i=1,2,3)$ as follows.

(10) When $k>5.23\times {10}^{-3}$ and $f>f_1^{\ast }$, $E_3(N_1^{\ast },L_1^{\ast },T_1^{\ast },\frac{s}{f})$ is locally asymptotically stable.

(11) When $1.23\times {10}^4<p<2.92\times {10}^4$, $0<k<\frac{1}{1.8\times {10}^{-2}p}$, and $0<f<f_3^{\ast }$, $E_3(N_2^{\ast },L_2^{\ast },T_2^{\ast },\frac{s}{f})$ is a locally asymptotically stable focus.

(12) When 0 < p < 5.06 $\times {10}^3$, 0 < k < k${}_1^{\ast }$, and $f>f_2^{\ast }$, $E_3(N_3^{\ast },L_3^{\ast },T_3^{\ast },\frac{s}{f})$ is a saddle point.

4. Parameter Sensitivity Analysis

In this section, sensitivity analysis is used to assess the impact of small parameter changes on treatment outcomes by simulating the progression of cell numbers over 25 days. We fix all the model parameter values except one parameter. By varying the parameter value to a certain range, we check the difference of model outcome to that of the model outcome using the original model parameters. Figure 2 gives the sensitivity analysis results.

In model (3), we set $k=0.1$, $f=0.1$, $s=10$, and parameter $p\in (0,1)$ with an increment of 0.1. The initial condition is $N_0=1\times {10}^5$, $L_0=1\times {10}^2$, $T_0=1\times {10}^7$ and $u_0=10$. The simulation time period is $t=30$ days. Figure 2a shows that the number of tumor cells decreases with the increase of the NK cell growth rate p. However, there are no significant changes in the number of tumor cells, which indicates that NK cells have minor effects on the output of the model. This suggests that NK cells are not as important as the CTLs to kill tumor cells.

Similarly, we fix $p=20$, $f=0.1$ and $s=10$, but use $k\in (0,0.2)$ with an increment of 0.01. Figure 2b suggests that the tumor cell number decreases if the fatality rate of tumor cells k by the chemotherapy drug increases. Figure 2b also suggests that a small variation in parameter k may have a substantial influence on the tumor cell number.

In the third test, we fix $p=20$, $k=0.1$ and $f=0.1$, but s takes a value of 0 to 100. The parameter s represents the constant input of the chemotherapeutic drug. Figure 2c shows that the tumor cell number changes slightly when the drug input constant is changed. This indicates that the model is not particularly sensitive to parameter s, which suggests that we cannot achieve the therapeutic effect by increasing the drug input dose alone.

Finally, we choose $p=20$, $k=0.1$, and $s=10$, but f takes a value of 0 to 0.5, with an increment of 0.01. Figure 2d gives the number of tumor cells at $t=30$ days. This figure indicates that the faster the rate of chemotherapy drug decay, the larger the number of tumor cells in the host.

Figure 2 suggests that the model is sensitive to the tumor mortality rate k induced by chemotherapy drug, as well as the drug decay rate f. These results imply that in addition to the decay rate of the chemotherapy drug, any improvement in the drug efficiency may change the treatment results substantially. Since $k=\frac{k_T}{k_N}$, we can increase the chemotherapy drug-induced tumor mortality rate $k_T$, or decrease the side effects of therapy drug $k_N$. However, the tumor cell number is not very sensitive to the NK cells growth rate p and the chemotherapeutic drug input constant s. Thus, the cytolytic activity of the NK cells alone may not be a vital factor for the treatment outcomes. The activity of CTLs should also be considered.

5. Numerical Simulations

5.1. Numerical Simulations of the Equilibrium States

This subsection gives numerical simulations for the stability properties of the equilibrium states that have been proved in the previous section. The parameter values are given in Table 1 except p, k, s, and f. We use the same initial immune system strengths. Here, the parameters p, k, s, and f of each case are shown in

Table 4. The values of the fixed parameter p, k, s, and f for each simulation graph.

No.	p	k	s	f	Equilibrium Points
1	10	$0.4$	10	$0.1$	$E_0$
2	10	10	1	$0.1$	$E_1$
3	20	$0.1$	10	$0.2$	$E_2$
4	20	$0.01$	10	$0.6$	$E_3$

. The initial conditions of all cases are $(N_0,L_0,T_0,u_0)=(1\times {10}^5,1\times {10}^2,1\times {10}^7,10)$.

In Figure 3a, except for the drug concentration, the number of all cell populations tends to zero, that is, the equilibrium point $E_0$ is stable, indicating the death of the host. In Figure 3b, after a period of time, only NK cells and drugs are present, and the tumor cells are zero. Thus, the equilibrium point $E_1$ is stable. This is the ideal case for treatment. Figure 3c shows that only drugs and tumor cells remain, but the immune cells tend to zero. Thus, the equilibrium point $E_2$ is stable, which indicates that tumor cells can not be controlled by drug therapy, leading to the death of the host. Figure 3d depicts the coexistence of the immune cells and tumor cells.

5.2. Simulations Using Different Immune Strengths

To demonstrate the effect of immune system strengths, we hypothesize that the drug has the same input constant and the same decay rate. However, different hosts have different immune system strengths. Therefore, we fix the parameter values $k=2$, $f=1$, and $s=10$. The value of parameter p is in the range of $[10,20]$. Thus, $p=20$ is the strongest immune system strength, but $p=10$ is the weakest one. The initial concentrations of the drug and tumor cell number are the same, namely $(T_0,u_0)=(1\times {10}^7,10)$. However, the initial numbers of immune cells are different. So the initial condition in Figure 4a is $N_0=1\times {10}^5$, and $L_0=1\times {10}^2$. However, the initial condition is $N_0=1\times {10}^4$ and $L_0=1\times 10$ in Figure 4b.

Figure 4 shows that different hosts with different immune system strengths have different influences on tumor cells treated with the same drug. Figure 4a shows that, after a period of treatment, only NK cells and chemotherapy drugs remain in the host. Figure 4b suggests that, under the same chemotherapeutic drug, a large number of tumor cells are present in the host when the immune system is weaker. These results suggest that a host with a stronger immune system can eliminate a small number of tumor cells under a single chemotherapy. Therefore, treatments that increase the strength of the immune system should be developed, which is a different treatment scheme (i.e., immunotherapy).

6. Discussion

Through the analysis of system equations, we determine the existence of the equilibrium points and the conditions of corresponding stability and symmetry. For a specific parameter set, we find the model may have four equilibrium states. The first one is the dead equilibrium state, which means that both immune cells and tumor cells approach zero. The second one is the tumor-free equilibrium state, in which there are only NK cells but the tumor cell number is zero. The third case is the tumor-present equilibrium state, in which there are only tumor cells in the host. The last one is the coexisting equilibrium point, which suggests that both immune and tumor cells coexist with nonzero populations. The purpose of this study is to explore the model conditions that lead to the elimination of tumor cells, namely the zero-tumor burden or tumor-free equilibrium points. Therefore, the stability of the tumor-free equilibrium is very important. The instability of the tumor-free equilibrium state means that any successful treatment must be able to alter system parameters in order to stabilize this equilibrium. However, the stability of the tumor equilibrium with a large value means that the reduction of tumor burden through chemotherapy alone is not enough to eliminate tumor cells. Once chemotherapy is stopped, the system eventually returns to the state with a large tumor number. According to the above results, the tumor-free equilibrium point $E_1$ is stable when parameters p and s satisfy the conditions (18) or (19). Therefore, we can change the stability of the equilibrium with zero tumor number by modifying the model parameters with chemotherapy. These results suggest a potential approach to improve the effects of chemotherapy and to design new therapy schemes.

We also analyze the sensitivity of model parameters. It should be noted that the chemotherapy drug-induced tumor mortality rate k and the drug decay rate f have the greatest influence on the treatment outcome, as seen in Figure 2. However, the number of NK cells does not have much influence on tumor growth. In addition, the cytolytic activity of the NK cells should be considered together with the CTL activity to eliminate tumor cells. Numerical simulations in Figure 3 show the importance of the host immune system and explain the side effects of chemotherapy drugs. As can be seen from the simulations in Figure 4, treatments that increase the strength of the host’s immune system should be implemented. Based on the qualitative conclusions, it is suggested to search for adjuvant treatments that increase the chemotherapy drug-induced tumor mortality rate, decrease the drug decay rate, and/or reduce the side effects of therapeutic drugs.

The main contribution of this paper is the application of two immune cell populations and chemotherapy drugs to decrease the tumor burden. However, the chemotherapeutic drug considered in our model is a constant input, which is a limitation to this study. Nevertheless, we believe that our simple model can lay the foundation for the development of more sophisticated models for tumor treatment. The future work will improve the model in this work by taking into account drug diffusion, drug-delivery methods, and drug resistance. More sophisticated models will provide more realistic simulations of tumor populations and the effects of different treatment schemes.

7. Conclusions

In conclusion, this work has extended previous mathematical models that describe tumor–immune interactions. We introduce chemotherapy into the model and propose a model to describe tumor–immune interactions under chemotherapy. We analyze the equilibrium states of the proposed model and also investigate the stability and symmetrical properties of these steady states. Sensitivity analysis is conducted to find important model parameters that have more influence on the model dynamics than other parameters. The theoretical results provide predictions for designing new treatment schemes.