Stability Analysis of a Targeted Chemotherapy-Cancer Delayed Model

Zahra Mousa Albariqi, Lama Shujaa Aljoufi, Rawan Ali Aldosary, Ashwaq Atallah Alharbi, Hanan Ali Batarfi

Department of Mathematics, Faculty of Science, King Abdulaziz University, Jeddah, Saudi Arabia.

DOI: 10.4236/ajcm.2025.152006                

Abstract

In this paper, we investigate the dynamic behaviour of a mathematical model of cancer that includes immune cells, tumor cells, and normal cells, and explore the effects of the introduction of a delayed term of targeted therapy on the model. This model was first proposed by Anusmita Das et al., numerous studies have attempted to model the interaction between tumours and the immune system using deterministic delay differential equations (DDEs) so a delay term was added, in this paper, on the basis to make the model more realistic. Also, the local and global stability of the equilibrium point of the model is analyzed by linearization and Lyapunov method, and the numerical simulation of MATLAB is used to verify the analysis results.

Keywords

Mathematical Model of Cancer, Delay, Stability Analysis, Chemotherapy, Numerical Simulation, Delay Stability Analysis

1. Introduction

Cancer is a large group of diseases that can start in almost any organ or tissue of the body when abnormal cells grow uncontrollably, go beyond their usual boundaries to invade adjoining parts of the body and/or spread to other organs. The latter process is called metastasizing and is a major cause of death from cancer. A neoplasm and malignant tumour are other common names for cancer. Characterized by the uncontrolled growth and spread of abnormal cells, is a major public health challenge and the second leading cause of death worldwide [1]. The cancer burden continues to grow globally, exerting tremendous physical, emotional and financial strain on individuals, families, communities and health systems. Many health systems in low- and middle-income countries are least prepared to manage this burden, and large numbers of cancer patients globally do not have access to timely quality diagnosis and treatment. In countries where health systems are strong, survival rates of many types of cancers are improving thanks to accessible early detection, quality treatment and survivorship care.

Extensive research efforts have been dedicated to understanding the complex mechanisms underlying cancer development, progression, and metastasis, with the ultimate goal of developing effective treatment strategies. Mathematical modeling has emerged as a powerful tool in cancer research, enabling researchers to gain insights into the intricate dynamics of tumor growth, interactions with the immune system, and responses to various therapeutic interventions. Those Mathematical models simulate complex systems in a relatively fast time without the enormous costs of laboratory experiments and the corresponding biological variations. In particular for oncology, such models can be calibrated using experimental or clinical data [2]-[4].

Cancer is the second leading cause of death in humans according to WHO. Many Medical research centers are trying to develop new treatments and strategies to fight cancer by understanding the dynamic of the growth and the interaction of the tumor growth with the cells and the immune system [5]-[9]. This study together with the different types of treatment was heavily modelled mathematically helping medics to develop treatments and ways to control the spreading of cancer. Over the years, numerous mathematical models have been proposed to study various aspects of cancer, including tumor evolution, angiogenesis, immune response, and the effects of different treatment modalities. These models have incorporated various factors such as tumor cell proliferation, immune cell interactions, and drug pharmacokinetics and pharmacodynamics. Some notable examples include the work of Allison et al. [10] and Liu et al. [11] on modeling tumor growth rate, carrying capacity, and cytolytic activity of effector cells; Li et al. [12] on the effects of angiogenic growth factors; and De Pillis et al. [13]-[15] on chemotherapy, immunotherapy, and combination therapies. The main problem in modelling the cancer treatment or spread is the necessity to obtain data that define the values of the biological parameters used within the mathematical model.

One of the key challenges in accurately modeling cancer dynamics is accounting for the time delays inherent in biological processes, such as the immune system’s response time and the lag between drug administration and therapeutic effect. Time delays can significantly impact the behavior of mathematical models and their stability properties. Several researchers have incorporated time delays into their cancer models to better reflect the real-world dynamics. For instance, Li et al. [12] explored the effects of time delays on angiogenic growth, while Arabameri et al. [16] incorporated time delays in a dendritic cell-based immunotherapy model.

At the beginning, the study was done within a time range with no information passed from the ‘past’. For the system to be realistic, numerous studies have attempted to model the interaction between tumours and the immune system using delay differential equations (DDEs) (see for example [8]), which has an important role in the mathematical modelling of multi-species interactions. The time-delay Kuznetsor and Taylor model was studied by Galach [8] to achieve better consistency with reality. Dehingia et al. [17] studied the tumor-immune cells’ interaction with a delayed time, followed by Anusmita et al. [18] who studied the dynamical behavior of a mathematical model of cancer including tumor-cells, immune-cells, and normal-cells is investigated when a delay term is induced.

In this research, we modified the work of Anusmita et al. [19] and analyzed the delay-induced mathematical model of cancer that incorporates time delays to the more realistic simulation of the immune system’s behavior and its interactions with tumor cells. By studying the effects of these delays on the model dynamics and stability properties through mathematical and numerical analysis, we aim to gain insights into the complex interplay between cancer progression, immune response, and potential therapeutic interventions.

The manuscript is structured as follows: Section 2, describes the work of Anusmita et al. [19] in the normal case where the model is constructed as a system of ordinary differential equations (ODEs) together with the basic assumptions and results. Section 3, we stated the suggested constant-delay model and examined the positivity and boundedness of the model’s solution. Section 4, the existence of equilibrium points and local stability for the system were discussed. Section 5, we proved the global stability of the tumor-free state. Finally, Section 6, the simulation of the solution is done and analysed.

2. The Ordinary Differential Equation (ODE) Model

The ODE model presented by Anusmita Das et al. [19] shows that a prescribed treatment can eradicate tumor cells from the body for a threshold value of tumor growth rate. The numerical simulations in their paper, also, showed 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 the model is that when the tumor size is large and, consequently, requires prolonged treatment involving a high dosage or variety of drugs, this can harm the patient’s body.

The system of Anusmita Das et al was described by the ODEs;

$I^{\prime }=s+\frac{\rho I\left(t\right)T\left(t\right)}{\sigma +T\left(t\right)}-c_{1}IT-d_{1}I-a_1\left(1-\eta \right)CI,$ (1a)

$T^{\prime }\mathrm{<mo>=</mo>}{r}_{1}T\left(1-b_{1}T\right)-c_{2}IT-c_{3}TN-a_{2}CT\mathrm{,}$ (1b)

$N^{\prime }=r_{2}N\left(1-N\right)-c_{4}TN-a_3\left(1-\eta \right)CN\mathrm{,}$ (1c)

$C^{\prime }=u-d_{2}C-kTC\mathrm{,}$ (1d)

with the initial conditions $I\left(0\right)=I_0>0$ , $T\left(0\right)=T_0>0$ , $N\left(0\right)=N_0>0$ and $C\left(0\right)=C_0>0$ . $I=I\left(t\right)$ , $T=T\left(t\right)$ , $N=N\left(t\right)$ , and $C=C\left(t\right)$ , for $t\ge 0$ . The system represents the change of densities with time $t$ of the effector cells $I\left(t\right)$ , tumor cells $T\left(t\right)$ , healthy-cells $N\left(t\right)$ and the amount of targeted chemo drug $C\left(t\right)$ administered at time $t$ . The number of normal-cells was normalized by taking their carrying capacity equal to one.

The first term in (1a), namely $s$ , represents the constant source rate of effector-cells. Effector-cells are recruited by tumor-cells through the Michaelis-Menten term $\rho IT/\left(\sigma +T\right)$ , where $\rho$ is the rate at which the effector-cells grow, and $\sigma$ represents the steepness of effector is response.

In the second Equation (1b), tumor-cells are assumed to grow logistically with an intrinsic growth rate $r_1$ and the maximum carrying capacity of $b_1$ in the absence of immune-cells and drug therapy. Tumor-cells are killed by the interaction with immune-cells, normal-cells and targeted chemo-drug are represented by the terms $-c_{2}IT$ , $-c_{3}TN$ and $-a_{2}CT$ .

In the third Equation (1c), normal-cells are also assumed to grow logistically with an intrinsic growth rate $r_2$ and the maximum carrying capacity of one. The second and third terms, $-c_{2}NT$ , $-a_3\left(1-\eta \right)NC$ , represent the kill rate of the normal-cell due to interaction with tumor-cells and targeted chemo-drug.

Finally, in the fourth Equation (1d), $u$ represents the chemotherapy-treatment. The second term, $d_2$ , represents decay rate of targeted chemo-drug, while the third term, $-KTC$ , represents the rate of the attachments of targeted chemo-drugs with tumor-cells. The fourth term indicates that effector-cells die off at a rate of $d_1$ per day and the fifth term $-a_1\left(1-\eta \right)CI$ , represents kill rate of effector-cell by targeted chemo-drug.

The authors proved that the model (1a-1d) solution is positive and bounded in the domain;

$\Delta =\left\{\left(I\mathrm{,}T\mathrm{,}N\mathrm{,}C\right)\in {ℝ}_{+}^4\mathrm{<mo>|</mo>}I\left(t\right)\le s/d_1\mathrm{,}T\left(t\right)\le 1/b_1\mathrm{,}N\left(t\right)\le \mathrm{1,}C\left(t\right)\le u/d_2\right\}\mathrm{.}$

The tumor free and co-axial equilibrium points $E_1$ and $E^{\mathrm{<mo>*</mo>}}$ (respectively) have been evaluated symbolically and proven to be locally stable. Then, the global stability of the tumor-free equilibrium point was investigated and showed that the targeted chemotherapeutic treatment kills the (small size) tumor cells under a certain condition.

3. The Time-Delay Model

Understanding why and how delays affect a model is critical, as they can significantly influence system stability and dynamic behavior. To describe the time lag by the effector system for developing a suitable response after recognizing the tumor-cells, we need to include the effect of time-delay $\tau$ (constant) into the Michaelis-Menten term [20]. Hence, a discrete-time delay is added to the second term of the first equation of system (1a) and, as a result, the term becomes $\rho I\left(t-\tau \right)T\left(t-\tau \right)/\left[\sigma +T\left(t-\tau \right)\right]$ . The third term, $-c_{1}IT$ , represents the kill rate of effector-cells due to interaction with tumor-cells, at time $t$ .

Our modification takes the following form: The model presented in this paper is a modified form of that proposed by Anusmita Das et al. [19]. The modification is carried out by inducing a delay term in the second term (Michaelis-Menten term) of the first equation of the system with proper biological justifications. The change is made to see the impact of the delay in the complete system. Our modification takes the following form:

$I^{\prime }=s+\frac{\rho I\left(t-\tau \right)T\left(t-\tau \right)}{\sigma +T\left(t-\tau \right)}-c_{1}IT-d_{1}I-a_1\left(1-\eta \right)CI,$ (2a)

$T^{\prime }=r_{1}T\left(1-b_{1}T\right)-c_{2}IT-c_{3}TN-a_{2}CT\mathrm{,}$ (2b)

$N^{\prime }=r_{2}N\left(1-N\right)-c_{4}TN-a_3\left(1-\eta \right)CN,$ (2c)

$C^{\prime }=u-d_{2}C-kTC\mathrm{,}$ (2d)

4. Positive Invariance and Boundedness

Using standard comparison theory [21], on the equations of system (2a-2d), we get;

$\frac{\text{d}I\left(t\right)}{\text{d}t}=s+\frac{\rho I\left(t-\tau \right)T\left(t-\tau \right)}{\sigma +T\left(t-\tau \right)}-c_{1}IT-d_{1}I-a_1\left(1-\eta \right)CI\le s-d_{1}I,$

integrating the above leads to;

$I\left(t\right)\le \frac{s}{d_1}+I\left(0\right){\text{e}}^{-d_{1}t},\text{i}\text{.e}.\underset{t\to \infty }{\lim }\sup \left[I\left(t\right)\right]\le \frac{s}{d_1}.$

Furthermore;

$\frac{\text{d}T\left(t\right)}{\text{d}t}=r_{1}T\left(t\right)-r_1{b}_1{T}^2-c_{2}IT-c_{3}TN-a_{3}CT\le r_{1}T-r_1{b}_1{T}^2\mathrm{.}$

Proceeding as above, we have

$T\left(t\right)\le \frac{1}{b_1+T\left(0\right)e^{-r_{1}t}}\mathrm{,}\text{i}\text{.e}\mathrm{.}\underset{t\to \infty }{\lim }\sup \left[T\left(t\right)\right]\le \frac{1}{b_1}\mathrm{.}$

The next equation;

$\frac{\text{d}N\left(t\right)}{\text{d}t}=r_{2}N-r_2{N}^2-c_{4}TN-a_3\left(1-\eta \right)CN\le r_{2}N-r_2{N}^2\mathrm{,}$

$N\left(t\right)\le \frac{1}{b_1+N\left(0\right){\text{e}}^{-r_{2}t}}\mathrm{,}\text{i}\text{.e}\text{.}\underset{t\to \infty }{\lim }\sup \left[N\left(t\right)\right]\le \mathrm{1,}$

and similarly, the last equation we find;

$\frac{\text{d}C\left(t\right)}{\text{d}t}=u-d_{2}C-kTC\le u-d_{2}C\mathrm{.}$

$C\left(t\right)\le \frac{u}{d_2}+C\left(0\right){\text{e}}^{-d_{2}t}\mathrm{,}\text{i}\text{.e}\mathrm{.}\underset{t\to \infty }{\lim }\sup \left[C\left(t\right)\right]\le \frac{1}{d_2}$

with the initial conditions $I\left(0\right)>0$ , $T\left(0\right)>0$ , $N\left(0\right)>0$ and $C\left(0\right)>0$ , for $t\ge 0$ . Using the considered initial values, we assume that $I\left(t\right)\ge 0$ , $T\left(t\right)\ge 0$ , $N\left(t\right)\ge 0$ and $C\left(t\right)\ge 0$ for all $t\ge 0$ . Consequently, the corresponding domain region for the system (2a-2d) is

$\Lambda =\left\{\left(I\mathrm{,}T\mathrm{,}N\mathrm{,}C\right)\in R_{+}^4\mathrm{<mo>|</mo>0}\le I\left(t\right)\le \frac{s}{d_1}\mathrm{,}0\le T\left(t\right)\le \frac{1}{b_1}\mathrm{,}0\le N\left(t\right)\le \mathrm{1,}0\le C\left(t\right)\le \frac{1}{d_2}\right\}\mathrm{.}$

The domain region Λ is a bounded set, and hence the solution of the system (2a-2d) is bounded. The model’s Equations (2a-2d) are subject to the following initial conditions:

$I\left(\xi \right)={\psi }_1\left(\xi \right),T\left(\xi \right)={\psi }_2\left(\xi \right),N\left(\xi \right)={\psi }_3\left(\xi \right),C\left(\xi \right)={\psi }_4\left(\xi \right),$

${\psi }_1\left(\xi \right)\ge \mathrm{0,}{\psi }_2\left(\xi \right)\ge \mathrm{0,}{\psi }_3\left(\xi \right)\ge \mathrm{0,}{\psi }_4\left(\xi \right)\ge 0\text{for}\xi \in \left[-\tau \mathrm{,0}\right]\mathrm{,}$ (3a)

and ${\psi }_1\left(0\right)\ge \mathrm{0,}{\psi }_2\left(0\right)\ge \mathrm{0,}{\psi }_3\left(0\right)\ge \mathrm{0,}{\psi }_4\left(0\right)\ge 0$ , where

$S_{+}=\left\{\left({\psi }_1\left(\xi \right)\mathrm{,}{\psi }_2\left(\xi \right)\mathrm{,}{\psi }_3\left(\xi \right)\mathrm{,}{\psi }_4\left(\xi \right)\right)\in S\left(\left[-\tau \mathrm{,0}\right]\mathrm{,}{R}_{+}^4\right)\right\}\mathrm{.}$

The delay differential system (2a-2d) can be written in the vector form as

$\begin{array}{c}{X}^{\prime }=\left(\begin{array}{c}{I}^{\prime }\left(t\right)\\ T^{\prime }\left(t\right)\\ N^{\prime }\left(t\right)\\ C^{\prime }\left(t\right)\end{array}\right)=\left(\begin{array}{c}s+\frac{\rho I\left(t-\tau \right)T\left(t-\tau \right)}{\sigma +T\left(t-\tau \right)}-c_{1}IT-d_{1}I-a_1\left(1-\eta \right)CI\\ r_{1}T\left(1-T\right)-c_{2}IT-c_{3}TN-a_{3}CT\\ r_{2}N\left(1-N\right)-c_{4}TN-a_3\left(1-\eta \right)CN\\ u-d_{2}C-kTC\end{array}\right)\\ =\left(\begin{array}{c}{V}_1\left(X\right)\\ V_2\left(X\right)\\ V_3\left(X\right)\\ V_4\left(X\right)\end{array}\right)=V\left(X\right),\end{array}$ (3b)

where $V\in C^{\infty }\left(R_{+}^4\right)$ is defined in the positive quadrant $R_{+}^4$ and represents a mapping $V\mathrm{<mo>:</mo>}{S}_{+}\to R_{+}^4$ . The right-hand side of the system (3b) is locally Lipchitz, meaning that the derivative is bounded and satisfies;

${V_i\left(X\right)|}_{Y_i\left(t\right)=0,X\in S_{+}}=V_i\left(0\right),\forall i=1,2,3,4.$

According to the second lemma by Yang et al. [22], every solution of system (3b) with the initial conditions (3a), where ${\psi }_i\left(t\right)\in S_{+}$ , say $Y\left(t\right)=Y\left(t\mathrm{<mo>;</mo>}Y\left(0\right)\right)$ , for all $t>0$ , remains positive throughout the domain $S_{+},\forall t>0$ . Therefore, the solution of (3b) is positively invariant in time $t$ .

5. Stability Analysis

5.1. Tumor-free equilibrium

Free the equilibrium point, the change with time is set to zero, i.e,

$\frac{\text{d}I}{\text{d}t}=\frac{\text{d}T}{\text{d}t}=\frac{\text{d}N}{\text{d}t}=\frac{\text{d}C}{\text{d}t}=0$ which gives the following non-linear algebraic non-homogenous equation;

$0=s+\frac{\rho I_1{T}_1}{\sigma +T_1}-c_1{I}_1{T}_1-d_1{I}_1-a_1\left(1-\eta \right)C_1{I}_1\mathrm{,}$ (4a)

$0=r_1{T}_1\left(1-b_1{T}_1\right)-c_2{I}_1{T}_1-c_3{T}_1{N}_1-a_2{C}_1{T}_1\mathrm{,}$ (4b)

$0=r_2{N}_1\left(1-N_1\right)-c_4{T}_1{N}_1-a_3\left(1-\eta \right)C_1{N}_1\mathrm{,}$ (4c)

$0=u-d_2{C}_1-k{T}_1{C}_1\mathrm{,}$ (4d)

and hence the tumor-free equilibrium: $E_1\left(I_1\mathrm{,}{T}_1\mathrm{,}{N}_1\mathrm{,}{C}_1\right)$ ; (where $T=0$ ) can be evaluated as follows; from Equation (4c);

$r_2{N}_1\left(1-N_1\right)-c_4{T}_1{N}_1-a_3\left(1-\eta \right)C_1{N}_1=\mathrm{0,}$

$r_2{N}_1\left(1-N_1\right)-a_3\left(1-\eta \right)C_1{N}_1=\mathrm{0,}$

$N_1\left[r_2\left(1-N_1\right)-a_3\left(1-\eta \right)C_1\right]=0$

which gives the equilibrium coordinate:

$N_1=0\text{or}{N}_1=\frac{r_2-a_3\left(1-\eta \right)C_1}{r_2}$ (4e)

For $N_1=0$ , implies that the patients will not be alive, so we do not consider those cases where $N_1=0$ . From Equation (4d) we get;

$C_1=\frac{u}{d_2}\mathrm{,}$ (5a)

from Equations (5a) and (4e);

$N_1=\frac{r_2-a_3\left(1-\eta \right)\frac{u}{d_2}}{r_2}\mathrm{.}$ (5b)

Now, taking Equation (4b);

$r_1{T}_1\left(1-b_1{T}_1\right)-c_2{I}_1{T}_1-c_3{T}_1{N}_1-a_2{C}_1{T}_1\mathrm{,}$

$T_1\left[r_1\left(1-b_1{T}_1\right)-c_2{I}_1-c_3{N}_1-a_2{C}_1\right]=\mathrm{0,}$

$T_1=0,T_1=\frac{r_1-\left(c_2{I}_1-c_3{N}_1-a_2{C}_1\right)}{r_1{b}_1},$

$0=\frac{r_1-\left(c_2{I}_1-c_3{N}_1-a_2{C}_1\right)}{r_1{b}_1}\mathrm{,}$

and hence;

$I_1=\frac{r_1-c_3{N}_1-a_2{C}_1}{c_2}\mathrm{.}$

By substituting Equation (5a);

$\begin{array}{c}{I}_1=\frac{1}{c_2}\left[r_1-\frac{c_3}{r_2}\left(r_2-a_3\left(1-\eta \right)\frac{u}{d_2}\right)-\frac{a_{2}u}{d_2}\right]\\ =\frac{1}{c_2}\left[r_1-c_3+\left(\frac{c_3{a}_3\left(1-\eta \right)}{r_2}-a_2\right)\frac{u}{d_2}\right],\end{array}$ (5c)

and hence from Equations (5a) the tumor-free equilibrium is

$E_1=\left(\frac{r_1-c_3{N}_1-a_2{C}_1}{c_2}\mathrm{,0,}\frac{r_2-a_3\left(1-\eta \right)\frac{u}{d_2}}{r_2}\mathrm{,}\frac{u}{d_2}\right)$

$E_1$ exists if $r_1+u\left[c_{3}a-3\left(1-\eta \right)-a_2\right]/d_2>c_3{r}_2$ and $r_2>a_{3}u\left(1-\eta \right)/d_2$ .

5.2. Co-Axial Equilibrium Point

Next, we compute the Co-axial equilibrium point $E^{\mathrm{<mo>*</mo>}}=\left(I^{\mathrm{<mo>*</mo>}}\mathrm{,}{T}^{\mathrm{<mo>*</mo>}}\mathrm{,}{N}^{\mathrm{<mo>*</mo>}}\mathrm{,}{C}^{\mathrm{<mo>*</mo>}}\right)$ . From (4c), we get

$C^{\mathrm{<mo>*</mo>}}=\frac{u}{d_2+k{T}^{\mathrm{<mo>*</mo>}}}\mathrm{,}$ (6a)

then Equation (4d) becomes;

$r_2{N}^{\mathrm{<mo>*</mo>}}\left(1-N^{\mathrm{<mo>*</mo>}}\right)-c_{4}T{N}^{\mathrm{<mo>*</mo>}}-a_3\left(1-\eta \right)C^{\mathrm{<mo>*</mo>}}{N}^{\mathrm{<mo>*</mo>}}=\mathrm{0,}$

$N^{\mathrm{<mo>*</mo>}}=0\text{or}{r}_2\left(1-N^{\mathrm{<mo>*</mo>}}\right)-c_4{T}^{\mathrm{<mo>*</mo>}}-a_3\left(1-\eta \right)C^{\mathrm{<mo>*</mo>}}=0$

which gives;

$N^{\mathrm{<mo>*</mo>}}=\frac{r_2-c_4{T}^{\mathrm{<mo>*</mo>}}-a_3\left(1-\eta \right)C^{\mathrm{<mo>*</mo>}}}{r_2}\mathrm{,}$

by substituting Equation (6a);

$N^{\mathrm{<mo>*</mo>}}=\frac{r_2-c_4{T}^{\mathrm{<mo>*</mo>}}-a_3\left(1-\eta \right)\frac{u}{d_2+k{T}^{\mathrm{<mo>*</mo>}}}}{r_2}\mathrm{.}$ (6b)

From Equation (4b);

$r_1{T}^{\mathrm{<mo>*</mo>}}\left(1-b_1{T}^{\mathrm{<mo>*</mo>}}\right)-c_2{I}^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}-c_3{T}^{\mathrm{<mo>*</mo>}}{N}^{\mathrm{<mo>*</mo>}}-a_2{C}^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}=\mathrm{0,}$

$T^{\mathrm{<mo>*</mo>}}\left[r_1\left(1-b_1{T}^{\mathrm{<mo>*</mo>}}\right)-c_2{I}^{\mathrm{<mo>*</mo>}}-c_3{N}^{\mathrm{<mo>*</mo>}}-a_2{C}^{\mathrm{<mo>*</mo>}}\right]=\mathrm{0,}$

leads to;

$T^{\mathrm{<mo>*</mo>}}=0\text{or}{r}_1\left(1-b_1{T}^{\mathrm{<mo>*</mo>}}\right)-c_2{I}^{\mathrm{<mo>*</mo>}}-c_3{N}^{\mathrm{<mo>*</mo>}}-a_2{C}^{\mathrm{<mo>*</mo>}}=\mathrm{0,}$

this gives;

$I^{\mathrm{<mo>*</mo>}}=\frac{1}{c_2}\left(r_1\left(1-b_1{T}^{\mathrm{<mo>*</mo>}}\right)-c_3{N}^{\mathrm{<mo>*</mo>}}-a_2{C}^{\mathrm{<mo>*</mo>}}\right)\mathrm{,}$

and

$I^{\mathrm{<mo>*</mo>}}=\frac{1}{c_2}\left[r_1\left(1-b_1{T}^{\mathrm{<mo>*</mo>}}\right)-c_3\left[\frac{r_2-c_{4}T-a_3\left(1-\eta \right)\frac{u}{d_2+K{T}^{\mathrm{<mo>*</mo>}}}}{r_2}\right]-a_2\frac{u}{d_2+k{T}^{\mathrm{<mo>*</mo>}}}\right]\mathrm{,}$ (6c)

evaluating $T^{\mathrm{<mo>*</mo>}}$ ;

$\begin{array}{l}{c}_1{I}^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>2}}-\left[s+\rho I^{\mathrm{<mo>*</mo>}}-c_1\sigma I^{\mathrm{<mo>*</mo>}}-d_1{I}^{\mathrm{<mo>*</mo>}}-a_1\left(1-\eta \right)C^{\mathrm{<mo>*</mo>}}{I}^{\mathrm{<mo>*</mo>}}\right]T^{\mathrm{<mo>*</mo>}}\\ +\left[d_1\sigma I^{\mathrm{<mo>*</mo>}}+a_1\sigma \left(1-\eta \right)C^{\mathrm{<mo>*</mo>}}{I}^{\mathrm{<mo>*</mo>}}-s\sigma \right]=0.\end{array}$

After substituting the values of $I^{\mathrm{<mo>*</mo>}}$ and $C^{\mathrm{<mo>*</mo>}}$ , we get

$\frac{A_1{T}^{\mathrm{<mo>*</mo>5}}+A_2{T}^{\mathrm{<mo>*</mo>4}}+A_3{T}^{\mathrm{<mo>*</mo>3}}+A_4{T}^{\mathrm{<mo>*</mo>2}}+A_5{T}^{\mathrm{<mo>*</mo>}}+A_6}{A_7{T}^{\mathrm{<mo>*</mo>2}}+A_8{T}^{\mathrm{<mo>*</mo>}}+A_9}=0$ (6d)

where;

$A_1=c_1{c}_3{c}_4{k}^2-b_1{c}_1{k}^2{r}_1{r}_2$

$\begin{array}{c}{A}_2=c_3{c}_4{d}_1{k}^2-c_3{c}_4{k}^2\rho -c_1{c}_3{k}^2{r}_2+c_1{k}^2{r}_1{r}_2+2c_1{c}_3{c}_4{d}_{2}k+c_1{c}_3{c}_4{k}^2\sigma ,\\ -b_1{d}_1{k}^2{r}_1{r}_2+b_1{k}^2{r}_1{r}_2\rho -2b_1{c}_1{d}_{2}k{r}_1{r}_2-b_1{c}_1{k}^2{r}_1{r}_2\sigma ,\end{array}$

$\begin{array}{c}{A}_3=c_1{c}_3{c}_4{d}_2^2-c_3{d}_1{k}^2{r}_2+c_3{k}^2{r}_2\rho -c_2{k}^2{r}_{2}s+d_1{k}^2{r}_1{r}_2-k^2{r}_1{r}_2\rho \\ +2c_3{c}_4{d}_1{d}_{2}k+a_1{c}_3{c}_{4}ku-2c_3{c}_4{d}_{2}k\rho +c_3{c}_4{d}_1{k}^2\sigma +a_3{c}_1{c}_{3}ku\\ -2c_1{c}_3{d}_{2}k{r}_2-a_2{c}_{1}k{r}_{2}u+2c_1{d}_{2}k{r}_1{r}_2-b_1{c}_1{d}_2^2{r}_1{r}_2-c_1{c}_3{k}^2{r}_2\sigma \\ +c_1{k}^2{r}_1{r}_2\sigma +2c_1{c}_3{c}_4{d}_{2}k\sigma -a_1{c}_3{c}_4\eta ku-a_3{c}_1{c}_3\eta ku-2b_1{d}_1{d}_{2}k{r}_1{r}_2\\ -a_1{b}_{1}k{r}_1{r}_{2}u+2b_1{d}_{2}k{r}_1{r}_2\rho -b_1{d}_1{k}^2{r}_1{r}_2\sigma -2b_1{c}_1{d}_{2}k{r}_1{r}_2\sigma +a_1{b}_1\eta k{r}_1{r}_{2}u,\end{array}$

$\begin{array}{c}{A}_4=c_3{c}_4{d}_1{d}_2^2-c_3{c}_4{d}_2^2\rho -c_1{c}_3{d}_2^2{r}_2+c_1{d}_2^2{r}_1{r}_2+a_1{c}_3{c}_4{d}_{2}u+c_1{c}_3{c}_4{d}_2^2\sigma \\ +a_3{c}_1{c}_3{d}_{2}u+a_3{c}_3{d}_{1}ku-2c_3{d}_1{d}_{2}k{r}_2-a_2{c}_1{d}_2{r}_{2}u-a_1{c}_{3}k{r}_{2}u-a_3{c}_{3}k\rho u\\ +2c_3{d}_{2}k{r}_2\rho -a_2{d}_{1}k{r}_{2}u-2c_2{d}_{2}k{r}_{2}s+2d_1{d}_{2}k{r}_1{r}_2+a_{1}k{r}_1{r}_{2}u+a_{2}k{r}_2\rho u\\ -2d_{2}k{r}_1{r}_2\rho -b_1{d}_1{d}_2^2{r}_1{r}_2-c_3{d}_1{k}^2{r}_2\sigma +b_1{d}_2^2{r}_1{r}_2\rho -c_2{k}^2{r}_{2}s\sigma \\ +d_1{k}^2{r}_1{r}_2\sigma -a_1{c}_3{c}_4{d}_2\eta u+2c_3{c}_4{d}_1{d}_{2}k\sigma +a_1{c}_3{c}_{4}k\sigma u-a_1{c}_3{c}_4\eta k\sigma u\\ -a_3{c}_1{c}_3{d}_2\eta u-a_3{c}_3{d}_1\eta ku+a_3{c}_1{c}_{3}k\sigma u-2c_1{c}_3{d}_{2}k{r}_2\sigma +a_1{c}_3\eta k{r}_{2}u\\ +a_3{c}_3\eta k\rho u-a_1{b}_1{d}_2{r}_1{r}_{2}u-a_2{c}_{1}k{r}_2\sigma u+2c_1{d}_{2}k{r}_1{r}_2\sigma -a_1\eta k{r}_1{r}_{2}u\\ -b_1{c}_1{d}_2^2{r}_1{r}_2\sigma -a_1{b}_{1}k{r}_1{r}_2\sigma u-a_3{c}_1{c}_3\eta k\sigma u+a_1{b}_1{d}_2\eta r_1{r}_{2}u\\ -2b_1{d}_1{d}_{2}k{r}_1{r}_2\sigma +a_1{b}_1\eta k{r}_1{r}_2\sigma u,\end{array}$

$\begin{array}{c}{A}_5=c_3{d}_2^2{r}_2\rho -a_1{a}_2{r}_2{u}^2-c_3{d}_1{d}_2^2{r}_2-c_2{d}_2^2{r}_{2}s+d_1{d}_2^2{r}_1{r}_2-d_2^2{r}_1{r}_2\rho \\ +c_3{c}_4{d}_1{d}_2^2\sigma +a_1{a}_3{c}_3{u}^2+a_3{c}_3{d}_1{d}_{2}u-a_1{c}_3{d}_2{r}_{2}u-a_3{c}_3{d}_2\rho u\\ -a_2{d}_1{d}_2{r}_{2}u+a_1{d}_2{r}_1{r}_{2}u+a_2{d}_2{r}_2\rho u-2a_1{a}_3{c}_3\eta u^2+a_1{a}_2\eta r_2{u}^2\\ -c_1{c}_3{d}_2^2{r}_2\sigma +c_1{d}_2^2{r}_1{r}_2\sigma +a_1{a}_3{c}_3{\eta }^2{u}^2+a_1{c}_3{c}_4{d}_2\sigma u-a_1{c}_3{c}_4{d}_2\eta \sigma u\\ -a_3{c}_3{d}_1{d}_2\eta u+a_3{c}_1{c}_3{d}_2\sigma u+a_1{c}_3{d}_2\eta r_{2}u+a_3{c}_3{d}_2\eta \rho u+a_3{c}_3{d}_{1}k\sigma u\\ -2c_3{d}_1{d}_{2}k{r}_2\sigma -a_2{c}_1{d}_2{r}_2\sigma u-a_1{d}_2\eta r_1{r}_{2}u-a_1{c}_{3}k{r}_2\sigma u-a_2{d}_{1}k{r}_2\sigma u\\ -2c_2{d}_{2}k{r}_{2}s\sigma +2d_1{d}_{2}k{r}_1{r}_2\sigma +a_{1}k{r}_1{r}_2\sigma u-b_1{d}_1{d}_2{}^2{r}_1{r}_2\sigma -a_1\eta k{r}_1{r}_2\sigma u\\ -a_3{c}_1{c}_3{d}_2\eta \sigma u-a_3{c}_3{d}_1\eta k\sigma u+a_1{c}_3\eta k{r}_2\sigma u-a_1{b}_1{d}_2{r}_1{r}_2\sigma u\\ +a_1{b}_1{d}_2\eta r_1{r}_2\sigma u,\end{array}$

$\begin{array}{c}{A}_6=a_1{a}_3{c}_3\sigma u^2-c_3{d}_1{d}_2^2{r}_2\sigma -a_1{a}_2{r}_2\sigma u^2-c_2{d}_2^2{r}_{2}s\sigma +d_1{d}_2^2{r}_1{r}_2\sigma \\ +a_1{a}_3{c}_3{\eta }^2\sigma u^2+a_3{c}_3{d}_1{d}_2\sigma u-a_1{c}_3{d}_2{r}_2\sigma u-a_2{d}_1{d}_2{r}_2\sigma u\\ +a_1{d}_2{r}_1{r}_2\sigma u-2a_1{a}_3{c}_3\eta \sigma u^2+a_1{a}_2\eta r_2\sigma u^2-a_3{c}_3{d}_1{d}_2\eta \sigma u\\ +a_1{c}_3{d}_2\eta r_2\sigma u-a_1{d}_2\eta r_1{r}_2\sigma u,\end{array}$

$A_7=c_2{k}^2{r}_2,$

$A_8=2c_2{d}_{2}k{r}_2,$

$A_9=c_2{d}_2^2{r}_2,$

The co-axial equilibrium $E^{\mathrm{<mo>*</mo>}}$ exists if the roots of the Equation (6d) is positive, i.e., $T^{*}>0$ and the following inequalities must hold:

$r_1\left(1-b_1{T}^{*}\right)-\frac{c_3}{r_2}\left[r_2-c_{4}T-a_3\left(1-\eta \right)\frac{u}{d_2+K{T}^{*}}\right]-a_2\frac{u}{d_2+K{T}^{*}}>0,$

$\text{and}{r}_2>c_4{T}^{\mathrm{<mo>*</mo>}}+a_3\left(1-\eta \right)\frac{u}{d_2+K{T}^{\mathrm{<mo>*</mo>}}}\mathrm{.}$

To investigate the linear stability of the system around the two above stability, one must compute the jacobian of the system;

$J=\left[\begin{array}{cccc}{J}_{11}& J_{12}& 0& J_{13}\\ J_{21}& J_{22}& J_{23}& J_{24}\\ 0& J_{32}& J_{33}& J_{34}\\ 0& J_{42}& 0& J_{44}\end{array}\right]$ (7a)

where

$J_{11}=\rho T/\left(\sigma +T\right){\text{e}}^{-\lambda \tau }-c_{1}T-d_1-a_1\left(1-\eta \right)C$ ,

$J_{12}=\sigma \rho I/{\left(\sigma +T\right)}^2{\text{e}}^{-\lambda \tau }-c_{1}I$ ,

$J_{13}=-a_1\left(1-\eta \right)I$ ,

$J_{21}=-c_{2}T$ ,

$J_{22}=r_1-2r_1{b}_{1}T-c_{2}I-c_{3}N-a_{2}C$ ,

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

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

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

$J_{33}=r_2-2r_{2}N-c_{4}T-a_3\left(1-\eta \right)C$ ,

$J_{34}=-a_3\left(1-\eta \right)N$ ,

$J_{42}=-kC$ ,

$J_{44}=-d_2-kT$ .

The jacobian matrix of system (2a-2d) at the equilibrium point $E_1$ :

$\begin{array}{l}{J}_{E_1}\\ =\left[\begin{array}{cccc}-d_1-a_1\left(1-\eta \right)C_1& \frac{\rho I_1}{\sigma }{\text{e}}^{-\lambda \tau }-c_1{I}_1& 0& -a_1\left(1-\eta \right)I_1\\ 0& r_1-c_2{I}_1-c_3{N}_1-a_2{C}_1& 0& 0\\ 0& -c_4{N}_1& r_2-2r_2{N}_1-a_3\left(1-\eta \right)C_1& -a_3\left(1-\eta \right)N_1\\ 0& -k{C}_1& 0& -d_2\end{array}\right]\mathrm{,}\end{array}$

and so, the eigenvalues of the jacobian matrix corresponding to the steady-state $E_1$ are:

${\lambda }_1=-d_1-a_1\left(1-\eta \right)C_1$ ,

${\lambda }_2=r_1-c_2{I}_1-c_3{N}_1-a_2{C}_1$ ,

${\lambda }_3=r_2-2r_2{N}_1-a_3\left(1-\eta \right)C_1$ ,

and:

${\lambda }_4=-d_2$ .

Clearly, ${\lambda }_1<0$ and ${\lambda }_4<0$ .

Therefore, $E_1$ is stable if ${\lambda }_2<0$ i.e., if $r_1<c_2{I}_1+c_3{N}_1+a_2{C}_1$ and ${\lambda }_3<0$ i.e., if $r_2<2r_2{N}_1+a_3\left(1-\eta \right)C_1$ , otherwise $E_1$ is unstable.

The jacobian matrix of system (2a-2d) at the co-axial equilibrium point $E^{\mathrm{<mo>*</mo>}}$ is

$J_{E^{\mathrm{<mo>*</mo>}}}=\left[\begin{array}{cccc}{X}_2+X_1{\text{e}}^{-\lambda \tau }& X_5{\text{e}}^{-\lambda \tau }-I^{\mathrm{<mo>*</mo>}}{c}_1& 0& X_6\\ -T^{\mathrm{<mo>*</mo>}}{c}_2& X_3& -T^{\mathrm{<mo>*</mo>}}{c}_3& -T^{\mathrm{<mo>*</mo>}}{a}_2\\ 0& -N^{\mathrm{<mo>*</mo>}}{c}_4& X_4& X_7\\ 0& -C^{\mathrm{<mo>*</mo>}}k& 0& -d_2-T^{\mathrm{<mo>*</mo>}}k\end{array}\right]\mathrm{,}$

where

$X_1=\rho T^{\mathrm{<mo>*</mo>}}\left(\sigma +T^{\mathrm{<mo>*</mo>}}\right)$ ,

$X_2=-c_1{T}^{\mathrm{<mo>*</mo>}}-d_1-a_1\left(1-\eta \right)C^{\mathrm{<mo>*</mo>}}$ ,

$X_3=r_1-2r_1{b}_1{T}^{\mathrm{<mo>*</mo>}}-c_2{I}^{\mathrm{<mo>*</mo>}}-c_3{N}^{\mathrm{<mo>*</mo>}}-a_2{C}^{\mathrm{<mo>*</mo>}}$ ,

$X_4=r_2-2r_2{N}^{\mathrm{<mo>*</mo>}}-c_4{T}^{\mathrm{<mo>*</mo>}}-a_3\left(1-\eta \right)C^{\mathrm{<mo>*</mo>}}$ ,

$X_5=\sigma \rho I^{\mathrm{<mo>*</mo>}}{\left(\sigma +T^{\mathrm{<mo>*</mo>}}\right)}^2$ ,

$X_6=-a_1\left(1-\eta \right)I^{\mathrm{<mo>*</mo>}}$ ,

$X_7=-a_3\left(1-\eta \right)N^{\mathrm{<mo>*</mo>}}$ .

The eigenvalues associated with the coexisting equilibrium point $E^{\mathrm{<mo>*</mo>}}$ are derived from the characteristic equation

${\lambda }^4+X_{11}{\lambda }^3+X_{12}{\lambda }^2+X_{13}\lambda +X_{14}+\left(Y_{11}{\lambda }^3+Y_{12}{\lambda }^2+Y_{13}\lambda +Y_{14}\right){\text{e}}^{-\lambda \tau }=0$ (7b)

where

$X_{11}=d_2-X_3-X_4-X_2+T^{*}k,$

$\begin{array}{c}{X}_{12}=X_2{X}_3+X_2{X}_4+X_3{X}_4-X_2{d}_2-X_3{d}_2-X_4{d}_2-T^{*}{X}_{2}k-T^{*}{X}_{3}k\\ -T^{*}{X}_{4}k-C^{*}{T}^{*}{a}_{2}k-I^{*}{T}^{*}{c}_1{c}_2-N^{*}{T}^{*}{c}_3{c}_4,\end{array}$

$\begin{array}{c}{X}_{13}=X_2{X}_3{d}_2-X_2{X}_3{X}_4+X_2{X}_4{d}_2+X_3{X}_4{d}_2+T^{*}{X}_2{X}_{3}k+T^{*}{X}_2{X}_{4}k\\ +T^{*}{X}_3{X}_{4}k+C^{*}{T}^{*}{X}_2{a}_{2}k+C^{*}{T}^{*}{X}_4{a}_{2}k+I^{*}{T}^{*}{X}_4{c}_1{c}_2-C^{*}{T}^{*}{X}_6{c}_{2}k\\ -C^{*}{T}^{*}{X}_7{c}_{3}k+N^{*}{T}^{*}{X}_2{c}_3{c}_4-I^{*}{T}^{*}{c}_1{c}_2{d}_2-N^{*}{T}^{*}{c}_3{c}_4{d}_2-I^{*}{T}^{*}{}^2{c}_1{c}_{2}k\\ -N^{*}{T}^{*}{}^2{c}_3{c}_{4}k,\end{array}$

$\begin{array}{c}{X}_{14}=-X_2{X}_3{X}_4{d}_2-T^{*}{X}_2{X}_3{X}_{4}k-C^{*}{T}^{*}{X}_2{X}_4{a}_{2}k+C^{*}{T}^{*}{X}_2{X}_7{c}_{3}k\\ +C^{*}{T}^{*}{X}_4{X}_6{c}_{2}k+I^{*}{T}^{*}{X}_4{c}_1{c}_2{d}_2+N^{*}{T}^{*}{X}_2{c}_3{c}_4{d}_2+I^{*}{T}^{*}{}^2{X}_4{c}_1{c}_{2}k\\ +N^{*}{T}^{*}{}^2{X}_2{c}_3{c}_{4}k,\end{array}$

$Y_{11}=-X_1,$

$Y_{12}=X_1{X}_3+X_1{X}_4-X_1{d}_2+T^{*}{X}_5{c}_2-T^{*}{X}_{1}k,$

$\begin{array}{c}{Y}_{13}=X_1{X}_3{d}_2-X_1{X}_3{X}_4+X_1{X}_4{d}_2+T^{*}{}^2{X}_5{c}_{2}k-T^{*}{X}_4{X}_5{c}_2+T^{*}{X}_1{X}_{3}k\\ +T^{*}{X}_1{X}_{4}k+T^{*}{X}_5{c}_2{d}_2+C^{*}{T}^{*}{X}_1{a}_{2}k+N^{*}{T}^{*}{X}_1{c}_3{c}_4,\end{array}$

$\begin{array}{c}{Y}_{14}=-X_1{X}_3{X}_4{d}_2-T^{*}{X}_1{X}_3{X}_{4}k-T^{*}{X}_4{X}_5{c}_2{d}_2-T^{*}{}^2{X}_4{X}_5{c}_{2}k\\ -C^{*}{T}^{*}{X}_1{X}_4{a}_{2}k+C^{*}{T}^{*}{X}_1{X}_7{c}_{3}k+N^{*}{T}^{*}{X}_1{c}_3{c}_4{d}_2+N^{*}{T}^{*}{}^2{X}_1{c}_3{c}_{4}k,\end{array}$

The classical Routh-Hurwitz criterion does not apply to the delay system (1-4), since this equation is transcendental and has an infinite number of solutions. Substituting $\lambda =\varphi i$ ($\varphi >0$ ) into Equation (7b) gives;

$\begin{array}{l}{\left(\varphi i\right)}^4+X_{11}{\left(\varphi i\right)}^3+X_{12}{\left(\varphi i\right)}^2+X_{13}\left(\varphi i\right)+X_{14}\\ +\left(Y_{11}{\left(\varphi i\right)}^3+Y_{12}{\left(\varphi i\right)}^2+Y_{13}\left(\varphi i\right)\right){\text{e}}^{-\tau \varphi i}=\mathrm{0,}\end{array}$

by separating the imaginary and real parts leads to

$X_{11}{\varphi }^3-X_{13}\varphi =\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)\cos \left(\varphi \tau \right)+\left(Y_{12}{\varphi }^2-Y_{14}\right)\sin \left(\varphi \tau \right)\mathrm{,}$ (7c)

${\varphi }^4-X_{12}{\varphi }^2+X_{14}=\left(Y_{12}{\varphi }^2-Y_{14}\right)\cos \left(\varphi \tau \right)-\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)\sin \left(\varphi \tau \right)\mathrm{,}$ (7d)

squaring Equations (7c-7d) and adding the resulting equations yields to;

${\varphi }^8+p_{11}{\varphi }^6+p_{12}{\varphi }^5+p_{13}{\varphi }^4+p_{14}{\varphi }^2+p_{15}\varphi +p_{16}=\mathrm{0,}$ (7e)

where

$p_{11}=X_{11}^2-Y_{11}^2$ ,

$p_{12}=-2X_{12}$ ,

$p_{13}=2X_{14}-Y_{12}^2-2X_{11}{X}_{13}+2Y_{13}{Y}_{11}$ ,

$p_{14}=X_{13}^2-Y_{13}^2+2Y_{12}{Y}_{14}+X_{12}^2$ ,

$p_{15}=-2X_{12}{X}_{14}$ ,

$p_{16}=X_{14}^2-Y_{14}^2$ .

$p_{11}$ and $p_{16}$ can be re-written as;

$\begin{array}{c}{p}_{11}=(d_1+d_2-r_1-r_2+C^{\mathrm{<mo>*</mo>}}{a}_2+I^{\mathrm{<mo>*</mo>}}{c}_2+N^{\mathrm{<mo>*</mo>}}{c}_3+T^{\mathrm{<mo>*</mo>}}{c}_1+T^{\mathrm{<mo>*</mo>}}{c}_4+T^{\mathrm{<mo>*</mo>}}k\\ +{2N^{\mathrm{<mo>*</mo>}}{r}_2+2T^{\mathrm{<mo>*</mo>}}{b}_1{r}_1+C^{\mathrm{<mo>*</mo>}}{a}_1\left(1-\eta \right)+C^{\mathrm{<mo>*</mo>}}{a}_3\left(1-\eta \right))}^2-\frac{T^{\mathrm{<mo>*</mo>}}{}^2{\rho }^2}{{\left(\sigma +T^{\mathrm{<mo>*</mo>}}\right)}^2}\text{,}\end{array}$

and

$\begin{array}{c}{p}_{16}=-(\frac{T^{\mathrm{<mo>*</mo>}}{}^{2}k\rho {\sigma }_1{\sigma }_3}{T^{\mathrm{<mo>*</mo>}}+\sigma }+\frac{T^{\mathrm{<mo>*</mo>}}{d}_2\rho {\sigma }_1{\sigma }_3}{T^{\mathrm{<mo>*</mo>}}+\sigma }+\frac{N^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{}^2{c}_3{c}_4{d}_2\rho }{T^{\mathrm{<mo>*</mo>}}+\sigma }+\frac{N^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{}^3{c}_3{c}_{4}k\rho }{T^{\mathrm{<mo>*</mo>}}+\sigma }\\ -\frac{C^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{}^2{a}_{2}k\rho {\sigma }_1}{T^{\mathrm{<mo>*</mo>}}+\sigma }-\frac{I^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{}^2{c}_{2}k\rho \sigma {\sigma }_1}{{\left(T^{\mathrm{<mo>*</mo>}}+\sigma \right)}^2}-\frac{I^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{c}_2{d}_2\rho \sigma {\sigma }_1}{{\left(T^{\mathrm{<mo>*</mo>}}+\sigma \right)}^2}\\ +{\frac{C^{\mathrm{<mo>*</mo>}}{N}^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{}^2{a}_3{c}_{3}k\rho \left(\eta -1\right)}{T^{\mathrm{<mo>*</mo>}}+\sigma })}^2+(d_2{\sigma }_2{\sigma }_1{\sigma }_3+T^{\mathrm{<mo>*</mo>}}k{\sigma }_2{\sigma }_1{\sigma }_3\\ +N^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{}^2{c}_3{c}_{4}k{\sigma }_2-C^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{a}_{2}k{\sigma }_2{\sigma }_1-I^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{}^2{c}_1{c}_{2}k{\sigma }_1+N^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{c}_3{c}_4{d}_2{\sigma }_2\\ -I^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{c}_1{c}_2{d}_2{\sigma }_1-C^{\mathrm{<mo>*</mo>}}{I}^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{a}_1{c}_{2}k\left(\eta -1\right){\sigma }_1+{C^{\mathrm{<mo>*</mo>}}{N}^{\mathrm{<mo>*</mo>}}{T}^{\mathrm{<mo>*</mo>}}{a}_3{c}_{3}k\left(\eta -1\right){\sigma }_2)}^2\text{,}\end{array}$

where ${\sigma }_1=r_2-T^{\mathrm{<mo>*</mo>}}{c}_4-2N^{\mathrm{<mo>*</mo>}}{r}_2+C^{\mathrm{<mo>*</mo>}}{a}_3\left(\eta -1\right)$ , ${\sigma }_2=d_1+T^{\mathrm{<mo>*</mo>}}{c}_1-C^{\mathrm{<mo>*</mo>}}{a}_1\left(\eta -1\right)$ , and ${\sigma }_3=C^{\mathrm{<mo>*</mo>}}{a}_2-r_1+I^{\mathrm{<mo>*</mo>}}{c}_2+N^{\mathrm{<mo>*</mo>}}{c}_3+2T^{\mathrm{<mo>*</mo>}}{b}_1{r}_1$ . Equation (7e) will have a positive root if $P_{11}>0$ and $P_{16}<0$ . Now, we can conclude that Equation (7e) has at least one nonnegative real root and so, the characteristic Equation (7b) has purely imaginary roots $\pm \varphi i$ (say). This implies that there is a stability switch at $E^{\mathrm{<mo>*</mo>}}$ as $\tau$ changes. Eliminating $\sin \left(\varphi \tau \right)$ from (7c) and (7d) we get;

$\begin{array}{c}{\varphi }^4-X_{12}{\varphi }^2+X_{14}=\left(Y_{12}{\varphi }^2-Y_{14}\right)\cos \left(\varphi \tau \right)\\ -\frac{\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)\left(\left(X_{11}{\varphi }^3-X_{13}\varphi \right)-\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)\cos \left(\varphi \tau \right)\right)}{Y_{12}{\varphi }^2-Y_{14}}\text{,}\end{array}$

$\begin{array}{l}\left({\varphi }^4-X_{12}{\varphi }^2+X_{14}\right)\left(Y_{12}{\varphi }^2-Y_{14}\right)\\ ={\left(Y_{12}{\varphi }^2-Y_{14}\right)}^2\cos \left(\varphi \tau \right)-\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)\left(X_{11}{\varphi }^3-X_{13}\varphi \right)\\ -{\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)}^2\cos \left(\varphi \tau \right)\text{,}\end{array}$

$\begin{array}{l}\cos \left(\varphi \tau \right)=\frac{\left({\varphi }^4-X_{12}{\varphi }^2+X_{14}\right)\left(Y_{12}{\varphi }^2-Y_{14}\right)+\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)\left(X_{11}{\varphi }^3-X_{13}\varphi \right)}{{\left(Y_{12}{\varphi }^2-Y_{14}\right)}^2{\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)}^2}\text{,}\hfill \end{array}$

From the above equation, the time lag ${\tau }_n^{\mathrm{<mo>*</mo>}}$ corresponding to ${\varphi }_0$ is given by

$\begin{array}{c}{\tau }_n^{\mathrm{<mo>*</mo>}}=\frac{2n\pi }{{\varphi }_0}\\ +\frac{1}{{\varphi }_0}\arccos \left[\frac{\left({\varphi }^4-X_{12}{\varphi }^2+X_{14}\right)\left(Y_{12}{\varphi }^2-Y_{14}\right)+\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)\left(X_{11}{\varphi }^3-X_{13}\varphi \right)}{{\left(Y_{12}{\varphi }^2-Y_{14}\right)}^2{\left(Y_{13}\varphi -Y_{11}{\varphi }^3\right)}^2}\right]\text{,}\end{array}$

where $n$ is an integer. The equilibrium point $E^{\mathrm{<mo>*</mo>}}$ is locally asymptotically stable for all $\left[\mathrm{0,}{\tau }_0\right)$ where ${\tau }_0={\tau }_0^{\mathrm{<mo>*</mo>}}$ (by putting $n=0$ in the above expression of ${\tau }_0^{\mathrm{<mo>*</mo>}}$ ) if $p_{11}>0$ and $p_{16}<0$ [?].

6. Global Stability at Tumor-Free Steady State E₁

Let;

$V\left(\psi \right)={\psi }_1\left(0\right)+{\psi }_2\left(0\right)+{\psi }_3\left(0\right)+{\psi }_4\left(0\right)+{\displaystyle {\int }_{-\tau }^0}\frac{\rho {\psi }_1\left(s\right){\psi }_2\left(s\right)}{\sigma +{\psi }_2\left(s\right)}\text{d}s\mathrm{,}$

$\begin{array}{c}\frac{V\left(\psi \right)}{\text{d}t}=\frac{\text{d}{\psi }_1\left(0\right)}{\text{d}t}+\frac{{\psi }_2\left(0\right)}{\text{d}t}+\frac{{\psi }_3\left(0\right)}{\text{d}t}+\frac{{\psi }_4\left(0\right)}{\text{d}t}+\frac{\text{d}}{\text{d}t}{\displaystyle {\int }_{-\tau }^0}\frac{\rho {\psi }_1\left(s\right){\psi }_2\left(s\right)}{\sigma +{\psi }_2\left(s\right)}\text{d}s,\\ =s+\frac{\rho {\psi }_1\left(-\tau \right){\psi }_2\left(-\tau \right)}{\sigma +{\psi }_2\left(-\tau \right)}-c_1{\psi }_1\left(0\right){\psi }_2\left(0\right)-d_1{\psi }_1\left(0\right)-a_1\left(1-\eta \right){\psi }_4\left(0\right){\psi }_1\left(0\right)\\ +r_2{\psi }_2\left(0\right)\left(1-b_1{\psi }_2\left(0\right)\right)-c_2{\psi }_1\left(0\right){\psi }_2\left(0\right)-c_3{\psi }_2\left(0\right){\psi }_3\left(0\right)-a_2{\psi }_4\left(0\right){\psi }_2\left(0\right)\\ +r_2{\psi }_3\left(0\right)\left(1-{\psi }_3\left(0\right)\right)-c_4{\psi }_2\left(0\right){\psi }_3\left(0\right)-a_3\left(1-\eta \right){\psi }_4\left(0\right){\psi }_3\left(0\right)\\ +u-d_2{\psi }_4\left(0\right)-k{\psi }_2\left(0\right){\psi }_4\left(0\right)+\frac{\rho {\psi }_1\left(0\right){\psi }_2\left(0\right)}{\sigma +{\psi }_2\left(0\right)}-\frac{\rho {\psi }_1\left(-\tau \right){\psi }_2\left(-\tau \right)}{\sigma +{\psi }_2\left(-\tau \right)}\mathrm{,}\\ =s-d_1{\psi }_1\left(0\right)-a_1\left(1-\eta \right){\psi }_4\left(0\right){\psi }_1\left(0\right)+r_2{\psi }_3\left(0\right)\left(1-{\psi }_3\left(0\right)\right)\\ -a_3\left(1-\eta \right){\psi }_4\left(0\right){\psi }_3\left(0\right)+u-d_2{\psi }_4\left(0\right)\mathrm{.}\end{array}$

Using the equilibrium condition: $s=d_1{\psi }_1\left(0\right)-a_1\left(1-\eta \right){\psi }_4\left(0\right){\psi }_1\left(0\right)$ ,

$r_2{\psi }_3\left(0\right)\left(1-{\psi }_3\left(0\right)\right)-a_3\left(1-\eta \right){\psi }_4\left(0\right){\psi }_3\left(0\right)=0$ , $u=d_2{\psi }_4\left(0\right)$ , and the fact that

$\frac{V\left(\psi \right)}{\text{d}t}=0$ , we get the set

$E_1=\left\{\psi \in R_{+}^4:I\left(\xi \right)={\psi }_1\left(\xi \right),T\left(\xi \right)={\psi }_2\left(\xi \right)=0,N\left(\xi \right)={\psi }_3\left(\xi \right),C\left(\xi \right)={\psi }_4\left(\xi \right)\right\}.$

The classical La Salle’s invariance principle implies that $E_1$ is globally attractive. This confirms the global asymptotically stability of $E_1$ when$r_2<2r_2{N}_1+a_3\left(1-\eta \right)C_1$ .

7. Numerical Simulation

In this section, we used Matlab to graphically verify our analytical results for the system (2a-2d), which is crucial from a practical standpoint. The parameter values listed in Table 1 have been used in all simulations [19]. We assume that the parameter values’ units are arbitrary.

Table 1. Model’s parameters definitions and values.

Parameters	Definition	Values	Ref.
$s$	constant population of effector cells present in the body	0.05	[19]
$\rho$	maximum recruitment of effector cells by tumor cells	1	[19]
$\sigma$	half saturation constant for the proliferation term	0.4	[19]
$d_1$	effector cells’ natural death rate	0.2	[19]
$r_1$	intrinsic growth rate of tumor cells	0.4	[19]
$r_2$	normal cells’ growth rate	0.35	[19]
$\frac{1}{b_1}$	maximum carrying capacity of tumor cells	$\frac{1}{1.5}$	[19]
$d_2$	decay rate of targeted chemo-drug	0.05	[19]
$a_1$	kill rate of effector cell by targeted chemo-drug	0.2	[19]
$a_2$	kill rate of tumor cell by targeted chemo-drug	0.5	[19]
$a_3$	kill rate of normal cell by targeted chemo-drug	0.25	[19]
$c_1$	effector cells’ growth rate due to tumor cells	0.2	[19]
$c_2$	tumor cells’ decay rate due to immune cells	0.3	[19]
$c_3$	tumor cells’ decay rate due to normal cells	0.2	[19]
$c_4$	normal cells’ decay rate due to tumor cells	0.25	[19]
$\eta$	effectiveness of the targeted chemo-drug	0.01	[19]
$k$	rate of attachments of targeted chemo-drug with tumor cells	0.01	[19]

The two-dimensional (2D) plot of time $t$ versus the density of the normal cells $N\left(t\right)$ using the same data mentioned in Table 1 and different values of the delay factor $\tau =\mathrm{0.01,2,15}$ , is demonstrated in Figure 1. As time increases, the density decreases then increases with the highest value at $\tau =15$ at a larger value of time $t$ compared with $\tau =0.01$ . As time increases more, the density $N\left(t\right)$ decreases and tends to $E_1=\left(\mathrm{0.1766,0,0.7101,0.41}\right)$ . The case of higher delay $\tau =15$ takes a slightly longer time in approaching the steady state value with a larger number of normal cells compared with $\tau =0.01$ and 2.

Figure 1. Two dimensional plot (2D) of the density of tumor $T\left(t\right)$ versus time $t$ for the initial values $I\left(0\right)=0.6$ ; $T\left(0\right)=0.4$ ; $N\left(0\right)=0.9$ and $C\left(0\right)=0.1$ , $u=0.0205$ and the rest of data stated in Table 1 for the delay lag $\tau =\mathrm{0.01,2}$ and 15.

As time increases, Figure 2 shows that the density of the effector cells $I\left(t\right)$ is increased, then decreases and stabilizes at the steady-state value. It further shows that when the delay term increases it takes longer time to stabilize towards the equilibrium point $E_1$ , indicating the effectiveness of targeted treatment.

Figure 2. 2D plot shows the density Effector cells $I\left(t\right)$ versus time $t$ for same data of Figure 1.

For the same value of the chemotherapy-treatment, Figure 2 shows that the density of effector cells increases slowly as time increases (compared to normal cells). As the delay value increases, $\tau ={10}^{-6}\mathrm{,10,15}$ we conclude a higher increase of the effector cells’ density in a higher value of the time. A decrease on the density is observed at a higher value of time as the delay gets higher $\tau ={10}^{-6}\mathrm{,10,15}$ and stabilized at the tumor-free equilibrium point $E_1$ . This shows that when the delay term increases $I\left(t\right)$ takes more time to stabilize towards the equilibrium point $E_1$ .

The density of Tumor cells $T\left(t\right)$ , Figure 3, decreases rapidly as time increases reaching the stable state point at a shorter time for a larger delay value $\tau =15$ . This shows the effectiveness of the medication provided in an earlier stage. The earlier effector cells recognize the tumor cells, the better effectiveness of the medication afterwards. For this reason, and before any chemo-treatment begins, a less harmful medica-tion, probably natural remedies, has to be provided that will enhance the power and health of the effector cells.

Figure 3. 2D plot shows the density Tumor cells $T\left(t\right)$ versus time $t$ for same data of Figure 1.

For the chemotherapy treatment and as time increases, the density for different delay values stays close, Figure 4, as all stabilized towards the tumor-free point. This is logical as the treatment will be the same for any delay in $I$ and $T$ .

Figure 4. 2D plot shows the density chemotherapy cells $C\left(t\right)$ versus time $t$ for same data of Figure 1.

Figures 5(a)-(c) show the effect of increasing the chemotherapy treatment rate gradually on the different densities of the system for the delay term $\tau =15$ . As seen in Figures 5(a)-(c), the intersection points of the density-curves $I\left(t\right)$ , the effector, and $C\left(t\right)$ , the chemotherapy cells, happened in an earlier time of $t$ where the density of the effector cells become less than the chemotherapy ones. All densities tend to the free-tumor equilibrium point as time increases.

Figure 5. System’s densities when time delay $\tau =15$ and increasing treatment rate (a) $u=0.0205$ , (b) $u=0.021$ and (c) $u=0.024$ for with the same data of Figure 1.

Next, Figure 6 is showing the behaviour of the three densities of $I\left(t\right)$ , $T\left(t\right)$ and $C\left(t\right)$ for different delay factors.

Figure 6. 3D for the density of Effector, Normal and Tumor cells for the initial values $I\left(0\right)=0.6$ ; $T\left(0\right)=0.4$ ; $N\left(0\right)=0.9$ and $C\left(0\right)=0.1$ when (a) $\tau =0.01$ (b) $\tau =2$ and (c) $\tau =15$ with the same data of Table 1.

Figure 7(a), Figure 7(b) show the effect of increasing the intrinsic growth rate of the tumor cells $r_1=\mathrm{0.4,0.5}$ respectively, with the rest of data stay the same, $u=0.0205$ and the delay-time $\tau =15$ . For $r_1=0.4$ , Figure 7(a), we observe that as $C\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}$ increases, the density of the tumor cells decreases and tends to the free-tumor equilibrium value, the density of the effector cells increases up-to a certain value of $C\left(t\right)$ then rapidly decreases (as the number of tumor cells is reduced) and the density of the normal cells approach the equilibrium $E_1=\left[\mathrm{0.18,0,0.71,0.41}\right]$ . As $r_1=0.5$ increases, an almost similar behaviour was observed as $C\left(t\right)$ increases, but the system loses its stability as $C\left(t\right)$ approaches 0.41, Figure 7(b), Figure 7(c), moving away from the coexisting equilibrium point $E_2$ , while the effector system takes a longer time to respond appropriately to the tumor cells for recognition in the case of tumors with a higher growth rate (which is in the unstable range).

From a biological perspective, we know that the immune system does not stop responding to tumors until all tumor cells are removed; otherwise, it stays active. Immune cells need more time to fully destroy tumor cells in this procedure because, as time lag is increased.

In Figure 8, both trajectories of densities, $C\left(t\right)$ and $T\left(t\right)$ , are displayed as the density of the normal cells, $N\left(t\right)$ increases/decreases. Figure 8(a), the tumor cells-density increases slowly as $N\left(t\right)$ increases and $C$ decreases, dashed red-arrow, then a sudden-jump appeared (switch-up) at around $N=0.82$ to the upper branch for a larger value of the tumor-denisity and continue increases. The

(a)

(b)

(c)

Figure 7. System’s densities with time delay $\tau =15$ , $u=0.0205$ and intrinsic tumor growth rate (a) $r_1=0.4$ , (b) $r_1=0.5$ for with the same data of Figure 1.

(a)

(b)

Figure 8. Parametric plot of $N\left(t\right)$ versus $T\left(t\right)$ and $C\left(t\right)$ cells for the similar initial values as Figure 1, $\tau =15$ , $u=0.0205$ and $r_1=0.4$ .

switch-down effect happened as $N$ decreases, at around $N=0.79$ where the density of $T$ increases and switches-down to a very small value and continue to approach zero (the equilibrium point). The same discussion can be done for Figure 8(b) for $C$ . This bistability behaviours show the range where the system is unstable.

8. Conclusions

In this paper, we introduce a delay term into the interaction term between the immune system and tumor to study the dynamic behavior of the nonlinear model first proposed by Anusmita Das et al. [19]. This is to make the model more realistic. We have examined the fundamental properties, such as positivity and boundedness, of the solutions of the model. To investigate the model’s dynamic behavior, we conducted a stability analysis of the system in question. Our findings indicate that the tumor-free steady state is locally stable, given the following condition $r_1<c_2{I}_1+c_3{N}_1+a_2{C}_1$ . Furthermore, this steady state is globally stable. However, when there is a significant delay in the immune system’s recognition of tumor cells, resulting in a delayed response, the tumor growth rate accelerates. As a consequence, the system loses stability and moves away from the tumor-free equilibrium point $E_1$ .

In conclusion, our findings would indicate that there is potential for further valuable research in, for example, an examination of real-life data or the compar-ison of different treatment regimes, understanding and advancing the field, for example, a researcher might examine the results using real data or/and modify the system to suggest different delay type (un-equal or variable), or a combination of two types of treatments that will be sequentially applied.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	WHO: Cancer. https://www.who.int/health-topics/cancer#tab=tab1
[2]	Wang, C.H., Rockhill, J.K., Mrugala, M., Peacock, D.L., Lai, A., Jusenius, K., et al. (2009) Prognostic Significance of Growth Kinetics in Newly Diagnosed Glioblastomas Revealed by Combining Serial Imaging with a Novel Biomathematical Model. Cancer Research, 69, 9133-9140. https://doi.org/10.1158/0008-5472.can-08-3863
[3]	Macklin, P., Edgerton, M.E., Thompson, A.M. and Cristini, V. (2012) Patient-calibrated Agent-Based Modelling of Ductal Carcinoma in Situ (DCIS): From Microscopic Measurements to Macroscopic Predictions of Clinical Progression. Journal of Theoretical Biology, 301, 122-140. https://doi.org/10.1016/j.jtbi.2012.02.002
[4]	Gao, X., McDonald, J.T., Hlatky, L. and Enderling, H. (2012) Acute and Fractionated Irradiation Differentially Modulate Glioma Stem Cell Division Kinetics. Cancer Research, 73, 1481-1490. https://doi.org/10.1158/0008-5472.can-12-3429
[5]	Unni, P. and Seshaiyer, P. (2019) Mathematical Modeling, Analysis, and Simulation of Tumor Dynamics with Drug Interventions. Computational and Mathematical Methods in Medicine, 2019, Article ID: 4079298. https://doi.org/10.1155/2019/4079298
[6]	Din, Q. and Jameel, K. (2020) Mathematical Modelling of Tumor-Immune Interaction. Lap Lambert Academic Publishing.
[7]	Beigmohammadi, F., Khorrami, M., Masoudi, A.A. and Fatollahi, A.H. (2023) Mathematical Modeling of Tumor Growth as a Random Process in the Presence of Interaction between Tumor Cells and Normal Cells. International Journal of Modern Physics C, 35, Article ID: 2450102.
[8]	Galach, M. (2003) Dynamics of the Tumour-Immune System Competition: The Effect of Time Delay. International Journal of Applied Mathematics and Computer Science, 13, 395-406.
[9]	Mahlbacher, G.E., Reihmer, K.C. and Frieboes, H.B. (2019) Mathematical Modeling of Tumor-Immune Cell Interactions. Journal of Theoretical Biology, 469, 47-60. https://doi.org/10.1016/j.jtbi.2019.03.002
[10]	Allison, E., Colton, A.D., Gorman, A.D., Kurt, R. and Shainheit, M. (2004) A Mathematical Model of the Effector Cell Response to Cancer. Mathematical and Computer Modelling, 39, 1313-1327. https://doi.org/10.1016/j.mcm.2004.06.010
[11]	Liu, P. and Liu, X. (2017) Dynamics of a Tumor-Immune Model Considering Targeted Chemotherapy. Chaos, Solitons & Fractals, 98, 7-13. https://doi.org/10.1016/j.chaos.2017.03.002
[12]	Li, D., Ma, W. and Guo, S. (2016) Stability of a Mathematical Model of Tumour-Induced Angiogenesis. Nonlinear Analysis: Modelling and Control, 21, 325-344. https://doi.org/10.15388/na.2016.3.3
[13]	De Pillis, L.G. and Radunskaya, A. (2003) The Dynamics of an Optimally Controlled Tumor Model: A Case Study. Mathematical and Computer Modelling, 37, 1221-1244. https://doi.org/10.1016/s0895-7177(03)00133-x
[14]	De Pillis, L.G. and Radunskaya, A. (2000) A Mathematical Tumor Model with Immune Resistance and Drug Therapy: An Optimal Control Approach. Computational and Mathematical Methods in Medicine, 3, 79-100. https://doi.org/10.1080/10273660108833067
[15]	de Pillis, L.G., Gu, W., Fister, K.R., Head, T., Maples, K., Murugan, A., et al. (2007) Chemotherapy for Tumors: An Analysis of the Dynamics and a Study of Quadratic and Linear Optimal Controls. Mathematical Biosciences, 209, 292-315. https://doi.org/10.1016/j.mbs.2006.05.003
[16]	Arabameri, A., Asemani, D. and Hajati, J. (2018) Mathematical Modeling of In-Vivo Tumor-Immune Interactions for the Cancer Immunotherapy Using Matured Dendritic Cells. Journal of Biological Systems, 26, 167-188. https://doi.org/10.1142/s0218339018500080
[17]	Dehingia, K., Sarmah, H.K., Alharbi, Y. and Hosseini, K. (2021) Mathematical Analysis of a Cancer Model with Time-Delay in Tumor-Immune Interaction and Stimulation Processes. Advances in Difference Equations, 2021, Article No. 473. https://doi.org/10.1186/s13662-021-03621-4
[18]	Das, A., Dehingia, K., Sarmah, H.K., Hosseini, K., Sadri, K. and Salahshour, S. (2022) Analysis of a Delay-Induced Mathematical Model of Cancer. Advances in Continuous and Discrete Models, 2022, Article No. 15. https://doi.org/10.1186/s13662-022-03688-7
[19]	Das, A., Dehingia, K., Ray, N. and Sarmah, H.K. (2023) Stability Analysis of a Targeted Chemotherapy-Cancer Model. Mathematical Modelling and Control, 3, 116-126. https://doi.org/10.3934/mmc.2023011
[20]	Villasana, M. and Radunskaya, A. (2003) A Delay Differential Equation Model for Tumor Growth. Journal of Mathematical Biology, 47, 270-294. https://doi.org/10.1007/s00285-003-0211-0
[21]	Abernathy, Z., Abernathy, K. and Stevens, J. (2020) A Mathematical Model for Tumor Growth and Treatment Using Virotherapy. AIMS Mathematics, 5, 4136-4150. https://doi.org/10.3934/math.2020265
[22]	Yang, X., Chen, L. and Chen, J. (1996) Permanence and Positive Periodic Solution for the Single-Species Nonautonomous Delay Diffusive Models. Computers & Mathematics with Applications, 32, 109-116. https://doi.org/10.1016/0898-1221(96)00129-0