Stability analysis of a targeted chemotherapy-cancer model

1.   Introduction

Cancer, the abnormal growth of tumor cells that invade other body parts, is now the second leading cause of death worldwide after cardiovascular disease. So, managing cancer disease by developing new treatment strategies is a new research area for researchers. However, it is essential to understand the dynamics of tumor cells' growth and their complex interactions with the immune system to develop new treatment strategies. To do this, researchers heavily relied on mathematical models.

Many scientists have studied mathematical modelling of tumor evolution, interaction with different cells, and tumor proliferation by developing various models over the last few decades ^{[1,2,3,4,5,6,7,8]}. Allison et al. ^[9] developed an experiment-based, mathematical model to measure the effect of tumor growth rate, carrying capacity, and cytolytic activity of the effector cells on the progression of a tumor. Li et al. ^[10] explored the effects of angiogenic growth factors secreted by the tumor associated with the angiogenic process on tumor growth using a nonlinear time-delay model. They observed that the time delay on angiogenic growth is harmless for the model's local and global dynamical properties. In ^[11], the authors observed the dynamical properties such as stable steady states, unstable steady states, and limit cycles of a tumor micro-environment mathematical model consisting of tumor cells, cytotoxic T cells, and helper T cells by varying the infiltration rates of cytotoxic and helper T cells. Their results showed that the tumor always escaped when the infiltration rates of cytotoxic and helper T cells were both low. Two types of bistability can occur at intermediate values of the infiltration rates of cytotoxic and helper T cells. One is bi-stability between tumor escape and elimination; the other is between tumor escape and tumor coexistence with the immune system, the coexistence state representing either a finite equilibrium state or a time-dependent periodic solution. The tumor is permanently eliminated when the infiltration rates of cytotoxic and helper T cells are high. Ghosh and Banerjee ^[12] investigated the role of antibodies in cancer elimination by considering the clinical fact that antibodies can directly kill cancerous cells. Their quantitative analysis reveals that the antibodies' effectiveness plays a significant role in killing cancer cells directly. Shu et al. ^[13] explored the interactions between tumor cells, M1 and M2 macrophages and suggested some new results, such as the dual role of M1 and M2 immune macrophages on the tumor and the balance between the immune system and tumors in the host, which are essential from a biological perspective.

There are several common treatment modalities for cancer, such as chemotherapy, targeted chemotherapy, radiotherapy, immunotherapy, and surgery. Among these cancer treatments, targeted chemotherapy is a systematic therapy that fights and kills the cancer cells in the tumor site without any significant effect on the effector-normal cells so that the tumor cells can not migrate to other parts of the body ^[14]. Mathematical modeling of tumor growth and treatment has been approached by several researchers using a variety of models over the past few decades ^[15,16,17]. De Pillis and his collaborators investigated several mathematical models to measure the impact of chemotherapy ^[18,19,20], immunotherapy ^[17,21], chemo-immunotherapy ^[22], and monoclonal antibody (mAb) therapy ^[23] on tumor growth and other healthy tissue. In ^[19], the authors observed "Jeff's phenomenon" for a tumor growth competition model. They obtained new optimal treatment protocols that lower tumor growth and stabilize the normal cell population. A phase-space analysis of a mathematical model of tumor growth with immune response and chemotherapy has been performed, and it has been observed that optimally administered chemotherapy drugs could drive the system into a desirable basin of attraction ^[18]. In ^[24], the importance of CD8+ T cell activation in cancer therapy was addressed. Moreover, Pillis et al. developed a cancer treatment model ^[21] in which they observed that combining chemotherapy and immunotherapy can eliminate the entire tumor rather than the therapies applied alone. An optimally controlled chemo-drug administration was presented in ^[20], which discussed the role of a quadratic control, a linear control, and a state constraint. In ^[22], it is observed that CD8+T cells play an active role in chemo-immunotherapy to kill tumor cells. The effectiveness of mAbs treatment on colorectal cancer has been discussed in ^[23]. An optimal feedback scheme was proposed based on ^[24] that aims to tumor regression under a better health indicator profile and to improve treatment strategies in the case of mixed immunotherapy and chemotherapy of tumor ^[25]. Ansarizadeh et al. ^[26] examined the effectiveness of chemotherapy on tumor regression by developing a model of partial differential equations. They observed that the closer the point of injection of the chemotherapeutic drug is to the invasive fronts, the more regression in tumor cells occurs, and they also observed a clinically verified phenomenon called "Jeff's phenomenon." Arabameri et al. ^[27] designed a dendritic cell-based immunotherapy model by fitting it with experimental data, which predicted that the successive injection of DC-based vaccines might be very effective in suppressing tumor cells. A combination of stem-cell and chemo-drug therapy models for cancer treatment was investigated in ^[28], and the results suggested that the combination treatment may cure cancer and improve the patient's life. Recently, Liu and Liu ^[14] proposed a targeted chemotherapy cancer treatment mathematical model that suggests that the effectiveness of targeted chemotherapy in killing tumor cells is better than regular chemotherapy. The effect and efficacy of the targeted chemotherapeutic drug were investigated in ^[29], which shows that an adequate dosage of the targeted chemotherapeutic drug of low molecular weight is necessary for removing tumor cells from the infected tumor system.

In the present study, we will examine the effectiveness of targeted chemotherapy on tumor and normal cells by modifying the model of de Pillis et al. ^[18], which describes the interaction between tumor cells and healthy tissue or normal cells. In the very next section, we formulate our model. In Section 3, we verified the positive invariance and boundedness of the models' solutions. A steady-state analysis has been performed in Section 4. We have checked the global stability of the tumor-free state in Section 5. A numerical simulation and a concluding remark have been carried out in Sections 6 and 7.

2.   Model formulation

We consider the model developed by de Pillis et al. ^[18], in which a set of ordinary differential equations is used to characterize the dynamic interplay between effector cells, tumor cells, and healthy cells. We employ a novel adaptation in the model proposed by de Pillis et al. ^[18] that a targeted chemotherapy approach slows tumor growth. Traditional chemotherapy drugs kill all cell types at varying rates, which can cause severe side effects like hair loss, fatigue, anaemia, etc. On the other hand, targeted chemotherapy mainly targets tumor cells so that it may have fewer unintended consequences. We describe this attachment by the term $–kTC$ in the last equation of the system (2.1), where $k$ denotes the rate of attachments of targeted chemo-drugs with tumor cells. Also, to measure the effectiveness of targeted chemotherapy on immune and normal cells, we introduce a parameter $\eta$ ^[14]. Thus, our modification takes the following model:

and the initial conditions are: $I(0) = I_{0} > 0$, $T(0) = T_{0} > 0$, $N(0) = N_{0} > 0$ and $C(0) = C_{0} > 0. \; I(t)$, $T(t)$, and $N(t)$ represent the densities of effector cells, tumor cells, and normal cells, respectively, at time $t$ and $C(t)$ represents the amount of targeted chemo drug administered at time $t$. In terms of the model's novelty, we normalized the number of normal cells by taking their carrying capacity equal to one.

3.   Positive invariance and boundedness

In this section, we will investigate whether the model (2.1) solutions are biologically feasible or not for the considered values of all parameters. To do this, we must show that the solutions of the model are positive and bounded. Using standard comparison theory ^[30], we get

Integration of the above leads to

Again,

Proceeding as above, we have

Similarly, we have

and

Using the considered initial values, we assume that $I(t)\geq 0, T(t) \geq 0, N(t) \geq 0$ and $C(t)\geq 0$ for all $t > 0$. Consequently, the corresponding domain region for the system (2.1) is

The domain region $\Delta$ is positively invariant, which verifies that the model system (2.1) is biologically feasible.

4.   Equilibrium and stability

In this section, we will study the stability of the system and the system's stability around the equilibrium. To do this, we compute $\frac{dI}{dt} = 0$, $\frac{dT}{dt} = 0$, $\frac{dN}{dt} = 0$, and $\frac{dC}{dt} = 0$ and get the following equilibrium:

● Tumor-free equilibrium: $E_{1} (I_{1}, T_{1}, N_{1}, C_{1})$ where,

$I_{1} = \frac{r_{1}-c_{3}r_{2}+\{c_{3}a_{3}(1-\eta)-a_{2}\}\frac{u}{d_{2}}}{c_{2}}$, $T_{1} = 0$, $N_{1} = \frac{r_{2}-a_{3}(1-\eta)\frac{u}{d_{2}}}{r_{2}}$,

$C_{1} = \frac{u}{d_{2}}.$

The tumor-free equilibrium $E_{1}$ exists if $r_{1}+ \{c_{3}a_{3}(1-\eta)-a_{2}\}\frac{u}{d_{2}} > c_{3}r_{2}$ and $r_{2} > a_{3}(1-\eta)\frac{u}{d_{2}}$.

● Co-axial equilibrium: $E^{*}(I^{*}, T^{*}, N^{*}, C^{*})$ where,

$I^{*} = \frac{r_{1}(1-b_{1}T^{*}-c_{3}r_{2}+c_{3}c_{4}T^{*})+\{c_{3}a_{3}(1-\eta)-a_{2}\}\frac{u}{d_{2}+kT^{*}}}{c_{2}}$,

$N^{*} = \frac{r_{2}-c_{4}T^{*}-a_{3}(1-\eta)\frac{u}{d_{2}+kT^{*}}}{r_{2}}$, $C^{*} = \frac{u}{d_{2}+kT^{*}}$, and from the following quadratic equation, we will calculate the value of $T^{*}$.

After substituting the values of $I^{*}$ and $C^{*}$, we get

where,

The co-axial equilibrium $E^{*}$ exists if the roots of the equation (4.1) is positive i.e., $T^{*} > 0$ and following inequalities must holds:

As $N = 0$ implies that the patients will not be alive, we do not consider those cases where $N = 0$. To investigate the linear stability of the system around the two above stability, we must compute the Jacobian of the system, and the Jacobian is

where

$m_{11} = \frac{\rho T}{\sigma + T}-c_{1}T-d_{1}-a_{1}(1-\eta)C$,

$m_{12} = \frac{\sigma\rho I}{(\sigma + T)^{2}}-c_{1}I$,

$m_{14} = -a_{1}(1-\eta)I$,

$m_{21} = -c_{2}T$, $m_{22} = r_{1}-2r_{1}b_{1}T-c_{2}I-c_{3}N-a_{2}C$,

$m_{23} = -c_{3}T$,

$m_{24} = -a_{2}T$,

$m_{32} = -c_{4}N$,

$m_{33} = r_{2}-2r_{2}N-c_{4}T-a_{3}(1-\eta)C$,

$m_{34} = -a_{3}(1-\eta)N$,

$m_{42} = -kC$,

$m_{44} = -d_{2}-kT$.

● The eigenvalues of the Jacobian matrix (4.2) corresponding to the steady-state $E_{1}$ are: $\lambda_{1}^{1} = -d_{1}-a_{1}(1-\eta)C_{1}$, $\lambda_{2}^{1} = r_{1}-c_{2}I_{1}-c_{3}N_{1}-a_{2}C_{1}$, $\lambda_{3}^{1} = r_{2}-2r_{2}N_{1}-a_{3}(1-\eta)C_{1}$, and $\lambda_{4}^{1} = -d_{2}$. Clearly, $\lambda_{1}^{1} < 0$ and $\lambda_{4}^{1} < 0$.

Therefore, $E_{1}$ will stable if $\lambda_{2}^{1} < 0 \; \implies \; r_{1} < c_{2}I_{1}+c_{3}N_{1}+a_{2}C_{1}$ and $\lambda_{3}^{1} < 0 \; \implies \; r_{2} < 2r_{2}N_{1}+a_{3}(1-\eta)C_{1}$; otherwise $E_{1}$ becomes unstable.

● The Jacobian corresponding to the stability $E^{*}$ is

where

$m_{11}^{*} = \frac{\rho T^{*}}{\sigma + T^{*}}-c_{1}T^{*}-d_{1}-a_{1}(1-\eta)C^{*}$, $m_{12}^{*} = \frac{\sigma\rho I^{*}}{(\sigma + T^{*})^{2}}-c_{1}I^{*}$, $m_{14}^{*} = -a_{1}(1-\eta)I^{*}$, $m_{21}^{*} = -c_{2}T^{*}$, $m_{22}^{*} = r_{1}-2r_{1}b_{1}T^{*}-c_{2}I^{*}-c_{3}N^{*}-a_{2}C^{*}$, $m_{23}^{*} = -c_{3}T^{*}$, $m_{24}^{*} = -a_{2}T^{*}$, $m_{32}^{*} = -c_{4}N^{*}$, $m_{33}^{*} = r_{2}-2r_{2}N^{*}-c_{4}T^{*}-a_{3}(1-\eta)C^{*}$, $m_{34}^{*} = -a_{3}(1-\eta)N^{*}$, $m_{42}^{*} = -kC^{*}$, and $m_{44}^{*} = -d_{2}-kT^{*}$.

At $E^{*}$, the eigen values of the corresponding Jacobian matrix (4.3) are the roots of the following equation

where

According to Routh-Hurwitz stability rule, the equilibrium point $E^{*}$ will be asymptotically stable if

5.   Global stability at tumor-free steady state $E_{1}$

In this section, we will analyze the global stability around $E_{1}$ in order to investigate the behavior of the system (2.1) far away from the equilibrium point $E_{1}$. Let us define a Lyapunov function for the model (2.1) at $E_{1}$ as:

Differentiating $V(t)$ with respect to $t$, we get

where,

$P = [I-I_{1}, T, N-N_{1}, C-C_{1}]$ and

$Q = [0, -r_{1} +c_{2}I_{1}+c_{3}N_{1}+a_{2}C_{1}, 0, 0]$.

By noting the second component of the vector $Q$, we must have,

so that $Q^{t} P \geq 0$.

Furthermore, by considering the values of parameters from Table 1, if $I = \frac{s}{d_{1}}, T = \frac{1}{b_{1}}$ and $C = \frac{u}{d_{2}}$, then all minors are positive in the matrix $R$ (all eigenvalues of $R$ are also positive), and so $P^{t}RP > 0$. Hence, it is clear that $\frac{dV}{dt} < 0$.

Therefore, the tumor-free equilibrium $E_{1}$ satisfies local stability conditions, making the point globally stable. In biological terms, it means that targeted chemotherapy will kill the tumor cells if

must hold.

6.   Numerical simulation

Analytical studies can only be completed with numerical verification of the derived results. In this section, we verified our analytical results of the considered system (2.1) graphically using MATHEMATICA, which is very important from a practical point of view. All the simulations have been carried out using the parameter values of Table 1 ^[14,18]. We take the units of the parameter values to be arbitrary.

To investigate the effects of drugs on effector cells, tumor cells, and normal cells, we considered three cases: (a) $u = 0.019$, (b) $u = 0.020$ and (c) $u = 0.021$. For $u = 0.019$, there exist only co-axial equilibrium $E^{*} = (0.198, 0.0101, 0.725, 0.379)$. The eigenvalues correspond to coaxial equilibrium $E^{*}$ are $(-0.254, -0.248, -0.0504, -0.0096)$, implies that $E^{*}$ is stable. Further, for $u = 0.020$, only the co-axial equilibrium $E^{*} = (0.184, 0.0032, 0.715, 0.3997)$ exist and $(-0.270326, -0.250733, -0.0501065, -0.00282614)$ are the eigenvalues correspond to it. Hence, $E^{*}$ is stable for this case also. Thus, for the cases: (a) $u = 0.019$ (b)$u = 0.020$, only the co-axial equilibria exist, and tumor cells are present at this equilibrium; the drug dose could not be able to eradicate the tumor cells from the body, and hence there is a chance of the rebirth of abnormal cells. However, for the case of $u = 0.021$, there exist only tumor-free equilibrium $E_{1} = (0.177, 0, 0.703, 0.42)$. Corresponding to the equilibrium $E_{1}$, the eigenvalues are $-0.283, -0.246, -0.05, -0.0037$; which indicates that $E_{1}$ is a stable point. Therefore, for $u = 0.021$, no tumor cells are present in the body, and the patient recovers from the disease.

The above scenarios can be justified in Figure 1. We observe that the density of fast proliferation tumor cells gets suppressed quickly as the amount of drugs increases from $u = 0.019$, $u = 0.020$ to $u = 0.021$ (see Figure 1) and gets stable at zero for $u = 0.021$, suggesting targeted chemotherapy's success. So, the prescribed drugs quickly affected the tumor cells, which is clinically reliable.

In Figure 2, we observe that the density of effector cells decreases slowly (compared to tumor cells) while the amount of applied drug doses increases and gets stable at a required level.

In Figure 3, the density of normal cells is reduced more slowly (compared to tumor cells) while the amount of drugs increases. Also, the normal cells become stable at the desired level.

The time series diagrams of tumor cells for different tumor growth rates, $r_{1}$ with drug dose $u = 0.021$, have been presented in Figure 4.

After investigating the effect of drugs on the cell population, we are now interested in examining the interaction of tumor growth and drug administration. To examine this, we fixed the drug parameter $u = 0.021$ (as at this level, the tumor cells can be eradicated) and varied the tumor growth rate $r_{1}$. Also, we considered tumor growth rate as $r_{1} = 0.4, 0.42, 0.44$ along with the parameter Table 1. For $r_{1} = 0.4$, there exist only tumor-free equilibrium $E_{1} = (0.177, 0, 0.703, 0.42)$ corresponding to the eigenvalues $(-0.283, -0.246, -0.05, -0.0037)$; indicating $E_{1}$ is stable. Further, for $r_{1} = 0.42$ and $r_{1} = 0.44$, only the co-axial equilibrium $E^{*} = (0.205, 0.0175, 0.691, 0.419)$ and $E^{*} = (0.240, 0.0355, 0.680, 0.417)$ exist respectively. The respective eigenvalues are $(-0.239+0.001i, -0.239-0.001i, -0.0508, -0.0178)$ and $(-0.241, -0.181, -0.0553, -0.0422)$; hence, $E^{*}$ is stable for both cases. Biologically, the above results can be explained that the prescribed drugs eradicate the density of tumor cells if the tumor growth rate $r_{1}$ is low. From Figure 4, it can be concluded that the increase in the tumor growth rate decreases the tumor-reducing capability of the system.

Next, we fixed the value of treatment parameter $u = 0.0205$ and tumor growth rate, $r_{1} = 0.4$, to examine the optimal effects of drugs on the cell population. In this case, the system (2.1) exhibits one biologically valid equilibrium point and this $E_{1} = (0.178, 0, 0.710, 0.41)$, which is tumor-free. The eigenvalues corresponds to $E_{1}$ are $-0.28118, -0.248475, -0.05, -0.0004$, which indicates that $E_{1}$ is asymptotically stable node. From Figure 5, it can be observed that the tumor cells die off over a long period, where as the normal and immune cells are stable at their required level. Also, it is noticed that the system could reduce the density of tumor cells if tumor growth rate, $r_{1} = 0.4$ and drug dose, $u = 0.0205$.

In Figure 6, it is seen that the trajectory converges to the tumor-free steady state $E_{1}$ with the basin of attraction in the treatment case, indicating that it is a globally stable point for the system. The steady-state $E_{1}$ is a stable node, implying that the cell population incorporated with the treatment can suppress the cancer growth to zero with time increase. Biologically, this indicates that the body is recovering from the tumor regardless of the initial condition, which includes tumor growth.

Overall, we observe that if the size of the tumor is small, i.e., the growth rate of tumor $r_{1}$ is tiny, then it is quite possible to eradicate the tumor from the body using a smaller amount of drugs with less harm to the other healthy cells. If not, we require a high drug dose that can increase the side effects of the drugs.

7.   Conclusions

We have investigated a modified ODE mathematical model for tumor growth, considering the effector-tumor-normal cells' interaction under targeted chemotherapy. We established the basic characteristics, such as positivity and boundedness of the solutions of the model. To explore the dynamic behavior of the model, we performed a stability analysis of the considered system. We found that the tumor-free steady state is globally stable under the conditions: $c_{2}I_{1}+c_{3}N_{1}+a_{2}C_{1} > r_{1}, I = \frac{s}{d_{1}}, T = \frac{1}{b_{1}}, C = \frac{u}{d_{2}}$, provided it is locally stable; which suggests that the prescribed treatment can eradicate tumor cells from the body for a threshold value of tumor growth rate. Numerically, it is also observed that the prescribed treatment can eradicate tumor cells from the body without much effect on other healthy cells if the tumor size is small. However, one limitation in our model is that if the tumor size is large, it requires a high amount of drugs and a long period of time, which can harm the patient's body. So, in this regard, we must need an optimum period and drug dose for which the tumor is eradicated ^[18,19]; and it will be carried out in our future study.

Conflict of interest

The authors declare that they have no conflict of interest.