Dynamical Behavior of a Cancer Growth Model with Chemotherapy and Boosting of the Immune System

Abstract

In this study, we set up and analyze a cancer growth model that integrates a chemotherapy drug with the impact of vitamins in boosting and strengthening the immune system. The aim of this study is to determine the minimal amount of treatment required to eliminate cancer, which will help to reduce harm to patients. It is assumed that vitamins come from organic foods and beverages. The chemotherapy drug is added to delay and eliminate tumor cell growth and division. To that end, we suggest the tumor-immune model, composed of the interaction of tumor and immune cells, which is composed of two ordinary differential equations. The model’s fundamental mathematical properties, such as positivity, boundedness, and equilibrium existence, are examined. The equilibrium points’ asymptotic stability is analyzed using linear stability. Then, global stability and persistence are investigated using the Lyapunov strategy. The occurrence of bifurcations of the model, such as of trans-critical or Hopf type, is also explored. Numerical simulations are used to verify the theoretical analysis. The Runge–Kutta method of fourth order is used in the simulation of the model. The analytical study and simulation findings show that the immune system is boosted by regular vitamin consumption, inhibiting the growth of tumor cells. Further, the chemotherapy drug contributes to the control of tumor cell progression. Vitamin intake and chemotherapy are treated both individually and in combination, and in all situations, the minimal level required to eliminate the cancer is determined.

1. Introduction

Mathematical modeling can be used in our life to calculate a wide range of issues such as ecology, biology, and epidemiology [1]. Mathematical models can predict epidemics by employing fundamental assumptions of statistics and mathematics to determine infectious disease parameters and quantify the impact of interventions such as mass vaccination campaigns [2]. In particular, in the USA, the National Cancer Institute estimated that cancer survivors will be about 24.4% of all cancer patients in 2023 [3].

Cancer is a complex disease involving various ways in which aberrant cells can interact with the environment around them [4]. Cancer develops when the body’s immune cells cannot stop the proliferation of aberrant cells. When this happens, the body cannot control the formation of abnormal cells [5]. Treating cancer with surgery, radiation, chemotherapy, immunotherapy, and hormone therapy is prominently used. These treatments can be administered singly or in combination with two or more of the modalities mentioned above [6,7,8,9,10,11]. Chemotherapy is a form of treatment that makes extensive use of medications that have been chemically formulated [12]. Recently, there has been a sharp increase in interest in the mathematical modeling of tumor-immune dynamics, and many modeling methodologies have been employed to characterize these phenomena. Many tumor-immune models have been created using various types of equations, such as delay differential equations, ordinary differential equations, and partial differential equations, to formulate cancer models and see how tumor development affects the dynamics of other cells [13]. These models have led to the development of novel cancer medicines [1,12,14]. For example, Alqudah used an autonomous system to formulate a model of chemotherapy stem cells to treat cancer. She concludes that the treatments’ effects might assist in increasing the pace of effector cells to affect the immune system, resulting in the decay of tumor cells in cancer patients [7].

Numerous epidemiological studies have demonstrated the link between elevated cancer death rates and environmental changes in diet, pollution, lifestyle, and other variables. [15,16,17,18]. The relationship between a healthy diet, vitamin group, and the strengthening of the immune system has been under the microscope recently. It has been demonstrated that vitamins play an essential part in the regulation of the activity of the immune system, which is necessary for the protection of tissues from damage [17,19,20,21]. Ku-Carrillo developed an obesity–tumor model that has highlighted the role played by obesity in the tumor’s resistance to chemotherapy [22]. Rambely and his collaborator created a healthy immune system model that dynamically depicts how the immune system inhibits the progression of aberrant cells into tumors to contrast the analytical results of an unhealthy immune system model [23]. The numerical simulation of the model in [24] reveals that the immune system can be strengthened when an individual consistently consumes vitamins, at a daily rate of 16%, as a result of their effect on the body’s response to the formation of aberrant cells in the tissue.

The rate at which tumor cells multiply and disseminate varies significantly from species to species. Consequently, the rate at which cancer cells are suppressed or eliminated by the immune cell response may also vary. Many different mathematical models have been developed describing the immune system’s reaction to the suppuration of tumor cells. Most of these models have used the linear type of functional response. Some unhealthy behaviors that people engage in regularly, such as eating unhealthy food and smoking, can contribute to a reduction in the immune system’s performance. Because of this, we cannot utilize the same functional response to describe how well the immune system deals with tumors [25]. In particular, Alharbi and her collaborator developed two mathematical models: the first model describes how tumor-immune interactions are affected when poor dietary habits compromise the immune system. The second model includes the beneficial effect of vitamin consumption on the immune system. They conclude that a patient’s immune system might improve by taking vitamins consistently, at a daily rate of 55% [26].

In this paper, a cancer–immune–chemotherapy–vitamins model (CICV) governed by systems of ordinary differential equations is suggested based on models in [26]. We have modified Alharbi and Rambely’s model by changing the linear type of functional response to the Holling type II in order to describe the elimination of tumor cells by the immune system to account for the fact that the immune system may be weak. In addition, we used the Holling type III rather than the Holling type II functional response when defining the immune cells’ capacity to eradicate tumors. In addition, the impact of chemotherapy treatment of cancer and the regular intake of vitamins to support the immune cells is considered.

The outline for the rest of the paper is arranged as follows. The CICV model is introduced in Section 2. The existence of equilibria is presented in Section 3. In Section 4 and Section 5, the local stability of equilibria and local bifurcations are investigated. In Section 6, the numerical simulations of the CICV model are performed. The results of this paper are discussed in Section 7. Finally, a conclusion and outlook are given in Section 8.

2. Description of the Model

Consider a system of differential equations (CICV) consisting of tumor cells $C(t)$ and immune cells $I(t),$ which are presented as

$$\left.\begin{array}{c}\frac{dC}{dt}=r_{1}C\left(1-r_{2}C\right)-\frac{{\alpha }_{1}CI}{k_1+C}-{\beta }_{1}C=f_1\left(C,I\right),\\ \frac{dI}{dt}=\delta -{\alpha }_{2}CI+\frac{{\alpha }_3{C}^{2}I}{C^2+k_2}-{\beta }_{2}I+\gamma =f_2\left(C,I\right),\end{array}\right\}$$

(1)

with the initial conditions $C\left(0\right)=C_0\ge 0$, $I\left(0\right)=I_0\ge 0$. In the first equation of system (1), the term $r_{1}C\left(1-r_{2}C\right)$ represents the logistic growth of tumor cells with the growth rate $r_1$ and cancer cell capacity $r_2$. The second term $\frac{{\alpha }_{1}CI}{k_1+C}$ of Michaelis–Menten form models the killing of tumor cells by immune cells. The final term ${\beta }_{1}C$ represents the killing of tumor cells due to chemotherapy. In the second equation, $\delta$ is the constant source at which immune cells are produced. ${\alpha }_{2}CI$ models the suppression of the activity of the immune cells due to the action and rapid division of the tumor cells. The term $\frac{{\alpha }_3{C}^{2}I}{C^2+k_2}$ of Holling type III refers to the presence of tumor cells that incite the immune system’s reaction. ${\beta }_{2}I$ signifies the decay of immune cells due to natural death and chemotherapy. The last term $\gamma$ represents the increase of immune cells due to taking supplement vitamins regularly at a constant rate. All parameters for the tumor–immune–chemotherapy drug–vitamins model (CICV) are assumed to be positive and are clearly described in

Table 1. Explanation of CICV system’s parameters.

Parameter	Explanation	Values	Source
$r_1$	Tumor cell growth rate	$0.4426$ (day−1)	[26]
$r_2$	Tumor cell capacity	0.4 (cell−1)	[26]
${\alpha }_1$	Tumor clearance rate	$0.1469$ (cell−1)	[26]
$k_1$	The half-life of tumor cells	$0.76$ (day−1)	estimated
${\beta }_1$	The decay rate of the tumor cells due to taking chemotherapy	$0.2$ (cell−1 day−1)	estimated
$\delta$	The number of immune cells that are regularly produced in the body	$0.7$ (cell−1)	[26]
${\beta }_2$	The natural death rate of immune cells and decay rate of the immune cells killed by chemotherapy	$0.57$ (cell−1 day−1)	[26]
${\alpha }_2$	Rate of suppression of immune cells due to the activity and rapid division of the tumor cells	$0.3634$ (cell−1)	[26]
${\alpha }_3$	Immune cells’ rate of reaction and capacity to obliterate tumor cells	$0.7829$ (cell−1)	[26]
$k_2$	The half-life of immune cells	$0.8620$ (day−1)	[26]
$\gamma$	The rate of effectiveness of vitamins on immune cells	$0.5463$ (dose)	[26]

.

Figure 1 exemplifies the schematic sketch of the CICV model under examination.

The right-hand side of the CICV system is entirely continuous and differentiable on $R_{+}^2=\left\{\left(C,I\right):C\ge 0,I\ge 0\right\}$ and hence locally Lipschitzian [27]; therefore, the solution $(C\left(t\right),I\left(t\right))$ of the CICV model with initial conditions $C\left(0\right)=C_0\ge 0$, $I\left(0\right)=I_0\ge 0$ exists, and it is unique.

Theorem 1. 

All the solutions$(C\left(t\right),I\left(t\right))$of system (1), which start in$R_{+}^2$, will remain in$R_{+}^2$.

Proof. 

By integrating the first function of the CICV model for $C(t)$ and with a positive initial condition ${(C}_0,I_0)$, we obtain

$$C\left(t\right)=C_0\mathit{exp}\left\{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds\right\}\ge 0$$

Then,

$$dI=\left(\frac{{\alpha }_3{C}^{2}I}{C^2+k_2}-{\alpha }_{2}CI-{\beta }_{2}I+\delta +\gamma \right)dt$$

$$I=\left(\frac{{\alpha }_3{\left[C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\right]}^{2}I}{{\left[C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\right]}^2+k_2}-{C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\alpha }_{2}I-{\beta }_{2}I+\delta +\gamma \right)dt$$

$$dI=\left(\frac{{\alpha }_3{\left[C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\right]}^{2}I}{{\left[C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\right]}^2+k_2}-{C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\alpha }_{2}I-{\beta }_{2}I+\delta +\gamma \right)dt$$

Consequently, after dropping the non-negative quantities, this yields

$$dI\ge \left(\frac{{\alpha }_3{\left[C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\right]}^{2}I}{{\left[C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\right]}^2+k_2}-{C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\alpha }_{2}I-{\beta }_{2}I\right)dt$$

Now, by integrating the above equation for $I(t)$, we obtain

$$I\left(t\right)\ge I_0\mathit{exp}\left\{{\int }_0^t\left[\frac{{\alpha }_3{\left[C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\right]}^2}{{\left[C_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}\right]}^2+k_2}-{\alpha }_2{C}_0{e}^{{\int }_0^t\left[r_1-{\beta }_1-r_1{r}_{2}C(s)-\frac{{\alpha }_{1}I(s)}{k_1+C(s)}\right]ds}-{\beta }_2\right]ds\right\}$$

Thus, from the definition of the exponential function, any solution $C\left(t\right),I\left(t\right))$ that starts inside of $R_{+}^2$ with positive initial conditions ${(C}_0,I_0)$ will remain in $R_{+}^2$. □

Theorem 2. 

Assume that the following conditions hold.

$$\left.\begin{array}{c}{r}_1>{\beta }_1\\ {\beta }_2>\frac{{\alpha }_3{\ell }_1^2}{{\ell }_1^2+k_2}\end{array}\right\}.$$

(2)

Then, all solutions$(C\left(t\right),I\left(t\right))$of system (1) with positive initial values${(C}_0,I_0)$which start in

$$\zeta =\left\{\left(C,I\right){\in R}_{+}^2,C\le {\ell }_1,I\le \frac{\delta +\gamma }{{\ell }_2}\right\}$$

where

$${\ell }_1=\frac{\left(r_1-{\beta }_1\right)}{r_1{r}_2}$$

and

$${\ell }_2={\beta }_2-\frac{{\alpha }_3{\ell }_1^2}{{\ell }_1^2+k_2}$$

are uniformly bounded.

Proof. 

From the first equation of the CICV model, we obtain

$$\frac{dC}{dt}\le \left(r_1-{\beta }_1\right)C-r_1{r}_2{C}^2$$

Using Bernoulli’s method, we obtain

$$\left(t\right)\le \frac{\left(r_1-{\beta }_1\right)}{r_1{r}_2\left[1-e^{-\left(r_1-{\beta }_1\right)t}\right]+\left(r_1-{\beta }_1\right)(1/C\left(0\right))e^{-\left(r_1-{\beta }_1\right)t}}$$

As $t\to \infty$, the following is obtained

$$\underset{t\to \infty }{lim}\sup C\left(t\right)\le \frac{\left(r_1-{\beta }_1\right)}{r_1{r}_2}={\ell }_1$$

Using the above bound for the tumor cell, the following is obtained by the procedure of separation of variables.

$$\underset{t\to \infty }{lim}\sup I\left(t\right)\le \frac{\delta +\gamma }{{\ell }_2}$$

Consequently, $C\left(t\right)$ and $I\left(t\right)$ will remain bounded. □

3. Existence of Equilibria

To find the equilibria of system (1), we set

$$\frac{dC}{dt}=\frac{dI}{dt}=0$$

This system has two non-negative solutions, i.e., steady states, namely,

• The cancer-free state $S_1=(0,\widehat{I})$, where

  $$\widehat{I}=\frac{\delta +\gamma }{{\beta }_2}$$
• The interaction state $S_2=(C^{\ast },I^{\ast })$, where

  $$I^{\ast }=\frac{\left(r_1-{\beta }_1-r_1{r}_2{C}^{\ast }\right)\left(k_1+C^{\ast }\right)}{{\alpha }_1}$$

$C^{\ast }$ is the positive root of the following fifth-order polynomial

$$P_1{C}^5+P_2{C}^4+P_3{C}^3+P_4{C}^2+P_{5}C+P_6=0$$

(3)

where

$P_1=r_1{r}_2{\alpha }_2>0$

$P_2={\alpha }_2\left(r_1{r}_2{k}_1+{\beta }_1-r_1\right)+r_1{r}_2\left({\beta }_2-{\alpha }_3\right)$

$P_3=\left({\beta }_2-{\alpha }_3\right)\left(r_1{r}_2{k}_1+{\beta }_1-r_1\right)+k_1{\alpha }_2\left({\beta }_1-r_1\right)+r_1{r}_2{\alpha }_2{k}_2$

$P_4=k_1\left({\beta }_1-r_1\right)\left({\beta }_2-{\alpha }_3\right)+{\alpha }_2{k}_2\left(r_1{r}_2{k}_1+{\beta }_1-r_1\right)+{\alpha }_1(\delta +\gamma ))+r_1{r}_2{\beta }_2{k}_2$

$P_5=k_1{k}_2{\alpha }_2\left({\beta }_1-r_1\right)+{\beta }_2{k}_2\left(r_1{r}_2{k}_1+{\beta }_1-r_1\right)$

$P_6=\left({\beta }_1-r_1\right)k_1{k}_2{\beta }_2+{\alpha }_1{k}_2\left(\delta +\gamma \right)$

The condition

$${\beta }_1<r_1(1-r_2{C}^{\ast })$$

guarantees that $I^{\ast }>0$. The above condition implies that the immune system breaks down if the tumor cells’ decay rate must be less than the tumor growth rate by some amount.

Now, by applying Descartes’ rule of signs [27], Equation (3) has a unique positive root if one of the following conditions is met

$$P_2>0,P_3>0,P_4>0,P_5>0,P_6<0$$

$$P_2>0,P_3>0,P_4>0,P_5<0,P_6<0$$

$$P_2>0,P_3>0,P_4<0,P_5<0,P_6<0$$

$$P_2>0,P_3<0,P_4<0,P_5<0,P_6<0$$

$$P_2<0,P_3<0,P_4<0,P_5<0,P_6<0$$

The fact that the CICV model has an interaction steady-state indicates a deficiency in the immune system.

4. Stability Analysis

We compute the Jacobian matrix in order to obtain the local stability of the equilibria above

$$J\left(C,I\right)=\left[\begin{array}{cc}{r}_1-2r_1{r}_{2}C-\frac{k_1{\alpha }_{1}I}{{\left(k_1+C\right)}^2}-{\beta }_1& \frac{-{\alpha }_{1}C}{k_1+C}\\ -{\alpha }_{2}I+\frac{2k_2{\alpha }_{3}CI}{{\left(k_2+C^2\right)}^2}& {-\alpha }_{2}C+\frac{{\alpha }_3{C}^2}{k_2+C^2}-{\beta }_2\end{array}\right]$$

(4)

After computing the Jacobian matrix, the local analyzing behavior of the equilibrium points of the CICV model is described in the following theorem.

Theorem 3. 

Assume that

$$r_1<\frac{{\alpha }_1\left(\delta +\gamma \right)}{k_1{\beta }_2}+{\beta }_1$$

(5)

Then,$S_1$is locally asymptotically stable.

Proof. 

The Jacobian matrix at $S_1$ is computed, and it is given as

$$J\left(S_1\right)=\left[\begin{array}{cc}{r}_1-\frac{{\alpha }_1\left(\delta +\gamma \right)}{k_1{\beta }_2}-{\beta }_1& 0\\ \frac{{-\alpha }_2\left(\delta +\gamma \right)}{{\beta }_2}& -{\beta }_2\end{array}\right]$$

Then, the eigenvalues of $J\left(S_1\right)$ are

$${\lambda }_{11}=r_1-\frac{{\alpha }_1\left(\delta +\gamma \right)}{k_1{\beta }_2}-{\beta }_1$$

and

$${\lambda }_{12}=-{\beta }_2<0$$

Clearly

$${\lambda }_{11}<0$$

and hence, $S_1$ is locally asymptotically stable if condition (5) is satisfied.

It can be realized from condition (5) that the immune system is functioning properly. Qualitatively, condition (5) means that the rate of suppression of tumor cells by immune cells plus the decay rate of tumor cells will be more than the growth rate of tumor cells $r_1$. □

Theorem 4. 

Suppose that

$$\left.\begin{array}{c}{r}_1<2r_1{r}_2{C}^{\ast }+\frac{k_1{\alpha }_1{I}^{\ast }}{{\left(k_1+C^{\ast }\right)}^2}+{\beta }_1+\frac{\delta }{I^{\ast }}\\ \left(\frac{{\alpha }_1{C}^{\ast }}{k_1+C^{\ast }}\right)\left(\frac{2k_2{\alpha }_3{C}^{\ast }{I}^{\ast }}{{\left(k_2+{C^{\ast }}^2\right)}^2}{-\alpha }_2{I}^{\ast }\right)>\left(r_1-2r_1{r}_2{C}^{\ast }-\frac{k_1{\alpha }_1{I}^{\ast }}{{\left(k_1+C^{\ast }\right)}^2}-{\beta }_1\right)\left(\frac{\delta }{I^{\ast }}\right)\end{array}\right\}$$

(6)

Then,$S_2$is locally asymptotically stable.

Proof. 

The Jacobian matrix at $S_2$ is computed, and it is given as

$$J\left(S_2\right)=\left[\begin{array}{cc}{r}_1-2r_1{r}_2{C}^{\ast }-\frac{k_1{\alpha }_1{I}^{\ast }}{{\left(k_1+C^{\ast }\right)}^2}-{\beta }_1& \frac{-{\alpha }_1{C}^{\ast }}{k_1+C^{\ast }}\\ {-\alpha }_2{I}^{\ast }+\frac{2k_2{\alpha }_3{C}^{\ast }{I}^{\ast }}{{\left(k_2+{C^{\ast }}^2\right)}^2}& -\frac{\delta }{I^{\ast }}\end{array}\right]$$

Then, computing $\left|J\left(S_2\right)-I\lambda \right|=0$ gives:

$${\lambda }^2-Tr\left(J\left(S_2\right)\right)\lambda +det\left(J\left(S_2\right)\right)=0$$

(7)

where

$$Tr\left(J\left(S_2\right)\right)=r_1-2r_1{r}_2{C}^{\ast }-\frac{k_1{\alpha }_1{I}^{\ast }}{{\left(k_1+C^{\ast }\right)}^2}-{\beta }_1-\frac{\delta }{I^{\ast }}$$

$$det\left(J\left(S_2\right)\right)=\left(r_1-2r_1{r}_2{C}^{\ast }-\frac{k_1{\alpha }_1{I}^{\ast }}{{\left(k_1+C^{\ast }\right)}^2}-{\beta }_1\right)\left(-\frac{\delta }{I^{\ast }}\right)+\left(\frac{{\alpha }_1{C}^{\ast }}{k_1+C^{\ast }}\right)\left(\frac{2k_2{\alpha }_3{C}^{\ast }{I}^{\ast }}{{\left(k_2+{C^{\ast }}^2\right)}^2}{-\alpha }_2{I}^{\ast }\right)$$

Therefore, $S_2$ is locally asymptotically stable if condition (6) is satisfied. □

The above analysis shows that the steady state $S_2$ of the CICV model is unstable if the immune system is weak. In this situation, tumor cells have the potential to divide and multiply at a rapid rate.

Global stability implies that all routes with positive initial conditions eventually drift to the system’s attractor. The following two theorems address the global dynamics of the CICV model.

Theorem 5. 

$S_1$is globally asymptotically stable if the following requirement is met

$$\left.\begin{array}{c}{r}_1<{\beta }_1\\ {\left(\frac{{\alpha }_{3}C}{C^2+k_2}{-\alpha }_2\right)}^2\le \frac{4r_1{\beta }_2}{I}\end{array}\right\}$$

(8)

Proof. 

Let ${\rm Y}\left(C,I\right)=C+\left[I-\widehat{I}-\widehat{I}\mathit{ln}\left(\frac{I}{\widehat{I}}\right)\right]$, which is a positive function on $D=\left\{\left(C,I\right){\in R}_{+}^2:C\ge 0,I>0\right\}$. Thus,

$$\frac{d{\rm Y}}{dt}=\frac{dC}{dt}+\left(\frac{I-\widehat{I}}{I}\right)\frac{dI}{dt}=\left(r_1-{\beta }_1\right)C-\frac{{\alpha }_{1}CI}{k_1+C}-\frac{{{\beta }_2\left(I-\widehat{I}\right)}^2}{I}+\left[\frac{{\alpha }_{3}C}{C^2+k_2}{-\alpha }_2\right]C\left(I-\widehat{I}\right)-r_1{r}_2{C}^2$$

i.e.,

$$\frac{d{\rm Y}}{dt}\le \left(r_1-{\beta }_1\right)C-\frac{{\alpha }_{1}CI}{k_1+C}-{\left[\sqrt{\frac{{\beta }_2}{I}}\left(I-\widehat{I}\right)-\sqrt{r_1{r}_2}C\right]}^2$$

Then, $\frac{d{\rm Y}}{dt}<0$ under condition (8). Hence, ${\rm Y}$ is a Lyapunov function [28]. Consequently, $S_1$ is globally asymptotically stable in $D$ if $\left(C,I\right)$ is restricted as in condition (8). □

Theorem 6. 

Suppose that one of the following conditions is satisfied

$$\left.\begin{array}{c}\frac{{\alpha }_1}{{\left(k_1+C\right)}^2}<\frac{\left(\delta +\gamma \right)}{C{I}^2}+\frac{r_1{r}_2}{I}\\ \frac{{\alpha }_1}{{\left(k_1+C\right)}^2}>\frac{\left(\delta +\gamma \right)}{C{I}^2}+\frac{r_1{r}_2}{I}\end{array}\right\}$$

(9)

Then,$S_2$is globally asymptotically stable whenever it exists.

Proof. 

For any initial value $(C,I)$ in the interior of $R_{+}^2$, let $Z\left(C,I\right)=\frac{1}{CI}$. Clearly, $Z\left(C,I\right)>0$, and it is a $C^1$ function for all $(C,I)$ in the interior of $R_{+}^2$. Assume that

$$Z_1\left(C,I\right)=r_{1}C\left(1-r_{2}C\right)-\frac{{\alpha }_{1}CI}{k_1+C}-{\beta }_{1}C$$

$$Z_2\left(C,I\right)=\delta -{\alpha }_{2}CI+\frac{{\alpha }_3{C}^{2}I}{C^2+k_2}-{\beta }_{2}I+\gamma$$

Now, since

$${ZZ}_1\left(C,I\right)=\frac{r_1}{I}-\frac{r_1{r}_{2}C}{I}-\frac{{\alpha }_1}{k_1+C}-\frac{{\beta }_1}{I}$$

$${ZZ}_2\left(C,I\right)=\frac{\delta +\gamma }{CI}-{\alpha }_2+\frac{{\alpha }_{3}C}{C^2+k_2}-\frac{{\beta }_2}{C}$$

we obtain,

$$\mathrm{\Delta }\left(C,I\right)=\frac{\partial {ZZ}_1}{\partial C}+\frac{\partial {ZZ}_2}{\partial I}=-\frac{r_1{r}_2}{I}+\frac{{\alpha }_1}{{\left(k_1+C\right)}^2}-\frac{\left(\delta +\gamma \right)}{C{I}^2}$$

It is obvious that $\mathrm{\Delta }\left(C,I\right)\ne 0$, and it does not change the sign if one of the conditions in Equation (9) is met. Then, according to the Bendixson–Dulac criterion [29], there is no periodic solution in $R_{+}^2$. Since all the solutions of the CICV model are bounded and $S_2$ is the only interior steady state, by using the Poincare–Bendixson theorem [28], $S_2$ is globally asymptotically stable.

Persistence denotes the future survival of all system populations. Now, the average Lyapunov function approach [30] is used to investigate the persistence of the CICV model. □

Theorem 7. 

Assume that the following conditions are satisfied

$$r_1>\frac{{\alpha }_1\left(\delta +\gamma \right)}{k_1{\beta }_2}+{\beta }_1$$

(10)

Then, system (1) is uniformly persistent.

Proof. 

Define $\omega \left(C,I\right)=C^x{I}^y,$ where $x,y$ are positive constants. Clearly $\omega \left(C,I\right)>0$ for all $\left(C,I\right)\in R_{+}^2$, and $\omega \left(C,I\right)\to 0$ when one of the variables $C$ or $I$ approaches zero. Consequently, direct computation gives:

$$\phi \left(C,I\right)=\frac{{\omega }^{\prime }\left(C,I\right)}{\omega \left(C,I\right)}=x\left[r_1\left(1-r_{2}C\right)-\frac{{\alpha }_{1}I}{k_1+C}-{\beta }_1\right]+y\left[\frac{\delta +\gamma }{I}-{\alpha }_{2}C+\frac{{\alpha }_3{C}^2}{C^2+k_2}-{\beta }_2\right]$$

Now, we have

$$\phi \left(S_1\right)=x\left[r_1-\frac{{\alpha }_1\left(\delta +\gamma \right)}{k_1{\beta }_2}-{\beta }_1\right]$$

Hence, according to condition (10), the CICV model is uniformly persistent. □

5. Local Bifurcation

This section investigates the local bifurcation conditions near stable steady states using Sotomayor’s theorem [30]. For this purpose, the CICV model can be rewritten in the following vector form.

$$\frac{d\vartheta }{dt}=F(\vartheta )$$

with

$$\vartheta =\left[\begin{array}{c}C\\ I\end{array}\right]$$

and

$$F=\left[\begin{array}{c}{f}_1\left(C,I\right)\\ f_2\left(C,I\right)\end{array}\right]$$

Now, the Jacobian matrix at any point is given by Equation (4). Then, for any nonzero vector $h={(h_1,h_2)}^T$:

$$DF\left(h,h\right)=\left[\begin{array}{c}\left(r_1-2r_1{r}_{2}C-\frac{k_1{\alpha }_{1}I}{{\left(k_1+C\right)}^2}-{\beta }_1\right)h_1-\frac{{\alpha }_1{h}_{2}C}{k_1+C}\\ \left({\alpha }_{2}I+\frac{2k_2{\alpha }_{3}CI}{{\left(k_2+C^2\right)}^2}\right)h_1+\left(-{\alpha }_{2}C+\frac{{\alpha }_3{C}^2}{k_2+C^2}-{\beta }_2\right)h_2\end{array}\right]$$

and

$$D^{2}F\left(h,h\right)=\left[\begin{array}{c}\left(-2r_1{r}_{2}C+\frac{2k_1{\alpha }_{1}I}{{\left(k_1+C\right)}^3}\right)h_1-\frac{k_1{\alpha }_1{h}_2}{{\left(k_1+C\right)}^2}-\frac{{\alpha }_1{h}_{2}C}{k_1+C}\\ \left({\alpha }_{2}I+\frac{2k_2{\alpha }_{3}CI}{{\left(k_2+C^2\right)}^2}\right)h_1+\left({-\alpha }_{2}C+\frac{{\alpha }_3{C}^2}{k_2+C^2}-{\beta }_2\right)h_2\end{array}\right].$$

(11)

Theorem 8. 

For the parameter value${\beta }_1^{\ast }=r_1-\frac{{\alpha }_1\left(\delta +\gamma \right)}{k_1{\beta }_2}$, the CICV model system, at$S_1$, has a trans-critical bifurcation.

Proof. 

At ${\beta }_1^{\ast }=r_1-\frac{{\alpha }_1\left(\delta +\gamma \right)}{k_1{\beta }_2}$, $J\left(S_1\right)$ has a zero eigenvalue ${\lambda }_{11}=0$. Therefore, $J\left(S_1\right)$ at ${\beta }_1={\beta }_1^{\ast }$ becomes

$$J^{\ast }\left(S_1\right)=\left[\begin{array}{cc}0& 0\\ \frac{{-\alpha }_2\left(\delta +\gamma \right)}{{\beta }_2}& -{\beta }_2\end{array}\right]$$

Let $h^{\left[1\right]}={(h_1^{\left[1\right]},h_2^{\left[1\right]})}^T$ be the eigenvector corresponding to ${\lambda }_{11}=0$. Then, $\left(J^{\mathrm{*}}\left(S_1\right)-{\lambda }_{11}I\right)h^{\left[1\right]}=0$ gives ${\left(h_1^{\left[1\right]},h_2^{\left[1\right]}\right)}^T={\left(1,\frac{-{\alpha }_2(\delta +\gamma )}{{\beta }_2^2}\right)}^T$.

Now, let $v^{\left[1\right]}={(v_1^{\left[1\right]},v_2^{\left[1\right]})}^T$ be the eigenvector corresponding to the eigenvalue ${\lambda }_{11}=0$ of the matrix ${{(J}^{\ast }\left(S_1\right))}^T$. Then, $({{(J}^{\ast }\left(S_2\right))}^T-{\lambda }_{11}I)v^{\left[1\right]}=0$. Then, the direct calculation gives ${(v_1^{\left[1\right]},v_2^{\left[1\right]})}^T={(\mathrm{1,0})}^T$.

Subsequently, the following is taken into account to verify that the requirements of Sotomayor’s theorem for trans-critical bifurcation are obtained:

$$\frac{\partial F}{\partial {\beta }_1}=F_{{\beta }_1}\left(\vartheta ,{\beta }_1\right)={\left(\frac{\partial f_1}{\partial {\beta }_1},\frac{\partial f_2}{\partial {\beta }_1}\right)}^T={\left(-C,0\right)}^T$$

Therefore, $F_{{\beta }_1}\left(S_1,{\beta }_1^{\ast }\right)={\left(\mathrm{0,0}\right)}^T,$ and hence, ${\left(v^{\left[1\right]}\right)}^T{F}_{{\beta }_1}\left(S_1,{\beta }_1^{\ast }\right)=0$. So, the first condition of trans-critical bifurcation is met.

Now,

$${DF}_{{\beta }_1}\left(\vartheta ,{\beta }_1\right)=\left[\begin{array}{cc}-1& 0\\ 0& 0\end{array}\right]$$

where ${DF}_{{\beta }_1}\left(\vartheta ,{\beta }_1\right)$ denotes the derivative of $F_{{\beta }_1}\left(\vartheta ,{\beta }_1\right)$ with respect to $\vartheta ={\left(C,I\right)}^T$.

Further,

$${DF}_{{\beta }_1}\left(S_1,{\beta }_1^{\ast }\right)h^{\left[1\right]}=\left[\begin{array}{cc}-1& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}1\\ \frac{-{\alpha }_2(\delta +\gamma )}{{\beta }_2^2}\end{array}\right]=\left[\begin{array}{c}-1\\ 0\end{array}\right]$$

$${\left(v^{\left[1\right]}\right)}^T\left[{DF}_{{\beta }_1}\left(S_1,{\beta }_1^{\ast }\right)h^{\left[1\right]}\right]=\left(\mathrm{1,0}\right){\left(-\mathrm{1,0}\right)}^T=-1\ne 0$$

Now, by substituting in (11), it is found that

$$D^{2}F\left(S_1,{\beta }_1^{\ast }\right)\left(h_1^{\left[1\right]},h_2^{\left[1\right]}\right)=\left[\begin{array}{c}\frac{2{\alpha }_1(\delta +\gamma )}{{{\beta }_{2}k}_1^2}+\frac{{\alpha }_1{\alpha }_2(\delta +\gamma )}{k_1{\beta }_2^2}\\ \frac{{\alpha }_2(\delta +\gamma )}{{\beta }_2}+\frac{{\beta }_2{\alpha }_2(\delta +\gamma )}{{\beta }_2^2}\end{array}\right]$$

Hence,

$$\begin{gathered} {{\left(v^{\left[1\right]}\right)}^{T}D}^{2}F\left(S_1,{\beta }_1^{\ast }\right)\left(h_1^{\left[1\right]},h_2^{\left[1\right]}\right)=\left(\mathrm{1,0}\right){\left(\frac{2{\alpha }_1(\delta +\gamma )}{{{\beta }_{2}k}_1^2}+\frac{{\alpha }_1{\alpha }_2(\delta +\gamma )}{k_1{\beta }_2^2},\frac{{\alpha }_2(\delta +\gamma )}{{\beta }_2}+\frac{{\beta }_2{\alpha }_2(\delta +\gamma )}{{\beta }_2^2}\right)}^T \\ =\frac{2{\alpha }_1(\delta +\gamma )}{{\beta }_2{k}_1^2}+\frac{{{\alpha }_1\alpha }_2(\delta +\gamma )}{k_1{\beta }_2^2}\ne 0 \end{gathered}$$

Due to Sotomayor’s local bifurcation theorem, the CICV model has a trans-critical bifurcation at $S_1$ with $\beta ={\beta }_1^{\ast }$. □

Theorem 9. 

Assume that the following conditions are satisfied

$$\left.\begin{array}{c}{r}_1>2r_1{r}_2{C}^{\ast }+{\beta }_1+\frac{\delta }{I^{\ast }},\\ det\left(J\left(S_2\right)\right)>0.\end{array}\right\}$$

(12)

Then the CICV model has a Hopf bifurcation at${\alpha }_1={\alpha }_1^{\ast }$.

Proof. 

Consider the characteristic equation at $S_2$ which is given in (7). To validate the conditions for a Hopf bifurcation, we need to verify that $Tr\left(J\left(S_2\right)\right)=0$ is satisfied. It is detected that $Tr\left(J\left(S_2\right)\right)=0$ gives:

$$r_1-2r_1{r}_2{C}^{\ast }-\frac{k_1{\alpha }_1{I}^{\ast }}{{\left(k_1+C^{\ast }\right)}^2}-{\beta }_1-\frac{\delta }{I^{\ast }},$$

$${\alpha }_1^{\ast }=\frac{{\left(k_1+C^{\ast }\right)}^2}{k_1{I}^{\ast }}\left[r_1-2r_1{r}_2{C}^{\ast }-{\beta }_1-\frac{\delta }{I^{\ast }}\right]$$

Clearly, ${\alpha }_1^{\ast }>0$ provided that the first inequality of condition (12) holds. Now, at ${\alpha }_1={\alpha }_1^{\ast },$ the characteristic equation given by Equation (6) is rewritten as

$${\lambda }^2+det\left(J\left(S_2\right)\right)=0,$$

which has two roots

$${\lambda }_{\mathrm{1,2}}=\pm i\sqrt{det\left(J\left(S_2\right)\right)}$$

Clearly, at ${\alpha }_1={\alpha }_1^{\ast }$, there are two purely imaginary eigenvalues ${\lambda }_1$ and ${\lambda }_2$ which are complex conjugates if the second inequality of condition (12) holds. Further, for all values of ${\alpha }_1$ in a neighborhood of ${\alpha }_1^{\ast }$, the roots generally are given by the following formula:

$${\lambda }_{\mathrm{1,2}}=\frac{trac\left(J\left(S_2\right)\right)\pm i\sqrt{det\left(J\left(S_2\right)\right)}}{2}$$

Further, due to the transversality condition

$$\frac{d}{d{\alpha }_1}{\left[Re\left({\lambda }_{\mathrm{1,2}}\right)\right]}_{{\alpha }_1={\alpha }_1^{\ast }}=\frac{d}{d{\alpha }_1}{\left[\frac{trac\left(J\left(S_2\right)\right)}{2}\right]}_{{\alpha }_1={\alpha }_1^{\ast }}=-\frac{k_1{I}^{\ast }}{{\left(k_1+C^{\ast }\right)}^2}\ne 0$$

the CICV model has a Hopf bifurcation at ${\alpha }_1={\alpha }_1^{\ast }$. □

6. Numerical Simulations

Numerical simulations are carried out in this section to show various dynamic situations. A fourth-order Runge–Kutta method is used via the ode45 command in MATLAB R2021b to attain stable or unstable equilibrium solutions or convergent solutions for the CICV model. The simulations of the CICV model were performed over a time interval of ninety days, with the parameters defined in Table 1.

Now, four cases will be taken into account to understand the dynamic behavior of the CICV model and evaluate the impact of chemotherapy treatment on tumor suppression. Then the results of the four cases will be compared. The four cases are

1. 

Dynamic behavior of the CICV model without vitamins and chemotherapy.

In this case, we study the dynamics of the interaction between tumor cells $C(t)$ and immune cells $I(t)$ when no external therapy is applied (${\beta }_1=\gamma =0)$. We compare various situations with treatment to the case without therapy. We also estimate the amount of minimum treatment required to eliminate cancer. Figure 2 depicts the behavior of the CICV model for the data given in Table 1, with ${\beta }_1=\gamma =0$. It shows that the CICV model has the tumor-free equilibrium point $S_1=(\mathrm{0,1.21})$ and the unique positive equilibrium $S_2=(\mathrm{2.24,0.97})$. Moreover, for a variety of initial values, the solution first begins to increase or decrease for a certain amount of time before it eventually settles down asymptotically to $S_2$ after about thirty days. Further, the number of immune cells gradually decreases as the number of tumor cells gradually grows. On the other hand, $S_1$ shows saddle behavior. In light of this, it is also abundantly evident from Case 1 that eliminating tumor cells is impossible without a treatment plan.

2. 

Dynamic behavior of the CICV model with vitamins.

In this case, we study the dynamics of the CICV model when regular vitamin consumption is implemented to strengthen the immune system. Figure 3 illustrates the behavior of the CICV model for the data given in Table 1, with ${\beta }_1=0$ (without chemotherapy drug). It shows the solutions for all initial conditions reach the interaction state $S_2=(\mathrm{1.88,1.97})$ after about forty days. Further, the number of immune cells, in this case, gradually increases as the number of tumor cells decreases. Even though there is a significant reduction in the number of tumor cells compared to Case 1, the immune system still cannot eliminate all tumor cells.

3. 

Dynamic behavior of the CICV model with chemotherapy.

In this case, we are going to discuss the dynamics of the CICV model if external therapy (chemotherapy) is applied without vitamin consumption $(\gamma =0)$. Figure 4 clearly describes the global stability behavior of the positive steady state $S_2=(\mathrm{0.62,1.97})$ of the CICV model. Further, it can be concluded that after about forty-five days, the positive steady state is reached. Tumor cells are significantly reduced in the body with the use of chemotherapy compared with the previous two cases. Additionally, chemotherapy harms the immune cell as well; we can see a reduction in the number of immune cells compared with Case 2. In view of the above, more doses are required to reach the tumor-free state.

4. 

Dynamic behavior of the CICV model with chemotherapy and vitamins.

Finally, in the last case, the impact of vitamins and chemotherapy is applied to the dynamics of the CICV model. We simulate the CICV model with the parameter values presented in Table 1. It is clear from Figure 5 that taking chemotherapy in combination with vitamins succeeds in the clearance of tumor cells after about fifty days. There exists only the free-tumor equilibrium point $S_1=(0,2.15)$ with nodal sink behavior. Additionally, it can be observed that the level of immune cells rises significantly after taking the immune system booster. As a result, the CICV model with the parameters given in Table 1, including ${\beta }_1=0.2$ and $\gamma =0.5463$, loses the persistence, and the tumor-free state shows an asymptotically stable behavior. See Figure 5. In view of the above, the patient can reach a healthy state when a combination of vitamin intake and chemotherapy is applied.

The second target of the computational simulation is to determine the minimum consumption of vitamins and chemotherapy required to reach a healthy state. Figure 6 shows the impact of varying the parameters $\gamma$ on the CICV model. It is clear that by varying the values of $\gamma$ and keeping the rest of the parameters as in Table 1, the solution asymptotically approaches the tumor-free state at $\gamma \ge 0.35$. In contrast, the solution of the CICV model converges to the interaction state at $\gamma \le 0.34$. As a result, $\gamma =0.34$ is the minimum number of vitamins that should be consumed each day to eliminate cancer.

Finally, the impact of altering the number of chemotherapy doses is determined in Figure 7. It is clear that the trajectory of the CICV model settles down asymptotically to the interaction state $S_2$ for ${\beta }_1\le 0.12$. On the other hand, the system loses persistence and approaches the tumor-free state $S_1$ for ${\beta }_1\ge 0.14$. Therefore, ${\beta }_1=0.14$ shows the lowest dose of chemotherapy required to reach a cancer-free state.

7. Discussion

The CICV model dynamics have been considered to investigate the impacts of regular intakes of vitamins and chemotherapy on the dynamics of tumor-immune interactions. The theoretical study showed that the CICV model has two main steady states: the free-tumor and interaction equilibria. Depending on the choice of parameter values, the two equilibrium points can show stable, unstable, or saddle point behavior. We derived results both on the local and global behavior of the equilibria. The numerical simulations confirm the analytical results. In particular, the threshold values for the trans-critical bifurcation are computed, which shows the transition between the persistence of cancer and its eradication. Let us now compare the simulation results of our CICV model with the numerical results presented in [26]. When comparing the results in [26] with our results, the critical rate of vitamin intake required to strengthen the immune system so much that it leads to the elimination of cancer is the most important metric. The outcome of the numerical simulation of the model for vitamin intake but without chemotaxis, presented in [26], suggested that a person’s immune system can be strengthened enough to eliminate cancer by taking vitamins consistently at a rate of about 55% each day. However, in our system, as shown in the numerical simulation (Case 2), the regular use of vitamins, equivalent to 55% of the recommended daily allowance, is insufficient to eliminate malignant cells. On the other hand, when vitamin intake is integrated with chemotherapy, only a rate of 35% per day of vitamin intake is required to eliminate the cancer cells in the body completely.

8. Conclusions

This study aims to discover the conditions that lead to eliminating tumor cells in cancer–immune–chemotherapy–vitamins model. Through the analysis of the CICV model, the existence of the equilibrium points and their corresponding stability conditions has been determined. For a specific parameter set, it has been found that the model may have two equilibrium states. The first is the tumor-free equilibrium state, meaning that tumor cells will be eliminated. The second is the coexisting equilibrium point, which proposes that tumor cells and immune cells will coexist with nonzero populations. In this context, the stability of the tumor-free equilibrium is critical. The stability of the tumor-free equilibrium point means that the treatment is successful since the time-dependent solution will reach a cancer-free state. Numerical simulations highlight the importance of boosting the immune cells and taking chemotherapy to eliminate tumor cells. We considered the intake of vitamins and chemotherapy both individually and in combination, and we established the thresholds required for reaching a healthy state. Mathematically, these thresholds correspond to a trans-critical bifurcation. In summary, this study shows that applying regular doses of chemotherapy and taking vitamins can promote the immune system and inhibit and delay tumor cell growth and division, respectively. The successful implementation of the results in this paper might lead to treatment strategies which will help oncologists in practicing cancer treatment. In the future, we plan to generalize the CICV model considered here in different directions. We will consider a delay differential equation to examine the delayed effects of treatments on boosting immune cells and suppressing tumor cell growth. In addition, we will add radiation therapy to the model by extending the model to a three-component system.