Tumor-Immune Interaction System with the Effect of Time Delay and Hyperglycemia on the Breast Cancer Cells

Abeer Hamdan Alblowy^1,2, Normah Maan^2*, Nor Aziran Awang²

¹Department of Mathematics, Faculty of Science, University of Ha’il, Ha’il 2440, Saudi Arabia.
²Department of Mathematical Sciences, Faculty of Science, Universiti Teknologi Malaysia, Johor Bahru 81310, Malaysia.

DOI: 10.4236/jamp.2023.114076                

Abstract

Breast cancer in women is a complicated and multifaceted disease. Studies have demonstrated that hyperglycemia is one of the most significant risk factors for breast cancer. Hyperglycemia is when the sugar level in human blood is too high, which means excess glucose. Glucose excess can encourage the growth, invasion, and migration of breast cancer cells at the cellular level. Though, the effects of glucose on the dynamics of breast cancer cells have been examined mathematically by a system of ordinary differential equations. However, the non-instantaneous biological occurrences leading to the secretion of immuno-suppressive cytokines by tumors to evade immune surveillance and the immune cells’ derivation of cytokines to attack the tumor cells are not yet discussed. Therefore, investigating the biological process involved in the dynamics of tumors, immune and normal cells with excessive glucose concentration is inviolable to determining the best procedure for controlling tumors’ uncontrollable growth. Time delay, denoted by τ, is used to describe the time tumor cells take to secrete immunosuppressive cytokines to evade immune surveillance and the time immune cells take to recognize and attack the tumor cells. We have studied the local stability analysis of the biological steady states in both delayed and non-delayed system. The Routh-Hurwitz stability criterion is used to analyze the dynamical equilibrium of the cells’ population. Hopf bifurcation was analyzed by using time delay s as a bifurcation parameter. The analytical results suggest an unstable scenario for a tumor-free equilibrium point as normal cells are bound to grow to their carrying capacity. The result predicts a stable system for coexisting equilibrium when the interaction is instantaneous (τ = 0). However, when τ > 0, the coexisting equilibrium point switches from stable to unstable. The numerical results not only validate all the analytical results but also show the case of possible situations when glucose concentration is varied, indicating that both tumor growth and immune system efficiency are highly affected by the level of glucose in the blood. This concluded that the delay in the secretion of cytokines by immune cells and derivation cytokines by the tumors helps to identify the possible chaotic situation under different glucose concentration and the extent to which such delay can have on restoration of the normal cells when glucose concentration is low.

Keywords

Breast Cancer, Stability, Bifurcation, Time Delay, Nonlinear Dynamical Model, Glucose Risk Factor

1. Introduction

A considerable number of studies stated that breast cancer (BC) is the second type of cancer which is common in women and is considered the second most common factor in cancer deaths among females of all ages [1] . Most cancers occur due to the DNA structure changes inside the cell, called the mutations in DNA [2] . Breast cancer is one of the most common types of malignant tumours that affect breast tissues, and lead to abnormal growth in its cells [3] . It begins to spread from the inner lining, and then move from it to the milk ducts and the lobules [4] . Breast cancer is distinguished by its ability to spread to other areas of the body. Some factors can cause breast cancer, such as estrogen excess, genetic factors, and some environmental factors or life style such as poor diet, excessive alcohol consumption and smoking [4] [5] . On the other hand, there is some factors enhance cancer cells to grow and spread, such as increase of glucose in the blood which can be called hyperglycemia.

Hyperglycemia is one of the risk factors which increase the population of cancer cells rapidly including breast cancer cells as stated in many recent biological studies [6] - [14] . Hyperglycemia means the high glucose level in the blood. At the cellular level, hyperglycemia (glucose excess) can encourage the growth, invasion, and migration of breast cancer cells. It can also promote anti-apoptotic reactions to boost tumors’ chemo-resistance by causing abnormal glucose metabolism [15] [16] . The immune system is also affected by hyperglycemia. High glucose means low immune system efficiency [17] [18] . The evidence is that people with diabetes are more prone to infections and delayed wound healing, as high blood sugar levels weaken the immune system. High glucose decreases the immune system response against the aggression and diseases by reducing the efficiency of some complement factors and polynuclear cells chemotactism and phagocytosis and by increasing the inflammatory response [19] [20] .

The dynamic of breast cancer cells has been described mathematically in many researches. Mathematical models assist to comprehending the dynamics of tumour growth over time [21] . Two of these models are ordinary differential equations (ODEs) and delay differential equations (DDEs). ODEs are considered a very appropriate tool for modeling cancer problem. ODEs can be used to recognize the effect of chemotherapy treatment on breast cancer cells, then suggest optimal medicine [22] . ODEs can be improved to the delay differential equations (DDEs) if there is a delay in completing any process. For example, Intra and intercellular reaction in human body undergo biological process to manifest, and this suggest non instantaneous reaction which made using DDEs more suitable than ODEs [23] . Besides, the DDEs differ from ODESs in that DDEs have derivatives depend on solving the problems in previous time while ODEs contain derivatives depend on the solution at present time [24] .

Many mathematical models have been developed to understand the behavior of breast cancer cells (BCCs) under the effect of time delay. For instance, in 2020 Khajanchi [25] improved De Pillis and Radunskaya model [26] by introducing the time delay to the interaction between tumour and immune cells in order to get a better compatibility with reality. There are clinical evidences that most anti-tumor’s therapeutics also undergo process of healing which varies from 2 - 10 weeks after initiation of treatment [27] . In another research, Udomchalermpat et al. [28] extended the tumor-obesity model which was studied in [29] by including the time delay effect into immune cells equation. They incorporated into their model time delay to model the non-instantaneous interaction between tumors and immune cells occasioned by the process involves in developing suitable response by immune system and modulation of immune cells’ attack by the tumor. Likewise, Rihan et al. [21] reduced Kuznetsov et al. [30] model from five equations to be only two equations which are tumor and immune cells equations. The simplified model explains the interactions between tumor cells and immune cells such as cytotoxic T-cells or natural killer cells. In order to improve the results, time delay was added to immune cells equation especially into the tumor-immune interaction term as well as the nonlinear growth term for the immune cells. In 2022, Awang et al. [31] illustrated the tumor-immune interaction using a mathematical model by considering the cell cycle. The delay was included to tumor cell reside inter-phase. Ibrahim et al. [32] used a discrete time delay to model treatments’ control against inhibitory agents and suppression in addition to modification of immune cells by the tumor. Their findings suggested that control is not enough to eliminate malignant tumors. All the above previous studies contributed to give more accurate results compared with the studies that used ODEs since they simulate the real situation of cells interactions which produce delay in the time. The reason behind including the time delay in the previous researches was that the mathematical models describes population dynamics which are cells dynamics, and this kind of study contains a certain hidden process that required time to be completed, such as cell cycle stages, the immune period, the time of healthy cells and cancer cells interactions [23] [33] [34] . Besides, the malignant transformation is not an instantaneous process, it is logical to consider how a time lag can affect the dynamics brought on by the interactions between these three-cell populations. [35] . As a result, cells interactions can be modeled by using DDEs.

Alblowy et al. [36] studied the dynamics of tumor, immune and normal cells with incorporation of hyperglycemia (excessive glucose). The interesting thing in their study is that the glucose risk factor has not been studied in mathematical cancer models before. Therefore, they were the first authors who investigated the effect of glucose excess on the cells dynamics mathematically. Their results showed that glucose excess caused instability and chaos in immune and tumor cells. High level of glucose caused an increase in the growth of cancer cells and a decrease in immune system efficiency, causing chaos in the cells dynamics. However, the secretion of cytokines by the effector cells and derivation of suppressive cytokines are instantaneous against the biological hypothesis that it undergoes biological process to manifest [35] . It is significant to consider how a delay may affect the dynamics resulting from the interaction between tumour and immune cells as the interaction of various types of cells is not an instantaneous process. Hence, the model with time delay will give better representation and closer to the real situation. Also, when the time delay is introduced to any system of ODEs, the ultimate motive behind that is to determine whether Hopf bifurcation occurs as time delay increases such that the stability switches from stable to unstable or from unstable to stable. This indicates that the system’s eigenvalues change from having negative real parts to having positive real parts which happens only when they cross the imaginary axis [37] . In this work, time delay is incorporated into the system proposed in [36] to investigate the biological process leading to suppression of immune cells and the derivation of cytokines by the tumor cells. In Section 2, the modified system of the system proposed in [36] with time delay is presented. Qualitative analysis of the model including positivity theorem, existence of the equilibrium points, stability analysis and occurrence of Hopf bifurcation are provided in Section 3. Section 4 contains the numerical simulations of the model to validate the analytical results. Section 5 briefly discusses the main findings of this study.

2. The Breast Cancer Model

Cancer cells live in an environmental system that forces them to interact with immune cells. The interaction between these cells is like the prey-predator relationship [38] . Immune cells, like natural killer, macrophages, and CD8+T cells, are the important elements in the immune system that destroy the tumor cells through a kinetic process by which the tumor cells come in contact with immune cells and making them functionally inactive. Cancer cells secrete immunosuppressive cytokines to evade immune surveillance [25] [33] . Thus, instantaneously, the immune cells are unable to destroy the tumor cells so there is a time delay between the deactivation of tumor cells by immune cells. This time delay can be regarded as an interaction delay. As a result, we propose to modify the model in [36] which studied the effect of glucose risk factor on a breast cancer model since they considered the interplay between cancer and immune cells as instantaneous but contrarily, it undergoes some biological process before such interaction can manifest. Therefore, in this study, a time delay is incorporated into the model in [36] to depict the biological process leading to delay occasioned by the interaction in the suppression of immune cells and the derivation of cytokines by the tumor cells as well as to investigate how the resulting dynamics are affected by the time delay. The proposed model which includes three equations: normal cells, $N\left(t\right)$ , tumor cells, $T\left(t\right)$ , and immune cells, $M\left(t\right)$ , is presented by delay differential system as given below;

$\begin{array}{l}\frac{\text{d}N}{\text{d}t}=N\left(t\right)\left({\alpha }_1-{\mu }_{1}N\left(t\right)\right)-{\varphi }_{1}N\left(t\right)T\left(t\right)\\ \frac{\text{d}T}{\text{d}t}=T\left(t\right)\left({\alpha }_2-{\mu }_{2}T\left(t\right)\right)+gT\left(t\right)-{\gamma }_{1}M\left(t-\tau \right)T\left(t-\tau \right)+{\varphi }_{2}N\left(t\right)T\left(t\right)\\ \frac{\text{d}M}{\text{d}t}=s+\frac{\rho M\left(t\right)T\left(t\right)}{\omega +T\left(t\right)}-{\gamma }_{2}M\left(t-\tau \right)T\left(t-\tau \right)-{\mu }_{3}M\left(t\right)-gM\left(t\right)\end{array}$ (1)

with the initial conditions

$N_0={\varphi }_1\left(\xi \right),T_0={\varphi }_2\left(\xi \right),M_0={\varphi }_3\left(\xi \right),\text{for}\xi \in \left[-\tau ,0\right].$

such that

$\tau$ is used to describe the time tumor cells take to secrete immunosuppressive cytokines to evade immune surveillance and the time immune cells take to recognize and attack the tumor cells. ${\gamma }_{1}M\left(t-\tau \right)T\left(t-\tau \right)$ represents the delay in the clearance of cancer cells by immune cells, and ${\gamma }_{2}M\left(t-\tau \right)T\left(t-\tau \right)$ describes the delay in the clearance of immune cells by tumour cells. Other parameters in the model (1) are positive and defined as follows:

• ${\alpha }_1$ is the growth rate of the normal breast cells.

• ${\mu }_1$ is the death rate of the normal cells.

• ${\alpha }_2$ is the growth rate of the tumor cells.

• ${\mu }_2$ is the inhibition rate of the breast tumor cells.

• ${\varphi }_1$ is the rate at which tumor cells impair the normal cells.

• ${\varphi }_2$ denotes the rate at which tumor cells infect the normal cells to expand.

• $g$ represents the rate of glucose excess.

• $s$ is the source rate of immune cells.

• ${\mu }_3$ is the natural death rate of the immune cells.

• $\rho$ is the immune response rate.

• $\omega$ is the immune threshold rate.

• ${\gamma }_1$ is the rate at which the tumor cells are reduced by the immune cells.

• ${\gamma }_2$ is the rate at which the immune cells are reduced by the tumor cells.

Figure 1 shows a schematic diagram for the cells’ population interaction in the proposed model. $\tau$ describes the interaction delay between tumor and immune cells.

Figure 1. Schematic diagram for the cells’ population competition in the proposed model, $\tau$ describes the interaction delay between tumor and immune cells.

3. Qualitative Analysis

3.1. Existence of Non-Negative Solutions

For the need to demonstrate the biological compliance of System (1), it is sacrosanct to show that the solution of the system are bounded and non-negative. Hence, the existence of non-negative solution System (1) is given in Theorem 1 below:

Theorem 1 Suppose $N\left(t\right)$ , $T\left(t\right)$ and $M\left(t\right)$ are bounded and positive under the initial conditions in System (1), then, there exist non-negative solutions $N\left(t\right)$ , $T\left(t\right)$ and $M\left(t\right)$ of System (1) $\forall t\in \left[t-\tau \mathrm{,}t\right]$ .

Proof. Solving for $N\left(t\right)\mathrm{,}T\left(t\right)$ and $M\left(t\right)$ in System (1), we have.

Firstly, for $N(\; t\; )$

$\frac{\text{d}N}{\text{d}t}+\left({\varphi }_{1}T\left(t\right)-{\alpha }_1\right)N\left(t\right)=-{\mu }_1{N}^2\left(t\right).$ (2)

Equation (2) is obviously Bernouli equation, solving this equation we obtain the solution of $N\left(t\right)$ as follows

$N\left(t\right)=\frac{N\left(0\right){\text{e}}^{-{\displaystyle \int }\left({\varphi }_{1}T\left(t\right)-{\alpha }_1\right)\text{d}t}}{{\text{e}}^{-{\displaystyle \int }\left({\varphi }_{1}T\left(t\right)-{\alpha }_1\right)\text{d}t}+N\left(0\right){\displaystyle {\int }_0^t}{\mu }_1{\text{e}}^{-{\displaystyle \int }\left({\varphi }_{1}T\left(t\right)-{\alpha }_1\right)\text{d}t}\text{d}t}.$ (3)

Applying the same method to find $T\left(t\right)$ and $M\left(t\right)$ in System (1) yields the following solutions

$T\left(t\right)=\frac{T\left(0\right){\text{e}}^{{\displaystyle \int }\left[{\alpha }_2+{\varphi }_{2}N\left(t\right)+g\right]\text{d}t}}{{\text{e}}^{{\displaystyle \int }\left[{\alpha }_2+{\varphi }_{2}N\left(t\right)+g\right]\text{d}t}+T\left(0\right)\left[{\displaystyle {\int }_0^t}{\mu }_2{\text{e}}^{{\displaystyle \int }\left[{\alpha }_2+{\varphi }_{2}N\left(t\right)+g\right]\text{d}t}\text{d}t+{\displaystyle {\int }_{-\tau }^0}{\gamma }_{1}M\left(t-\tau \right)T\left(t-\tau \right){\left[T\left(t\right)\right]}^{-2}{\text{e}}^{{\displaystyle \int }\left[{\alpha }_2+{\varphi }_{2}N\left(t\right)+g\right]\text{d}t}\text{d}t\right]}$ (4)

$M\left(t\right)=\frac{M\left(0\right){\text{e}}^{-{\displaystyle \int }\left[\frac{\rho T\left(t\right)}{\omega +T\left(t\right)}-{\mu }_3-g\right]\text{d}t}+{\displaystyle {\int }_0^t}s\text{d}t-{\displaystyle {\int }_{-\tau }^0}{\gamma }_{2}M\left(t-\tau \right)T\left(t-\tau \right)\text{d}t}{{\text{e}}^{-{\displaystyle \int }\left[\frac{\rho T\left(t\right)}{\omega +T\left(t\right)}-{\mu }_3-g\right]\text{d}t}}.$ (5)

Obviously, if $N\left(0\right)>0$ , $T\left(0\right)>0$ and $M\left(0\right)>0$ , System (1) has a non-negative solution. $⊡$

3.2. The Existence of Equilibrium Points

Theorem 2 Assume there exist positively invariant points in system (1), then there exist two biologically relevant equilibrium points, namely tumor-free equilibrium and coexistence equilibrium.

Proof. To find the equilibrium points of system (1) set L.H.S of the system equal to zero as follows:

$\begin{array}{l}N\left({\alpha }_1-{\mu }_{1}N\right)-{\varphi }_{1}NT=0\\ T\left({\alpha }_2-{\mu }_{2}T\right)+gT-{\gamma }_{1}MT+{\varphi }_{2}NT=0\\ s+\frac{\rho MT}{\omega +T}-{\gamma }_{2}MT-{\mu }_{3}M-gM=0\end{array}$ (6)

solving System (6) yields the following two biologically feasible equilibrium points:

1) Tumor-Free Equilibrium:

$E_0=\left(N^{\ast },T^{\ast },M^{\ast }\right)=\left(\frac{{\alpha }_1}{{\mu }_1},0,\frac{s}{g+{\mu }_3}\right).$ (7)

This state represents the health equilibrium. Only normal cells and immune cells are alive.

It is clear that $N^{\ast }\mathrm{,}{T}^{\ast }$ and $M^{\ast }$ are positive since the parameters ${\alpha }_1\mathrm{,}{\mu }_1\mathrm{,}s\mathrm{,}g$ and ${\mu }_3$ are positive real values then $E_0$ lies in the feasible region.

2) Coexisting Equilibrium:

$E_1=\left(N^{\ast },T^{\ast },M^{\ast }\right)=\left(\frac{1}{{\mu }_1}\left({\alpha }_1-{\varphi }_1{T}^{\ast }\right),T^{\ast },\frac{1}{{\gamma }_1{\gamma }_2{\mu }_1}\left({\gamma }_2{l}_2-k_1{T}^{\ast }\right)\right)$ (8)

$T^{\ast }$ can be found by solving the equation:

$k_1{\left(T^{\ast }\right)}^3+k_2{\left(T^{\ast }\right)}^2+k_3{T}^{\ast }+k_4=0,$ (9)

$k_1={\gamma }_2\left({\mu }_1{\mu }_2+{\varphi }_1{\varphi }_2\right)>0,$

$k_2=\frac{k_1}{{\gamma }_2}{l}_1-{\gamma }_2{l}_2,$

$k_3=\frac{\omega }{{\gamma }_2}{k}_1\left(g+{\mu }_3\right)-l_1{l}_2+s{\gamma }_1{\mu }_1,$

$k_4=-\omega \left[\left(g+{\mu }_3\right)l_2-s{\gamma }_1{\mu }_1\right],$ (10)

and

$l_1=\omega {\gamma }_2+g+{\mu }_3-\rho ,$

$l_2={\mu }_1\left(g+{\alpha }_2\right)+{\alpha }_1{\varphi }_2.$

From Equation (8), we conclude that $E_1$ can exist under the condition:

$0<T^{\ast }<\min \left\{\frac{{\alpha }_1}{{\varphi }_1},\frac{{\gamma }_2{l}_2}{k_1}\right\}.$ (11)

This completes the proof.

3.3. Local Stability of Equilibrium Points

Here, we seek to determine the long-time behaviors of the interacting cells around the equilibrium points. This is done by linearizing System (1) at the equilibrium points. The Jacobian matrix of System (1) at the equilibrium point $E^{\ast }=\left(N^{\ast },T^{\ast },M^{\ast }\right)$ becomes

$J\left(E^{\ast }\right)=\left[\begin{array}{ccc}{\alpha }_1-2{\mu }_1{N}^{\ast }-{\varphi }_1{T}^{\ast }& -{\varphi }_1{N}^{\ast }& 0\\ {\varphi }_2{T}^{\ast }& -{\gamma }_1{M}^{\ast }{\text{e}}^{-\lambda \tau }+{\varphi }_2{N}^{\ast }-2{\mu }_2{T}^{\ast }+g+{\alpha }_2& -{\gamma }_1{T}^{\ast }{\text{e}}^{-\lambda \tau }\\ 0& M^{\ast }\left(\frac{\rho \omega }{{\left(\omega +T^{\ast }\right)}^2}-{\gamma }_2{\text{e}}^{-\lambda \tau }\right)& T^{\ast }\left(\frac{\rho }{\omega +T^{\ast }}-{\gamma }_2{\text{e}}^{-\lambda \tau }\right)-\left({\mu }_3+g\right)\end{array}\right].$ (12)

For the time delay $\tau$ , the characteristic polynomial around $E^{\ast }=\left(N^{\ast },T^{\ast },M^{\ast }\right)$ will be on the form:

$P\left(\lambda \mathrm{,}\tau \right)=A\left(\lambda \right)+B\left(\lambda \right){\text{e}}^{-\lambda \tau }$

Hence, the stability conditions for both the tumor-free and co-existence equilibrium points are given by Theorem (3 & 4) below:

Theorem 3 System (1) is locally stable at $E_0$ , $\forall \tau \ge 0$ , if $\frac{{\gamma }_{1}s}{{\mu }_3+g}>\frac{{\phi }_2{\alpha }_1}{{\mu }_1}+g+{\alpha }_2$ , otherwise it is unstable.

Proof. Substituting $\left(N^{\ast },T^{\ast },M^{\ast }\right)=\left(\frac{{\alpha }_1}{{\mu }_1},0,\frac{s}{g+{\mu }_3}\right)$ , Equation (12) becomes

$J\left(E_0\right)=\left[\begin{array}{ccc}-{\alpha }_1& \frac{-{\alpha }_1{\varphi }_1}{{\mu }_1}& 0\\ 0& -\frac{{\gamma }_{1}s}{{\mu }_3+g}{\text{e}}^{-\lambda \tau }+\frac{{\varphi }_2{\alpha }_1}{{\mu }_1}+g+{\alpha }_2& 0\\ 0& \frac{s}{{\mu }_3+g}\left(\frac{\rho }{\omega }-{\gamma }_2{\text{e}}^{-\lambda \tau }\right)& -\left({\mu }_3+g\right)\end{array}\right]$ (13)

then the characteristic equation of $E_0$ becomes

$\left(-{\alpha }_1-\lambda \right)\left(\frac{-{\gamma }_{1}s}{g+{\mu }_3}{\text{e}}^{-\lambda \tau }+\frac{{\alpha }_1{\phi }_2}{{\mu }_1}+g+{\alpha }_2-\lambda \right)\left(-\left({\mu }_3+g\right)-\lambda \right)=0.$ (14)

• When $\tau =0$ , obviously

${\lambda }_1=-{\alpha }_1$

${\lambda }_2=-\left({\mu }_3+g\right)$

and the third root of Equation (14) is

${\lambda }_3=\frac{-{\gamma }_{1}s}{g+{\mu }_3}+\frac{{\alpha }_1{\phi }_2}{{\mu }_1}+g+{\alpha }_2$

we conclude that ${\lambda }_3$ is negative if and only if the following condition is satisfied

$\frac{{\gamma }_{1}s}{{\mu }_3+g}>\frac{{\phi }_2{\alpha }_1}{{\mu }_1}+g+{\alpha }_2,$ (15)

Biologically, $\frac{{\alpha }_1{\phi }_2}{{\mu }_1}+g+{\alpha }_2$ depicts the rate at which normal cell growth to capacity in the present of tumor’s proliferation and glucose concentration, while $\frac{{\gamma }_{1}s}{g+{\mu }_3}$ describes the rate at which immune cells annex their source to reach carrying capacity in the midst of glucose concentration. It is biologically meaningful to assume that in the absence of the tumor cells normal cells will attain their carry capacity and there will be relatively no response to tumors’ antigen. Thus, we assume $\frac{{\gamma }_{1}s}{g+{\mu }_3}<\frac{{\phi }_2{\alpha }_1}{{\mu }_1}+g+{\alpha }_2$ and opine that $E_0$ is locally unstable for $\tau =0$ .

• When $\tau >0$ : the stability of $E_0$ is determined by the negativity of the root;

$\frac{-{\gamma }_{1}s}{g+{\mu }_3}{\text{e}}^{-\lambda \tau }+\frac{{\alpha }_1{\phi }_2}{{\mu }_1}+g+{\alpha }_2-\lambda =0$ (16)

of Equation (14).

Hence, Equation (16) takes the form

$\lambda -\left(\frac{{\alpha }_1{\phi }_2}{{\mu }_1}+g+{\alpha }_2\right)+\frac{{\gamma }_{1}s}{g+{\mu }_3}{\text{e}}^{-\lambda \tau }=0$ (17)

for the simplicity, Equation (17) can be written as

$\lambda +D+E{\text{e}}^{-\lambda \tau }=0.$ (18)

where

$D=-\left(\frac{{\alpha }_1{\phi }_2}{{\mu }_1}+g+{\alpha }_2\right)\mathrm{.}$

$E=\frac{{\gamma }_{1}s}{g+{\mu }_3}.$

Defining $\lambda =a+i\theta$ where $\theta >0$ and $a=0$ . Substituting by the value of in $\lambda =a+i\theta$ in Equation (18). We obtain

$i\theta +D+E{\text{e}}^{-\left(i\theta \right)\tau }=0$ (19)

Rewriting the exponential function in terms of trigonometric functions, we obtain

$i\theta +D+E\left(\cos \left(\theta \tau \right)-i\sin \left(\theta \tau \right)\right)=0$ (20)

Separating real and imaginary parts of Equation (20) yields

$D=-E\cos \left(\theta \tau \right).$ (21)

$\theta =E\sin \left(\theta \tau \right).$ (22)

Squaring both sides of Equations (21) and (22) then adding them together

${\theta }^2+D^2=E^2\left[{\cos }^2\left(\theta \tau \right)+{\sin }^2\left(\theta \tau \right)\right].$ (23)

Hence,

${\theta }^2=E^2-D^2.$ (24)

Obviously from (24), ${\theta }^2<0$ , if $E^2<D^2$ . Since $E^2<D^2$ replicates the instability condition when $\tau =0$ , Equation (14) has no positive root when $\tau >0$ . In other words, there is no stability switch, and there exists no occurrence of Hopf bifurcation. Hence, System (1) is unstable at $E_0\mathrm{,}\forall \tau \ge 0$ , if, and only if,

$\frac{{\gamma }_{1}s}{{\mu }_3+g}<\frac{{\phi }_2{\alpha }_1}{{\mu }_1}+g+{\alpha }_2$

□

Biologically, the delay does not affect the system at tumor free equilibrium point because when tumor is free i.e. $T=0$ , the interaction between the immune cells and tumor cells disappear, so the system will undergo the same behavior whether there is a delay or there is not.

Theorem 4 System (1) is locally stable at $E_1$ for $\tau =0$ if and only if ${\mu }_2>{\varphi }_1$ and ${\mu }_3+{\alpha }_1>{\varphi }_2{\alpha }_1-{\alpha }_2$ , and undergoes bifurcation at $\tau ={\tau }^{\ast }$ where stability switches from stable to unstable.

Proof. Substituting $\left(N^{\ast },T^{\ast },M^{\ast }\right)=\left(\frac{1}{{\mu }_1}\left({\alpha }_1-{\varphi }_1{T}^{\ast }\right),T^{\ast },\frac{1}{{\gamma }_1{\gamma }_2{\mu }_1}\left({\gamma }_2{l}_2-k_1{T}^{\ast }\right)\right)$ , Equation (12) takes the form of

$J\left(E_1\right)=\left[\begin{array}{ccc}-{\alpha }_1+{\varphi }_1{T}^{\ast }& -{\varphi }_1{a}_1& 0\\ {\varphi }_2{T}^{\ast }& -a_2{\gamma }_1{\text{e}}^{-\lambda \tau }+{\varphi }_2{a}_1-2{\mu }_2{T}^{\ast }+g+{\alpha }_2& -{\gamma }_1{T}^{\ast }{\text{e}}^{-\lambda \tau }\\ 0& \frac{\rho a_2\omega }{{\left(\omega +T^{\ast }\right)}^2}-a_2{\gamma }_2{\text{e}}^{-\lambda \tau }& \frac{\rho T^{\ast }}{\omega +T^{\ast }}-{\gamma }_2{T}^{\ast }{\text{e}}^{-\lambda \tau }-\left({\mu }_3+g\right)\end{array}\right],$ (25)

where

$a_1=\frac{{\alpha }_1-{\varphi }_1{T}^{\ast }}{{\mu }_1}.$

$a_2=\frac{{\gamma }_2{l}_2-k_1{T}^{\ast }}{{\gamma }_1{\gamma }_2{\mu }_1}.$

$k_1={\gamma }_2\left({\mu }_1{\mu }_2+{\varphi }_1{\varphi }_2\right)\mathrm{.}$

$l_2={\mu }_1\left(g+{\alpha }_2\right)+{\alpha }_1{\varphi }_2$ .

Then, the simplified characteristic Equation of (25) becomes

${\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3+\left(B_1{\lambda }^2+B_2\lambda +B_3\right){\text{e}}^{-\lambda \tau }=0.$ (26)

where,

$\begin{array}{c}{A}_1=\frac{1}{\omega +T}(\left(2{\mu }_2-{\varphi }_1\right)T^2+\left(\left(2{\mu }_2-{\varphi }_1\right)\omega -{\varphi }_2{a}_1-{\alpha }_2-\rho +{\mu }_3+{\alpha }_1\right)T\\ -\omega \left({\varphi }_2{a}_1-{\alpha }_1+{\alpha }_2-{\mu }_3\right))\end{array}$

$\begin{array}{c}{A}_2=\frac{1}{\omega +T}(-2{\mu }_2{\varphi }_1{T}^3+(-2{\varphi }_1\omega {\mu }_2+\left(2{\mu }_2-{\varphi }_1\right){\mu }_3+2g{\mu }_2+2{\alpha }_1{\mu }_2+{\alpha }_2{\varphi }_1\\ +\left(2{\varphi }_2{a}_1+\rho \right){\varphi }_1-2\rho {\mu }_2)T^2+((\left(2{\mu }_2-{\varphi }_1\right){\mu }_3+2g{\mu }_2+2{\alpha }_1{\mu }_2\\ +{\varphi }_1\left(2{\varphi }_2{a}_1+{\alpha }_2\right))\omega +\left(-{\varphi }_2{a}_1-g+{\alpha }_1-{\alpha }_2\right){\mu }_3-g^2\\ +\left(-{\varphi }_2{a}_1+\rho -{\alpha }_2\right)g+\left(-{\varphi }_2{a}_1-{\alpha }_2\right){\alpha }_1+\rho \left({\varphi }_2{a}_1+{\alpha }_2\right))T\\ -\omega \left(\left({\varphi }_2{a}_1+g-{\alpha }_1+{\alpha }_2\right){\mu }_3+g^2+\left({\varphi }_2{a}_1+{\alpha }_2\right)g+{\alpha }_1\left({\varphi }_2{a}_1+{\alpha }_2\right)\right))\end{array}$

$\begin{array}{c}{A}_3=-\frac{1}{\omega +T}(2(T^2{\varphi }_1{\mu }_2+\left(-{\alpha }_1{\mu }_2-\frac{{\varphi }_1\left(2a_1{\varphi }_2+g+{\alpha }_2\right)}{2}\right)T\\ +\frac{{\alpha }_1\left(a_1{\varphi }_2+g+{\alpha }_2\right)}{2})\left(\left(g-\rho +{\mu }_3\right)T+\omega \left(g+{\mu }_3\right)\right))\end{array}$

$B_1=T^{\ast }{\gamma }_2+{\gamma }_1{a}_2$

$\begin{array}{c}{B}_2=\frac{1}{{\left(\omega +T\right)}^2}(2{\gamma }_2\left(-\frac{{\varphi }_1}{2}+{\mu }_2\right)T^4+(4{\gamma }_2\left(-\frac{{\varphi }_1}{2}+{\mu }_2\right)\omega +(-{\varphi }_2{a}_1-g+{\alpha }_1\\ \stackrel{}{}-{\alpha }_2){\gamma }_2-{\varphi }_1{\gamma }_1{a}_2)T^3+(2{\gamma }_2\left(-\frac{{\varphi }_1}{2}+{\mu }_2\right){\omega }^2+((-2{\varphi }_2{a}_1-2g+2{\alpha }_1\\ \stackrel{}{}-2{\alpha }_2){\gamma }_2-2{\varphi }_1{\gamma }_1{a}_2)\omega +{\gamma }_1{a}_2\left({\alpha }_1+g-\rho +{\mu }_3\right))T^2-\omega ((({\varphi }_2{a}_1+g-{\alpha }_1\\ \stackrel{}{}+{\alpha }_2){\gamma }_2+{\varphi }_1{\gamma }_1{a}_2)\omega -2{\gamma }_1{a}_2\left({\alpha }_1+g+{\mu }_3\right))T+{\gamma }_1{\omega }^2{a}_2\left({\alpha }_1+g+{\mu }_3\right))\end{array}$ (27)

$\begin{array}{c}{B}_3=\frac{1}{{\left(\omega +T\right)}^2}(-2T^5{\gamma }_2{\mu }_2{\varphi }_1+2\left(-2{\mu }_2{\varphi }_1\omega +{\alpha }_1{\mu }_2+{\varphi }_1\left(a_1{\varphi }_2+\frac{{\alpha }_2}{2}+\frac{g}{2}\right)\right){\gamma }_2{T}^4\\ +(-2{\varphi }_1{\gamma }_2{\omega }^2{\mu }_2+4\left({\alpha }_1{\mu }_2+{\varphi }_1\left(a_1{\varphi }_2+\frac{{\alpha }_2}{2}+\frac{g}{2}\right)\right){\gamma }_2\omega -{\gamma }_2\left(a_1{\varphi }_2+g+{\alpha }_2\right){\alpha }_1\\ \stackrel{}{}-{\varphi }_1{\gamma }_1{a}_2\left(g-\rho +{\mu }_3\right))T^3+(2\left({\alpha }_1{\mu }_2+{\varphi }_1\left(a_1{\varphi }_2+\frac{{\alpha }_2}{2}+\frac{g}{2}\right)\right){\gamma }_2{\omega }^2\\ \underset{}{}+\left(-2{\gamma }_2\left(a_1{\varphi }_2+g+{\alpha }_2\right){\alpha }_1-2{\varphi }_1{\gamma }_1{a}_2\left(g+{\mu }_3\right)\right)\omega +{\alpha }_1{\gamma }_1{a}_2\left(g-\rho +{\mu }_3\right))T^2\\ -\omega \left(\left({\gamma }_2\left(a_1{\varphi }_2+g+{\alpha }_2\right){\alpha }_1+{\varphi }_1{\gamma }_1{a}_2\left(g+{\mu }_3\right)\right)\omega -2{\alpha }_1{\gamma }_1{a}_2\left(g+{\mu }_3\right)\right)T\\ \underset{}{\overset{}{}}+{\alpha }_1{\gamma }_1{\omega }^2{a}_2\left(g+{\mu }_3\right))\end{array}$

When $\tau =0$ , the characteristic Equation (26) becomes

${\lambda }^3+\left(A_1+B_1\right){\lambda }^2+\left(A_2+B_2\right)\lambda +\left(A_3+B_3\right)=0.$ (28)

According to Routh Hurwitz criteria, sufficient and necessary conditions for system (1) to be asymptotically stable at $E_1$ are:

1) $A_1+B_1>0$ .

2) $A_3+B_3>0$ .

3) $\left(A_1+B_1\right)\left(A_2+B_2\right)-\left(A_3+B_3\right)>0$ .

Inspecting (27), $A_1+B_1>0$ if ${\mu }_2>{\varphi }_1$ and ${\mu }_3+{\alpha }_1>{\varphi }_2{\alpha }_1-{\alpha }_2$ , $A_3+B_3>0$ since $B_3>A_3$ and $\left(A_1+B_1\right)\left(A_2+B_2\right)-\left(A_3+B_3\right)>0$ since $\left(A_1+B_1\right)\left(A_2+B_2\right)$ is of higher polynomial. Hence, the conditions (1-3) are hold, then all eigenvalues in Equation (28) are negative which verifies the stability of the system for coexistence equilibrium.

When $\tau >0$ , here we analyze how the stability is affected by the time delay $\tau$ by considering $\tau$ as the bifurcation parameter. To study the delay induced instability, we assume a purely imaginary root of Equation (26) and we substitute $\lambda =i\sigma$ , ( $\sigma >0$ ) into Equation (26)

${\left(i\sigma \right)}^3+A_1{\left(i\sigma \right)}^2+A_2\left(i\sigma \right)+A_3+\left(B_1{\left(i\sigma \right)}^2+B_2\left(i\sigma \right)+B_3\right){\text{e}}^{-i\sigma \tau }=0$ (29)

Rewriting the exponential in terms of trigonometric functions we get

$-i{\sigma }^3-A_1{\sigma }^2+A_{2}i\sigma +A_3+\left(-B_1{\sigma }^2+B_{2}i\sigma +B_3\right)\left(\cos \left(\sigma \tau \right)-i\sin \left(\sigma \tau \right)\right)=0$ (30)

Separating real and imaginary parts of Equation (30) gives

$-A_1{\sigma }^2+A_3=\left(B_1{\sigma }^2-B_3\right)\cos \left(\sigma \tau \right)-B_2\sigma \sin \left(\sigma \tau \right)$ (31)

$-{\sigma }^3+A_2\sigma =\left(B_3-B_1{\sigma }^2\right)\sin \left(\sigma \tau \right)-B_2\sigma \cos \left(\sigma \tau \right)$ (32)

By squaring both sides of Equations (31) and (32) then adding these equations together, we get

${\sigma }^6+\left(A_1^2-2A_2-B_1^2\right){\sigma }^4+\left(A_2^2-2A_1{A}_3-B_2^2+2B_1{B}_3\right){\sigma }^2+\left(A_3^2-B_3^2\right)=0$ (33)

let $\nu ={\sigma }^2$ in Equation (33), then we have

${\nu }^3+\left(A_1^2-2A_2-B_1^2\right){\nu }^2+\left(A_2^2-2A_1{A}_3-B_2^2+2B_1{B}_3\right)\nu +\left(A_3^2-B_3^2\right)=0$ (34)

Since the leading coefficient of Equation (34) is positive, and inspecting the constant term $A_3^2-B_3^2$ from the definition in (27) suggests that $A_3^2-B_3^2<0$ . Hence Equation (34) has one positive root. Therefore, Hopf bifurcation occurs at $\tau ={\tau }^{*}$ near co-existence equilibrium.

For simplicity, we define

$p=A_1^2-2A_2-B_1^2$

$q=A_2^2-2A_1{A}_3-B_2^2+2B_1{B}_3$

$r=A_3^2-B_3^2$

then Equation (34) has the new form

${\nu }^3+p{\nu }^2+q\nu +r=0.$ (35)

The conditions for Equation (35) to have positive roots are stated in the following lemma [39] .

Lemma 1

1) If $r<0$ , then Equation (35) has at least one positive root.

2) If $r\ge 0$ , and $\Delta \le 0$ , then Equation (35) has no positive roots.

3) If $r\ge 0$ and $\Delta >0$ , then (35) has positive roots if and only if ${\nu }_1=\frac{1}{3}\left(-p+\sqrt{\Delta }\right)>0$ and $h\left({\nu }_1\right)\le 0$ .

Obviously, $r<0$ since $A_3^2-B_3^2<0$ , and there exist at least one positive root for Equation (35). Since condition (1) of Lemma 1 is satisfied, we conclude that Equation (35) has at least one positive root and there exists occurrence of Hopf bifurcation.

To determine the corresponding value of $\tau$ at which Hopf bifurcation occurs and the transversality condition therein. Doing this, we eliminate $\sin \left(\sigma \tau \right)$ from both Equations (31) and (32) as follows:

Multiply Equation (31) by $\left(B_3-B_1{\sigma }^2\right)$ and Equation (32) by $B_2\sigma$ then we add both equations as follows

$\left(B_3-B_1{\sigma }^2\right)\left(-A_1{\sigma }^2+A_3\right)=\left(B_3-B_1{\sigma }^2\right)\left[\left(B_1{\sigma }^2-B_3\right)\cos \left(\sigma \tau \right)-B_2\sigma \sin \left(\sigma \tau \right)\right]$ (36)

$\left(B_2\sigma \right)\left(-{\sigma }^3+A_2\sigma \right)=\left(B_2\sigma \right)\left[\left(B_3-B_1{\sigma }^2\right)\sin \left(\sigma \tau \right)-B_2\sigma \cos \left(\sigma \tau \right)\right]$ (37)

then, adding the equations in Equations (36) and (37), we obtain

$\begin{array}{l}\left(A_1{B}_1-B_2\right){\sigma }^4+\left(A_2{B}_2-A_1{B}_3-A_3{B}_1\right){\sigma }^2+A_3{B}_3\\ =\left[-{\left(B_1{\sigma }^2-B_3\right)}^2-B_2^2{\sigma }^2\right]\cos \left(\sigma \tau \right)\end{array}$ (38)

making $\cos \left(\sigma \tau \right)$ the subject in Equation (38), yields

$\cos \left(\sigma \tau \right)=\frac{\left(A_1{B}_1-B_2\right){\sigma }^4+\left(A_2{B}_2-A_1{B}_3-A_3{B}_1\right){\sigma }^2+A_3{B}_3}{-{\left(B_1{\sigma }^2-B_3\right)}^2-B_2^2{\sigma }^2}$ (39)

so the value of $\tau$ will be on the form

${\tau }^{\ast }=\frac{1}{\sigma }\arccos \left[\frac{\left(A_1{B}_1-B_2\right){\sigma }^4+\left(A_2{B}_2-A_1{B}_3-A_3{B}_1\right){\sigma }^2+A_3{B}_3}{-{\left(B_1{\sigma }^2-B_3\right)}^2-B_2^2{\sigma }^2}\right]$ (40)

or

$\begin{array}{l}{\tau }_n^{\ast }=-\frac{1}{{\sigma }_0}\arccos \left[\frac{\left(A_1{B}_1-B_2\right){\sigma }_0^4+\left(A_2{B}_2-A_1{B}_3-A_3{B}_1\right){\sigma }_0^2+A_3{B}_3}{-{\left(B_1{\sigma }_0^2-B_3\right)}^2-B_2^2{\sigma }_0^2}\right]\\ +\frac{2n\pi }{{\sigma }_0},n=0,1,2,\cdots \end{array}$ (41)

□

Occurrence of Hopf Bifurcation

To verify the direction of the switch when the Hopf bifurcation occur i.e. when the delay term $\tau$ passes through the imaginary axis at critical value $\tau ={\tau }^{\ast }$ , we desire to show that

$\frac{\text{d}\left(Re\lambda \right)}{\text{d}\tau }\ne 0$

which is called the transversality condition for the Hopf bifurcation at ${\tau }_n$ , so to prove the transversality condition, Equation (26) is differentiated with respect to $\tau$ to give

$\begin{array}{l}\frac{\text{d}\lambda }{\text{d}\tau }\left[3{\lambda }^2+2A_1\lambda +A_2+2B_1\lambda {\text{e}}^{-\lambda \tau }-B_1{\lambda }^2\tau {\text{e}}^{-\lambda \tau }-B_2\lambda \tau {\text{e}}^{-\lambda \tau }-B_3\tau {\text{e}}^{-\lambda \tau }\right]\\ =B_1{\lambda }^3{\text{e}}^{-\lambda \tau }+B_2{\lambda }^2{\text{e}}^{-\lambda \tau }+B_3\lambda {\text{e}}^{-\lambda \tau }\end{array}$ (42)

$\begin{array}{l}\frac{\text{d}\lambda }{\text{d}\tau }\left[3{\lambda }^2+2A_1\lambda +A_2+2B_1\lambda {\text{e}}^{-\lambda \tau }-B_1{\lambda }^2\tau {\text{e}}^{-\lambda \tau }-B_2\lambda \tau {\text{e}}^{-\lambda \tau }-B_3\tau {\text{e}}^{-\lambda \tau }\right]\\ =\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda {\text{e}}^{-\lambda \tau }\end{array}$

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{3{\lambda }^2+2A_1\lambda +A_2+2B_1\lambda {\text{e}}^{-\lambda \tau }-B_1{\lambda }^2\tau {\text{e}}^{-\lambda \tau }-B_2\lambda \tau {\text{e}}^{-\lambda \tau }-B_3\tau {\text{e}}^{-\lambda \tau }}{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda {\text{e}}^{-\lambda \tau }}$

$\begin{array}{c}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{3{\lambda }^2+2A_1\lambda +A_2}{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda e^{-\lambda \tau }}+\frac{2B_1\lambda {\text{e}}^{-\lambda \tau }}{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda {\text{e}}^{-\lambda \tau }}\\ -\frac{B_1{\lambda }^2\tau {\text{e}}^{-\lambda \tau }+B_2\lambda \tau {\text{e}}^{-\lambda \tau }+B_3\tau {\text{e}}^{-\lambda \tau }}{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda {\text{e}}^{-\lambda \tau }}\end{array}$

$\begin{array}{c}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{3{\lambda }^2+2A_1\lambda +A_2}{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda {\text{e}}^{-\lambda \tau }}+\frac{2B_1\lambda {\text{e}}^{-\lambda \tau }}{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda {\text{e}}^{-\lambda \tau }}\\ -\frac{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\tau {\text{e}}^{-\lambda \tau }}{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda {\text{e}}^{-\lambda \tau }}\end{array}$

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{3{\lambda }^2+2A_1\lambda +A_2}{\left(B_1{\lambda }^2+B_2\lambda +B_3\right)\lambda {\text{e}}^{-\lambda \tau }}+\frac{2B_1}{B_1{\lambda }^2+B_2\lambda +B_3}-\frac{\tau }{\lambda }$ (43)

rewriting the characteristic Equation (26) to have ${\text{e}}^{-\lambda \tau }$ as a subject, we have

${\text{e}}^{-\lambda \tau }=-\frac{{\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3}{B_1{\lambda }^2+B_2\lambda +B_3}$ (44)

Substituting (44) into (43), we obtain

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{3{\lambda }^2+2A_1\lambda +A_2}{-\lambda \left({\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3\right)}+\frac{2B_1}{B_1{\lambda }^2+B_2\lambda +B_3}-\frac{\tau }{\lambda }$ (45)

Therefore

$\begin{array}{l}sign{\left\{\frac{\text{d}\left(Re\lambda \right)}{\text{d}\tau }\right\}}_{\tau ={\tau }_n}=sign{\left\{Re{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}}_{\lambda =i\sigma }\\ =sign\left\{Re{\left[\frac{3{\lambda }^2+2A_1\lambda +A_2}{-\lambda \left({\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3\right)}+\frac{2B_1}{B_1{\lambda }^2+B_2\lambda +B_3}-\frac{\tau }{\lambda }\right]}_{\lambda =i\sigma }\right\}\\ =sign\{Re{\left[\frac{3{\lambda }^2+2A_1\lambda +A_2}{-\lambda \left({\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3\right)}\right]}_{\lambda =i\sigma }+Re{\left[\frac{2B_1}{B_1{\lambda }^2+B_2\lambda +B_3}\right]}_{\lambda =i\sigma }\\ -Re{\left[\frac{\tau }{\lambda }\right]}_{\lambda =i\sigma }\}\end{array}$ (46)

Now, substituting by $\lambda =i\sigma$ into (46), we obtain

$\begin{array}{l}sign{\left\{Re{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}}_{\lambda =i\sigma }\\ =sign\{Re{\left[\frac{3{\left(i\sigma \right)}^2+2A_1\left(i\sigma \right)+A_2}{-\left(i\sigma \right)\left({\left(i\sigma \right)}^3+A_1{\left(i\sigma \right)}^2+A_2\left(i\sigma \right)+A_3\right)}\right]}_{\lambda =i\sigma }\\ +Re{\left[\frac{2B_1}{\mathrm{<mo\; stretchy="false">\; (\; </mo>}{B}_1{\left(i\sigma \right)}^2+B_2\left(i\sigma \right)+B_3\mathrm{<mo\; stretchy="false">\; )\; </mo>}}\right]}_{\lambda \mathrm{<mo>\; =\; </mo>}i\sigma }-Re{\left[\frac{\tau }{i\sigma }\right]}_{\lambda =i\sigma }\}\end{array}$ (47)

$=sign\left\{\frac{-3{\sigma }^6+\left(2A_2+2A_1^2\right){\sigma }^4+\left(2A_1{A}_3+A_2^2\right){\sigma }^2}{{\left(A_2{\sigma }^2-{\sigma }^4\right)}^2+{\left(A_3\sigma -A_1{\sigma }^3\right)}^2}-\frac{2B_1\left(B_3-B_1{\sigma }^2\right)}{{\left(B_3-B_1{\sigma }^2\right)}^2+{\left(B_2\sigma \right)}^2}\right\}$ (48)

$=sign\left\{\frac{-3{\sigma }^4+\left(2A_2+2A_1^2\right){\sigma }^2+\left(2A_1{A}_3+A_2^2\right)}{{\sigma }^6+\left(A_1^2-2A_2\right){\sigma }^4+\left(A_2^2-2A_1{A}_3\right){\sigma }^2+A_3^2}-\frac{2B_1\left(B_3-B_1{\sigma }^2\right)}{{\left(B_3-B_1{\sigma }^2\right)}^2+{\left(B_2\sigma \right)}^2}\right\}$ (49)

from Equation (33), we know that

$\begin{array}{l}{\sigma }^6+\left(A_1^2-2A_2\right){\sigma }^4+\left(A_2^2-2A_1{A}_3\right){\sigma }^2+A_3^2\\ =B_1^2{\sigma }^4+\left(B_2^2+2B_1{B}_3\right){\sigma }^2+B_3^2={\left(B_3-B_1{\sigma }^2\right)}^2+{\left(B_2\sigma \right)}^2\end{array}$

substituting in Equation (49)

$\begin{array}{l}sign{\left\{Re{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}}_{\lambda =i\sigma }\\ =sign\{\frac{-3{\sigma }^4+\left(2A_2+2A_1^2\right){\sigma }^2+\left(2A_1{A}_3+A_2^2\right)}{{\left(B_3-B_1{\sigma }^2\right)}^2+{\left(B_2\sigma \right)}^2}-\frac{2B_1\mathrm{<mo\; stretchy="false">\; (\; </mo>}{B}_3-B_1{\sigma }^2\mathrm{<mo\; stretchy="false">\; )\; </mo>}}{{\left(B_3-B_1{\sigma }^2\right)}^2+{\left(B_2\sigma \right)}^2}\}\\ =sign\left\{\frac{-3{\sigma }^4+\left(2A_2+2A_1^2\right){\sigma }^2+\left(2A_1{A}_3+A_2^2\right)-2B_1\left(B_3-B_1{\sigma }^2\right)}{{\left(B_3-B_1{\sigma }^2\right)}^2+{\left(B_2\sigma \right)}^2}\right\}\end{array}$

therefore

$\begin{array}{l}sign{\left\{Re{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}}_{\lambda =i\sigma }\\ =sign\left\{-\left[\frac{3{\sigma }^4-\left(2A_2+2A_1^2+2B_1^2\right){\sigma }^2-\left(2A_1{A}_3+A_2^2-2B_1{B}_3\right)}{{\left(B_3-B_1{\sigma }^2\right)}^2+{\left(B_2\sigma \right)}^2}\right]\right\}\end{array}$ (50)

Since $Sign\left\{Re{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}\right\}\ne 0$ then the transversality condition holds. Therefore, Hopf-bifurcation occurs at ${\tau }_n^{\ast }$ where the stability switches from stable to unstable.

4. Numerical Simulation

The numerical simulation is applying by using a Matlab version R2021b with DDE23 solver. DDE-Biftool, which is a Matlab package is used for the bifurcation analysis of this work. To calculate the eigenvalues numerically with the parameters values in Table 1, we use Maple software. Referring to [36] , we have two cases: tumour free equlibrium point and coexisting equilibrium point. Each case has its own parameters values.

4.1. Case 1: Tumor-Free Equilibrium Point E₀

For tumor-free equilibrium point, we set the following parameters values as follows: ${\alpha }_1=0.7$ , ${\mu }_1=0.1$ , ${\varphi }_1=0.1$ , ${\alpha }_2=0.98$ , ${\mu }_2=0.4$ , ${\gamma }_1=0.8$ , $\rho =0.3$ , $\omega =0.3$ , ${\gamma }_2=0.29$ , ${\mu }_3=0.15$ , ${\varphi }_2=0.5$ , and $s=0.6$ , $g=0.15$ . Substituting by these values we obtain the tumor free equilibrium point $E_0=\left(0.875,0,1.966\right)$ . If $\tau =0$ , system (1) is unstable as shown in Figure 2 since the condition (15) is not satisfied and the eigenvalues which obtained are ${\lambda }_1=-0.7$ , ${\lambda }_2=-0.3$ and ${\lambda }_3=0.870$ are obtained.

Figure 2. Cells population when $\tau =0$ at the tumour-free-equilibrium-point $E_0=\left(0.875,0,1.966\right)$ with the initial conditions $\left(\mathrm{0.01,0,0.01}\right)$ when (a) $g=0.01$ , (b) $g=0.15$ , (c) $g=0.5$ .

Table 1. Model parameters values.

When time delay is introducing to the system, we observed that there is no change in cells behavior which means that the delay will not make any change at the tumor free equilibrium point.

Also, the numerical substitution verified that there is no positive roots of $\theta$ since Equation (24) gives $E^2-D^2=-18.8769<0$ , then $\tau$ will not switch its stability, and as result Hopf bifurcation does not occur.

It is obvious that from Figure 2 and Figure 3, the delay did not impact the behavior of the cells in the system. Delay and non-delay have exactly the same graphs. Figure 2 and Figure 3 depicts the interacting cells dynamics for tumor-free equilibrium with variation of glucose concentration g at $\tau \ge 0$ . (a) is when $g=0.01$ , the normal cells grow until it reaches carrying capacity, and the immune cells only seen progressing initially before attaining a steady-state. (b) is when $g=0.15$ , the scenario is relatively as in case (a) except lower growth of immune cells. (c) is when $g=0.5$ , the normal cell exhibits the same scenario as in (a & b), but immune cells is relatively dormant comparing to (a & b). The growth of one of the interacting cells (normal cells) to their carrying capacity suggest an unstable scenario and this is tally with the analytical results. Also, the immune cells growth is affected by glucose rate change. This is because the immune efficiency are highly affected by glucose rate, the higher level of glucose in the blood, the lower efficiency of the immune system, as stated in the biological researches [17] [18] .

4.2. Case 2: Coexisting Equilibrium Point E₁

For coexisting equilibrium point, we use the following parameters values in the numerical simulation: ${\alpha }_1=0.7$ , ${\mu }_1=0.8$ , ${\varphi }_1=0.1$ , ${\alpha }_2=0.98$ , ${\mu }_2=0.4$ , ${\gamma }_1=0.8$ , $\rho =0.8$ , $\omega =0.3$ , ${\gamma }_2=0.29$ , ${\mu }_3=0.15$ , ${\varphi }_2=0.1$ , and $s=0.4$ , $g=0.3992$ . Substituting by these values we obtain the coexisting equilibrium point $E_1=\left(0.8377,0.2983,1.6889\right)$ .

If $\tau =0$ , system (1) is stable since the conditions in Theorem 4 are satisfied, as shown in Figure 4(b). When $\tau$ increases and becomes greater than zero, the system switches its case from stability to unstability.

Figure 3. Cells population when $\tau >0$ at the tumour-free-equilibrium-point $E_0=\left(0.875,0,1.966\right)$ with the initial conditions $\left(\mathrm{0.01,0,0.01}\right)$ when (a) $g=0.01$ , (b) $g=0.15$ , (c) $g=0.5$ .

Figure 4. Cells population when $\tau =0$ at the coexisting equilibrium-point $E_1=\left(\mathrm{0.8377,0.2983,1.6889}\right)$ with the initial conditions $\left(\mathrm{0.83,0.29,1.68}\right)$ when (a) $g=0.01$ , (b) $g=0.3922$ , (c) $g=0.5$ .

The values of $\tau$ are obtained by solving Equation (34) as follows:

${\nu }^3-0.1010209194{\nu }^2-0.3495525239\nu -0.04561671822=0$ (51)

The solutions of this equation are: ${\nu }_1=0.6967132133$ , ${\nu }_2=-0.1454054735$ , ${\nu }_3=-0.45028682046$ , so there exists one positive root. Since $\sigma =\pm \sqrt{\nu }$ , we will get the following values of $\sigma$ :

$\sigma =\pm \sqrt{0.6967132133}=\pm 0.8346934846$ .

Then substituting by the values of $\sigma$ in Equation (41) using the parameters values generates the following $\tau$ values: ${\tau }_0=0.574$ , ${\tau }_1=8.101$ , ${\tau }_2=15.629$ , ${\tau }_3=23.156$ , ${\tau }_4=30.684$ , ${\tau }_5=38.211$ , and ${\tau }_6=45.739$ .

Figure 4 depicts the interacting cells dynamics for coexistence equilibrium for $\tau =0$ with variation of glucose concentrations g. (a) is when $g=0.01$ , tumors cells regress as time continuous to zero, the immune cells show progressing and undulate to attain stable state. Normal cells demonstrate stable and steady state throughout. (b) is when $g=0.3922$ , the interacting cells are seen stable all through. (c) is when $g=0.5$ , tumor cells is seen progressing towards its carrying capacity, while immune cells regress to overlapped the regressing normal cells as a result of the glucose excess which is verified by the biological results [6] - [15] .

Figure 5 depicts the dynamics of the interacting cells for coexistence equilibrium when $\tau =0.574$ with variation of glucose concentration g. We observe from the figure the stability of the system starts to change and the cells behavior begin to be unstable. (a) is when $g=0.01$ , the stability switches from stable to unstable. Both the tumor and immune curves slide to zero but normal cell curve keep progressing. (b) is when $g=0.3922$ , the lines of stability began to ripple, which means that there is a change in the behavior of cells from stability to instability. (c) is when $g=0.5$ , the interacting cells dynamics show that the tumour grow up while immune cells and normal began to decline as a result of the glucose excess which is consistent with the biological results [6] - [15] .

Figure 6 depicts the dynamics of the interacting cells for coexistence equilibrium when $\tau =8.101$ with variation of glucose concentration g. (a) is when $g=0.01$ , stability of the interacting cells characterized with a fall of immune cells’ curve to zero from the steadiness obtained when $\tau =0$ (b) is when $g=0.3922$ , delay induces unregulated in immune cells resurgence, the normal cells curve is apparently flatten, and tumor cell curve exhibits progression to a dormant state before descending to zero, suggesting a delay-induced change of stability (c) is when $g=0.5$ , the interacting cells dynamics is relatively the same with (b) with exception of higher concentration of the interacting cells, that is a delay-induced change of stability also characterized this case as in (b).

Figures 7-11 have the same behavior when $\tau =8.101$ because the stability of the system will continue to be unstable and will not switch again to be stable for all delay values.

Figure 12 describes numerically, the bifurcation analysis of the coexistence equilibrium point for system (1). When $\tau ={\tau }_0$ , the stability switching from stable to unstable where Hopf bifurcation occurs as shown in Figure 12(a) and condition (50) is hold. However, when $\tau$ increases, the system continues to be unstable as shown in Figure 12(b). Figure 12(a) provides stability details of the coexistence equilibrium with the green and the red lines in representing the stable and unstable parts of the branches respectively. Figure 12(b) indicates all the critical values of $\tau$ . Figure 12(c) depicts the stable branch emanating from Hopf bifurcation (in green) for co-existence equilibrium.

Figure 5. Cells population when $\tau =0.574$ the coexisting equilibrium-point $E_1=\left(\mathrm{0.8377,0.2983,1.6889}\right)$ with the initial conditions $\left(0.83,0.29,1.68\right)$ when (a) $g=0.01$ , (b) $g=0.3922$ , (c) $g=0.5$ .

Figure 6. Cells population when $\tau =8.101$ at the coexisting equilibrium-point $E_1=\left(\mathrm{0.8377,0.2983,1.6889}\right)$ with the initial conditions $\left(0.83,0.29,1.68\right)$ when (a) $g=0.01$ , (b) $g=0.3922$ , (c) $g=0.5$ .

Figure 7. Cells population when $\tau =15.629$ at the coexisting equilibrium-point $E_1=\left(\mathrm{0.8377,0.2983,1.6889}\right)$ with the initial conditions $\left(0.83,0.29,1.68\right)$ when (a) $g=0.01$ , (b) $g=0.3922$ , (c) $g=0.5$ .

Figure 8. Cells population when $\tau =23.156$ at the coexisting equilibrium-point $E_1=\left(\mathrm{0.8377,0.2983,1.6889}\right)$ with the initial conditions $\left(0.83,0.29,1.68\right)$ when (a) $g=0.01$ , (b) $g=0.3922$ , (c) $g=0.5$ .

Figure 9. Cells population when $\tau =30.684$ at the coexisting equilibrium-point $E_1=\left(\mathrm{0.8377,0.2983,1.6889}\right)$ with the initial conditions $\left(0.83,0.29,1.68\right)$ when (a) $g=0.01$ , (b) $g=0.3922$ , (c) $g=0.5$ .

Figure 10. Cells population when $\tau =38.211$ at the coexisting equilibrium-point $E_1=\left(\mathrm{0.8377,0.2983,1.6889}\right)$ with the initial conditions $\left(0.83,0.29,1.68\right)$ when (a) $g=0.01$ , (b) $g=0.3922$ , (c) $g=0.5$ .

Figure 11. Cells population when $\tau =45.739$ at the coexisting equilibrium-point $E_1=\left(\mathrm{0.8377,0.2983,1.6889}\right)$ with the initial conditions $\left(0.83,0.29,1.68\right)$ when (a) $g=0.01$ , (b) $g=0.3922$ , (c) $g=0.5$ .

Figure 12. (a) Stability switching at $\tau =0.57$ such that green line represents stable and red line represents unstable, respectively, at coexisting equilibrium-point $E_1=\left(0.8377,0.2983,1.6889\right)$ when $g=0.3922$ . (b) All critical values of $\tau$ are shown in this figure. (c) Hopf bifurcation occurred at $\tau =0.57$ and a doubling bifurcation occurred at $\tau =0.605$ .

5. Discussion and Conclusions

The main objective of this work is to examine the effect of time delay occasioned by biological process in the derivation of cytokines for immune-suppression by the tumor, and the secretion of cytokines by the immune in the dynamics of interaction among tumor, immune and normal cells in the breast cancer model. To form basis for biological meaningfulness of System (1) outcomes, no-negativity and boundlessness of the solution are obtained, and so also, biological equilibrium points, namely; tumor-free and coexisting equilibrium point were found. The local stability and Hopf bifurcation analysis for these points were studied to obtain both stability conditions and possible occurrence of stability switch. The numerical simulation validated our theoretical results.

The analytical results suggest that the tumor-free equilibrium point is always unstable as normal cells are bound to grow and reach their carrying capacity regardless of delays. Also, it was observed that there is no change in cells behavior when the delay was introduced to the system and when it was not. Biologically, since the tumor equals to zero which means no existence of disease, the interaction between immune cells and tumor cells disappears, so the delay will not affect the system.

For coexistence equilibrium point, the result suggests as follows:

1) For $\tau =0$ , coexistence equilibrium point is stable scenario if the tumor death rate is greater than tumor impairment on normal cells growth.

2) For $\tau >0$ , coexistence equilibrium point, there exists occurrence of Hopf bifurcation with switch of stability from stable to unstable at a critical value of $\tau ={\tau }^{*}$ .

Numerical results validate all the analytical results and further provided graphical exhibition of the interacting cells under variation of glucose concentrations. For tumor-free equilibrium point, the normal cell is seen growing to reach their carrying capacity, regardless of increase in glucose concentrations (see Figure 2). The immune cells apparently less active and the concentration reduce as glucose increases (see Figure 2). For coexistence equilibrium point with $\tau =0$ , the interacting cells are relatively stable when the glucose concentration is 0.01 and 0.3922, but unstable when the glucose concentration surges to 0.5 as in Figure 4. When $\tau =0.574$ for the coexistence equilibrium point, the delay induces change stability from stable to unstable when the glucose concentration is 0.01 and 0.3922 (see Figure 5(a), Figure 5(b)), but the delay has no effect when the glucose concentration is 0.5 as shown in Figure 5(c). When $\tau =8.101$ and other values of $\tau$ , the cells behavior has changed and become completely different from when $\tau =0.574$ and it is obvious from Figures 6-10 the system is unstable at some g (glucose) values. Also, the variation of $\tau$ values does not affect the stability when the glucose concentration is 0.01, but the delay induces change stability from stable to unstable when the glucose concentration is 0.3922 and 0.5. Biologically, the delay in secretion of cytokines by immune cells and derivation of cytokines by the tumors cells depicts by $\tau =0.574$ can revive the normal cells while both the tumor and immune cells remain stable (see Figure 5(a)). However, the continuous in delay might affect no change when the glucose level is small (see Figures 6(a)-11(a)), but it might brew a chaotic unregulated growth of immune cells (see Figures 6(b)-11(b) & Figures 6(c)-11(c)).

In conclusion, the delay in the secretion of cytokines by immune cells and derivation cytokines by the tumors helps to identify the possible chaotic situation under different glucose concentration and the extent to which such delay can have on restoration of the normal cells when glucose concentration is low. The numerical results also hint at possible delay-induced chaotic situation when glucose concentration increases on the immune cells, that is the unregulated growth of immune cells. In our next work, we will look at incorporating an immunotherapy into our proposed model to avert such un-regulatory growth of the immune cells when there is increase in both delay and glucose concentration.

Acknowledgements

The authors are thankful to Universiti Teknologi Malaysia for providing the facilities in this research.

Funding

This research received the Research Management Center (UTM) for financial support through research grants of vote Q.J130000.2554.21H19.

Availability of Data and Materials

The data used to support the findings of this study are included in the article.

Authors’ Contributions

Conceptualization, A.A.; methodology A.A.; software, A.A., N.A; validation, A.A., N.M.; writing-original draft preparation, A.A.; writing-review and editing, N.M.; supervision, N.M.; project administration, N.M.; funding acquisition, A.A. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest

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

References

[1]	Monticciolo, D.L. (2022) Invited Commentary: The Challenges of Early-Onset Breast Cancer. RadioGraphics, 42, E16-E17. https://doi.org/10.1148/rg.210191
[2]	Simms, K.T. (2012) Mathematical Models of Cell Cycle Progression: Applications to Breast Cancer Cell Lines. The University of Adelaide, Adelaide.
[3]	Yousef, A., Bozkurt, F. and Abdeljawad, T. (2020) Mathematical Modeling of the Immune-Chemotherapeutic Treatment of Breast Cancer under some Control Parameters. Advances in Difference Equations, 2020, Article No. 696. https://doi.org/10.1186/s13662-020-03151-5
[4]	Oke, S.I., Matadi, M.B. and Xulu, S S. (2018) Optimal Control Analysis of a Mathematical Model for Breast Cancer. Mathematical and Computational Applications, 23, Article No. 21. https://doi.org/10.3390/mca23020021
[5]	Sun, Y.-S., Zhao, Z., Yang, Z.-N., Xu, F., Lu, H.-J., Zhu, Z.-Y., et al. (2017) Risk Factors and Preventions of Breast Cancer. International Journal of Biological Sciences, 13, 1387-1397. https://doi.org/10.7150/ijbs.21635
[6]	Santos, J.M. and Hussain, F. (2020) Higher Glucose Enhances Breast Cancer Cell Aggressiveness. Nutrition and Cancer, 72, 734-746. https://doi.org/10.1080/01635581.2019.1654527
[7]	Fadaka, A., et al. (2017) Biology of Glucose Metabolization in Cancer Cells. Journal of Oncological Sciences, 3, 45-51. https://doi.org/10.1016/j.jons.2017.06.002
[8]	Vander Heiden, M.G., Cantley, L.C. and Thompson, C.B. (2009) Understanding the Warburg Effect: The Metabolic Requirements of Cell Proliferation. Science, 324, 1029-1033. https://doi.org/10.1126/science.1160809
[9]	Chen, M.-C., et al. (2020) ROS Mediate xCT-Dependent Cell Death in Human Breast Cancer Cells under Glucose Deprivation. Cells, 9, Article No. 1598. https://doi.org/10.3390/cells9071598
[10]	Wardi, L., et al. (2014) Glucose Restriction Decreases Telomerase Activity and Enhances Its Inhibitor Response on Breast Cancer Cells: Possible Extra-Telomerase Role of BIBR 1532. Cancer Cell International, 14, Article No. 60. https://doi.org/10.1186/1475-2867-14-60
[11]	Barbosa, A.M. and Martel, F. (2020) Targeting Glucose Transporters for Breast Cancer Therapy: The Effect of Natural and Synthetic Compounds. Cancers, 12, Article No. 154. https://doi.org/10.3390/cancers12010154
[12]	Sun, S., Sun, Y., Rong, X. and Bai, L. (2019) High Glucose Promotes Breast Cancer Proliferation and Metastasis by Impairing Angiotensinogen Expression. Bioscience Reports, 39, Article ID: BSR20190436. https://doi.org/10.1042/BSR20190436
[13]	Krętowski, R., et al. (2016) Low Glucose Dependent Decrease of Apoptosis and Induction of Autophagy in Breast Cancer MCF-7 Cells. Molecular and Cellular Biochemistry, 417, 35-47. https://doi.org/10.1007/s11010-016-2711-4
[14]	O’Mahony, F., et al. (2012) Estrogen Modulates Metabolic Pathway Adaptation to Available Glucose in Breast Cancer Cells. Molecular Endocrinology, 26, 2058-2070. https://doi.org/10.1210/me.2012-1191
[15]	Qiu, J., Zheng, Q. and Meng, X. (2021) Hyperglycemia and Chemoresistance in Breast Cancer: From Cellular Mechanisms to Treatment Response. Frontiers in Oncology, 11, Article 628359. https://doi.org/10.3389/fonc.2021.628359
[16]	Duan, W., Shen, X., Lei, J., Xu, Q., Yu, Y., Li, R., Wu, E. and Ma, Q. (2014) Hyperglycemia, a Neglected Factor during Cancer Progression. BioMed Research International, 2014, Article ID: 461917. https://doi.org/10.1155/2014/461917
[17]	Shomali, N., et al. (2020) Harmful Effects of High Amounts of Glucose on the Immune System: An Updated Review. Biotechnology and Applied Biochemistry, 68, 404-410. https://doi.org/10.1002/bab.1938
[18]	Von Ah Morano, A.E., Dorneles, G.P., Peres, A. and Lira, F.S. (2020) The Role of Glucose Homeostasis on Immune Function in Response to Exercise: The Impact of Low or Higher Energetic Conditions. Journal of Cellular Physiology, 235, 3169-3188. https://doi.org/10.1002/jcp.29228
[19]	Jeandidier, N. and Boullu-Sanchis, S. (2006) Hyperglycemia and Acute Illness. Annales d’Endocrinologie, 67, 224-232. https://doi.org/10.1016/S0003-4266(06)72590-2
[20]	Berbudi, A., Rahmadika, N., Tjahjadi, A.I. and Ruslami, R. (2020) Type 2 Diabetes and Its Impact on the Immune System. Current Diabetes Reviews, 16, 442-449. https://doi.org/10.2174/1573399815666191024085838
[21]	Rihan, F.A., Abdelrahman, D.H., Al-Maskari, F., Ibrahim, F. and Abdeen, M.A. (2014) Delay Differential Model for Tumour-Immune Response with Chemoimmunotherapy and Optimal Control. Computational and Mathematical Methods in Medicine, 2014, Article ID: 982978. https://doi.org/10.1155/2014/982978
[22]	Chapman, M.P. and Tomlin, C.J. (2016) Ordinary Differential Equations in Cancer Biology. BioRxiv: 071134. https://doi.org/10.1101/071134
[23]	Newbury, G. (2007) A Numerical Study of a Delay Differential Equation Model for Breast Cancer. Virginia Polytechnic Institute and State University, Virginia.
[24]	Shampine, L.F., Thompson, S. and Kierzenka, J. (2000) Solving Delay Differential Equations with dde23. URL. https://scholar.google.com/scholar?hl=en&as_sdt=0%2C5&q=Solving+delay+differential+equations+with+dde23&btnG=
[25]	Khajanchi, S. (2020) Chaotic Dynamics of a Delayed Tumor—Immune Interaction Model. International Journal of Biomathematics, 13, Article ID: 2050009. https://doi.org/10.1142/S1793524520500096
[26]	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
[27]	Topalian, S.L., Hodi, F.S., Brahmer, J.R., Gettinger, S.N., Smith, D.C., McDermott, D.F., et al. (2012) Safety, Activity, and Immune Correlates of Anti-PD-1 Antibody in Cancer. New England Journal of Medicine, 366, 2443-2454. https://doi.org/10.1056/NEJMoa1200690
[28]	Udomchalermpat, S., Koonprasert, S. and Kunnawuttipreechachan, E. (2020) Dynamics of the Tumor-Obesity with Time Delay Effect. International Journal of Engineering Research and Technology, 13, 1854-1865. https://doi.org/10.37624/IJERT/13.8.2020.1854-1865
[29]	Ku-Carrillo, R.A., Delgadillo, S.E. and Chen-Charpentier, B.M. (2016) A Mathematical Model for the Effect of Obesity on Cancer Growth and on the Immune System Response. Applied Mathematical Modelling, 40, 4908-4920. https://doi.org/10.1016/j.apm.2015.12.018
[30]	Kuznetsov, V.A., Makalkin, I.A., Taylor, M.A. and Perelson, A.S. (1994) Nonlinear Dynamics of Immunogenic Tumors: Parameter Estimation and Global Bifurcation Analysis. Bulletin of Mathematical Biology, 56, 295-321.
[31]	Awang, N.A., Maan, N. and Sulain, M.D. (2022) Tumour-Natural Killer and CD8+ T Cells Interaction Model with Delay. Mathematics, 10, Article No. 2193. https://doi.org/10.3390/math10132193
[32]	Ibrahim, A.A., Maan, N., Jemon, K. and Abidemi, A. (2022) Global Stability and Thermal Optimal Control Strategies for Hyperthermia Treatment of Malignant Tumors. Mathematics, 10, Article No. 2188. https://doi.org/10.3390/math10132188
[33]	Khajanchi, S., Perc, M. and Ghosh, D. (2018) The Influence of Time Delay in a Chaotic Cancer Model. Chaos: An Interdisciplinary Journal of Nonlinear Science, 28, Article ID: 103101. https://doi.org/10.1063/1.5052496
[34]	Rihan, F.A., Tunc, C., Saker, S.H., Lakshmanan, S. and Rakkiyappan, R. (2018) Applications of Delay Differential Equations in Biological Systems. Complexity, 2018, Article ID: 4584389. https://doi.org/10.1155/2018/4584389
[35]	Khajanchi, S. (2015) Bifurcation Analysis of a Delayed Mathematical Model for Tumor Growth. Chaos, Solitons & Fractals, 77, 264-276. https://doi.org/10.1016/j.chaos.2015.06.001
[36]	Alblowy, A.H., Maan, N. and Alharbi, S.A. (2022) Role of Glucose Risk Factors on Human Breast Cancer: A Nonlinear Dynamical Model Evaluation. Mathematics, 10, Article No. 3640. https://doi.org/10.3390/math10193640
[37]	Forde, J. and Nelson, P. (2004) Applications of Sturm Sequences to Bifurcation Analysis of Delay Differential Equation Models. Journal of Mathematical Analysis and Applications, 300, 273-284. https://doi.org/10.1016/j.jmaa.2004.02.063
[38]	Banerjee, S. and Sarkar, R.R. (2008) Delay-Induced Model for Tumor-Immune Interaction and Control of Malignant Tumor Growth. Biosystems, 91, 268-288. https://doi.org/10.1016/j.biosystems.2007.10.002
[39]	Song, Y. and Yuan, S. (2006) Bifurcation Analysis in a Predator-Prey System with Time Delay. Nonlinear Analysis: Real World Applications, 7, 265-284. https://doi.org/10.1016/j.nonrwa.2005.03.002
[40]	Mufudza, C., Sorofa, W. and Chiyaka, E.T. (2012) Assessing the Effects of Estrogen on the Dynamics of Breast Cancer. Computational and Mathematical Methods in Medicine, 2012, Article ID: 473572. https://doi.org/10.1155/2012/473572
[41]	Das, A., Sarmah, H.K., Bhattacharya, D., Dehingia, K. and Hosseini, K. (2022) Combination of Virotherapy and Chemotherapy with Optimal Control for Combating Cancer. Mathematics and Computers in Simulation, 194, 460-488. https://doi.org/10.1016/j.matcom.2021.12.004