Analysis of a delay-induced mathematical model of cancer

In this paper, 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. Though the model was originally proposed by De Pillis et al. (Math. Comput. Model. 37:1221–1244, 2003), to make the model more realistic, we have added a delay term into the model, and it has incorporated novelty in our present work. The stability of existing equilibrium points in the delay-induced system is studied in detail. Global stability conditions of the tumor-free equilibrium point have been found. It is shown that due to this delay effect, the coexisting equilibrium point may lose its stability through a Hopf bifurcation. The implicit function theorem is applied to characterize a complex function in a neighborhood of delay terms. Additionally, the presence of Hopf bifurcation is demonstrated when the transversality conditions are satisfied. The length of delay for which the solutions preserve the stability of the limit cycle is estimated. Finally, through a series of numerical simulations, the theoretical results are formally examined.


Introduction
Cancer is one of the most difficult diseases to treat and is considered one of the leading causes of death. Fighting cancer is important for public health. For this and other economic reasons, extensive research is going on to understand the mechanism involved in the growth of cancer and to predict the impact of this growth on the system [2].
Mathematical modelling is very important in epidemiology because it can provide insight into the key mechanisms that contribute to the spread of disease and suggest control strategies [3]. It is possible to describe a specific disease with mathematical models of differential equations by considering the most important factors which are assumed to be responsible for it and then derive useful information by solving the equations of the model with suitable techniques. Researchers are using mathematical modelling for different diseases like Covid, mumps, rabies, etc., to name a few. Interested readers are referred to the recent works [4][5][6][7].
Mathematical models of tumor growth are widely studied to understand the pathogenesis and predict its future function. One of the most challenging tasks in creating a mathematical model for cancer treatment is calculating biological parameters from data. This task becomes more difficult when the model incorporates many cell types and treatment modalities [8]. Over the past few decades, several experiments and interventions have helped us understand the growth process and its interaction with the immune system.
Many models in the literature represent the dynamics of the systems that vary without time delay [3,[8][9][10][11]. To show the realistic nature of the models, it makes sense to incorporate time delays into the systems. Inclusion of delays in mathematical models makes them more realistic by explaining the effects of disease latency or immunity [12][13][14]. Specifically, the human immune system contains necessary detection systems and weapons to counterattack foreign bodies. This system plays an important role in the prevention of many diseases [15]. When abnormal cells appear in the body, the immune system is activated, and it tries to identify and destroy those abnormal cells. In such a process, the immune system will take some time to provide the right response after discovering the non-self-cells [16][17][18]. As a result, incorporating the delay term makes the model more realistic. Generally, when the system loses stability due to delays, other powerful nonlinear factors such as oscillatory behavior and bifurcation may occur [19][20][21][22]. Accordingly, the analysis of nonlinear dynamics has attracted the attention of many authors in recent years.
Time delay plays an important role in modeling the dynamics of multispecies interactions. Many researchers studied the delay-induced dynamics in predator-prey systems. Beretta and Takeuchi [23] discussed the global stability of the Lotka-Volterra autonomous model with diffusion and delay of time. The delay in the model proposed by Kuznetsov and Taylor was studied by Galach [13], where the delay in time was introduced to achieve better consistency with reality. In 2006, Yafia [24] explored the effect of delays on Kuznetsov's model. Xu and Ma [25] studied the SEIRS epidemic model with a saturation incidence rate and time delay defining the latent period. In 2014, Rihan [26] explored a delay differential model, analyzed it numerically, and demonstrated an effective way of combining chemotherapy with therapeutic immunotherapy. Xu et al. [27] studied a mathematical model with the delay effect on tumor growth due to periodic treatment such that their model was based on reaction-diffusion dynamics and mass conservation law which took into account cell proliferation delays. Yua et al. [28] incorporated delays into the model for activation of the effector cells by helper T-cells, and after analyzing the results, it was shown that delayed immune response could play a major role in regulating tumor growth. In [29], Malinzi studied a delayed model to report delays resulting from virus-cell infection and chemotherapy responses; based on the analysis, they proposed three druginjection methods: constant, single bolus, and periodic treatment. Laaroussi et al. [30] presented a mathematical model of oncolytic virotherapy that included a long-term delay representing multiple periods of the lytic cycle.
In the present study, we investigate a mathematical model that provides an important connection between tumor and immune cells. The model is represented by a system of three differential equations and was originally proposed by De Pillis et al. [1]. To make the model more realistic, a time delay is incorporated into the prey (immune cells) specific growth function due to tumor cells. This consideration has incorporated some novelty into our present investigation. We have analyzed the effect of delay on different dynamical behaviors shown by the model. The paper is structured as follows: Sect. 2 describes the proposed model and the basic assumptions. In Sect. 3, we examine the positive invariance and boundedness of the model solutions. In Sect. 4, the existence of equilibrium points for the system is discussed. Section 5 deals with finding the nature of stability of these equilibrium points, along with the occurrence of Hopf bifurcation. In Sect. 6, some simulation results are presented. In Sect. 7, we have our conclusion.

Model formulation
The model presented in this paper is a modified form of that proposed by De Pillis et al. [1]. The modification is carried out by inducing delay terms in the second term viz. Michaelis-Menten term of the first equation of the system with proper biological justifications. Biochemical reactions involving a single substrate are often assumed to follow Michaelis-Menten kinetics, without regarding the model's underlying assumptions. The change is done to see the impact of delay in the complete system.
The equations representing the system are as follows: where M(t), H(t), and R(t) are the densities of the immune cells, tumor cells, and normal cells, respectively. The first equation of system (1) represents the rate of change of immune cells with respect to time. The first term, namely μ, represents the constant source rate of immune cells. Immune cells are recruited by tumor cells through the Michaelis-Menten term ρMH/(σ + H), where ρ is the rate at which the immune cells grow, and σ represents the steepness of immune response. To describe the time lag by the immune system for developing a suitable response after recognizing the tumor cells, we need to include the effect of time-delay τ into the Michaelis-Menten term [31]. Hence, a discrete-time delay is added to the second term of the first equation of system (1) and, as a result, the term becomes ρM(tτ )H(tτ )/(σ + H(tτ )). The third term, -m 1 MH, represents the kill rate of immune cells due to interaction with tumor cells, while the fourth term indicates that immune cells die off at a rate of d 1 per day.
In the second equation of system (1), 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 interaction with immune cells and normal cells as shown by the terms -m 2 MH and -h 1 HR.
In the third equation of system (1), normal cells are also assumed to grow logistically with an intrinsic growth rate r 2 and the maximum carrying capacity of one. The second term, -h 2 RH, represents the kill rate of the normal cell due to interaction with tumor cells. Below, we have presented the parameter values which we have used for our investigation.
In the next section, we discuss the boundedness and positive invariance of the solutions of the considered model and related subtopics.

Boundedness and positive invariance
By following the standard comparison theory, from the equations of system (1), we get Integration of the above leads to Furthermore, Proceeding as above, we have and similarly, we find Therefore, the numbers of immune-tumor-normal cells are always bounded. Consequently, we get the bounded set So, the solutions of system (1) are bounded. The equations representing system (1) are subject to the following initial conditions: The delay differential system (1) can be written in the vector form aṡ with X = (M, H, R) T ∈ R 3 + and where F ∈ C ∞ (R 3 + ) is defined in the positive quadrant R 3 + and represents a mapping F : The right-hand side of system (3) is locally Lipchitz, meaning that the derivative is bounded and satisfies According to the second lemma in [32], every solution of system (1) with initial conditions ), for all t > 0, remains positive throughout the domain S + , ∀t > 0. Therefore, the solutions of (1) are positively invariant in time t.

Existence of equilibrium points
The equilibrium points of model (1) are found by equating each first derivative to zero. Those are found to be: (i) P 1 (M 1 , 0, R 1 ), the tumor-free equilibrium point. Here, M 1 ≥ 0 and R 1 ≥ 0 which are found to be , the coexisting equilibrium point. Here, M 2 ≥ 0, H 2 ≥ 0, and R 2 ≥ 0 and Note that H 2 can be found by solving the equation The equilibrium point P 2 exists if

Analysis of the stability of the equilibrium points
In this section, we are interested in studying the nature of stability of the tumor-free equilibrium and the coexisting equilibrium of the system.

Local stability analysis of the tumor-free equilibrium point P 1
The Jacobian matrix of system (1) at the equilibrium point P 1 ( μ d 1 , 0, 1) is At the equilibrium point P 1 ( μ d 1 , 0, 1), the eigenvalues of the Jacobian matrix J P 1 are By applying the standard result [3,8], the necessary condition for asymptotic stability of equilibrium point P 1 is found to be r 1 < m 2 μ d 1 + h 1 , and this point will be unstable when r 1 ≥ m 2 μ d 1 + h 1 . Thus, all eigenvalues have negative real parts if r 1 < (m 2 μ/d 1 ) + h 1 for all τ ≥ 0. Accordingly, the equilibrium point P 1 is locally asymptotically stable.

Global stability analysis of the tumor-free equilibrium point P 1
In this section, using a Lyapunov function, we show that the tumor-free equilibrium point P 1 is also globally asymptotically stable when r 1 < (m 2 μ/d 1 ) + h 1 .
To this end, the work [30] is followed to study the dynamics of system (1) when τ ≥ 0.
, and ψ i are continuous functions on the interval [-τ , 0]. Considering a Lyapunov function given as we find that [30] implies that P 1 is globally attractive. This confirms the global asymptotically stability of P 1 when r 1 < (m 2 μ/d 1 ) + h 1 .

Simulation results
To verify the above analytical results, we carry out the following simulations. In doing so, we have considered the parameter values as given in Table 1. Figure 1 shows that when τ ≥ 0 and r 1 < (m 2 μ/d 1 ) + h 1 , the number of tumor cells decreases to zero (approx.) and the growth of immune and normal cells will be stabilized at the tumor-free equilibrium point P 1 . It further shows that when the delay term increases it takes more time to stabilize towards the equilibrium point P 1 .
From the vector field plot (Fig. 2), it is seen that an arbitrary trajectory (in the basin of attraction) starting with different time delays converges to the tumor-free equilibrium point P 1 (indicated by the black dot) indicating that it is globally stable for system (1) as long as the condition r 1 < (m 2 μ/d 1 ) + h 1 is satisfied. Figure 2 further shows (as in Fig. 1) that, as time lag is increased, complete elimination of the tumor cells takes more time. Biologically, we know that unless the tumor cells get eradicated, the immune system remains active, and it eventually drops to its origi-  nal state when all the tumor cells get eliminated. When time lag increases, the immune system will take some time to provide the right response after discovering the non-selfcells, so in this process, immune cells require more time to completely eradicate the tumor cells.
After verifying our analytical results numerically for r 1 < (m 2 μ/d 1 ) + h 1 , we proceed with the numerical verification of the other case.
Case 2: r 1 > m 2 μ d 1 + h 1 . In the case of tumors with a high growth rate (r 1 = 0.4) (which is in the unstable range), Fig. 3 shows that when the immune system takes a long time to provide an adequate response to the tumor cells for recognition, the tumor growth rate gets faster and the system loses its stability and moves away from the coexisting equilibrium point P 2 .
Therefore, we need to consider the effect of delay in the proliferation of the cells induced by the tumor cells when r 1 > (m 2 μ/d 1 ) + h 1 . In the event of a long delay, the system will lose stability since other strong factors can occur such as the periodic oscillatory behavior or bifurcation.
According to the Routh-Hurwitz criterion, all the roots of the equation have negative real parts if and (X 11 + Y 11 )(X 12 + Y 12 ) -(X 13 + Y 13 ) Thus, if the above conditions are satisfied then the coexisting equilibrium point P 2 (M 2 , H 2 , R 2 ) will be stable.
Case (ii): τ > 0. The immune system needs some time to develop a suitable response by recognizing the tumor cells and, therefore, the time delay is introduced into the model.

Occurrence of Hopf bifurcation
We now investigate the occurrence of a Hopf bifurcation when the delay term τ passes through the critical value τ = τ 0 .

Verification of existence of limit cycle through simulation
In this subsection, we verify the above analytical results with numerical simulations for the parameter values listed in Table 1. The singular point P 2 is located at (0.579465, 0.0572258, 0.959124), and we have p 11 = 0.15238135, p 12 = 0.00642814197, and p 13 = -0.0000266. There exists a single pure real root λ 0 = 0.06156 for which we found that τ 0 = 18.91866. Figure 4(a) shows that the equilibrium point P 2 remains a stable spiral for τ < τ 0 = 18.9186. Physically, it means that the growth of the cancerous tumor is slowed down for a while but not eradicated. From Fig. 4(b), it is seen that a limit cycle, i.e., an isolated periodic solution, appears when the value of τ is increased to τ 0 = 18.91866. Physically, it means that with the delay of time, tumorous cells begin to show a dominant nature and, therefore, the immune response is insufficient to reduce the rapid growth of the number of tumor cells. After this point, administration of the drug is approved; otherwise, the tumorous growth becomes faster at the expense of immune and normal cells and can eventually lead the patient to a fatal condition.

Stability of limit cycle: length of time lag estimation
In this section, we investigate the stability of bifurcating periodic solutions and estimate the length of time lag preserving the stability of the period-1 limit cycle. For this, system (1) is first linearized around the interior equilibrium point P 2 (M 2 , H 2 , R 2 ) which gives Applying the Laplace transformation to system (11) leads to Now, using the theory provided by the classical Nyquist criteria [15] and Freedman et al. [33], the equilibrium point P 2 will be asymptotically stable if for the expression the following conditions hold: Re P(iζ 0 ) = 0, Im P(iζ 0 ) = 0, where ζ 0 is the minimal nonnegative root of (12) and (13). From (12), we have Using the inequalities | cos(ζ 0 τ )| ≤ 1 and | sin(ζ 0 τ )| ≤ 1, we get From (15), we find Also, from (13), we have Using (14) and (16) gives Through the inequality (16) for τ = 0, (17) becomes Now, the first and second terms of the left-hand side of (18) can be respectively written as Therefore, from (18), we find where Now, expression (18) gives Therefore, the period-1 limit cycle arising in the system preserves stability around the equilibrium P 2 for the maximum length of the time lag τ + .

Other simulation results
In this section, we demonstrate a series of numerical simulations associated with the theoretical results obtained in the previous sections. In particular, we are interested in the scenarios where the induced time delay changes the stability, i.e., situations where, for the same set of parameters, a stable steady state becomes unstable with a change in the delays or vice versa. The parameter values have been selected from Table 1.
Case 1: r 1 < m 2 μ d 1 + h 1 . It has already been shown in the figures of Sect. 5.3. Case 2: r 1 ≥ m 2 μ d 1 + h 1 . From Figs. 5 and 6, it is seen that when discrete-time delay τ increases, the growth of tumor cells also increases, and the system loses its stability and moves away from the coexisting equilibrium point P 2 . When the time lag τ < τ 0 , the equilibrium point P 2 is asymptotically stable and the solutions converge to a coexisting equilibrium point. The solutions of the system for τ < τ 0 can be seen in Figs. 5(a) and 6(a) where the system has damped oscillations and the coexisting equilibrium point is locally asymptotically stable. Figures 5(b) and 6(b) show a stable limit cycle with time delay τ 0 = 18.91866, i.e., high competition among the three cell populations, namely tumor, normal, and immune cells.
Furthermore, from Figs. 5(c) and 6(c), it is observed that when the time lag τ > τ 0 , the system is unstable. Determination of the bifurcation point is important for the control of tumor cell populations. After the bifurcation point, system (1) will exhibit unstable behavior, which in turn leads to the uncontrolled growth of the tumor (Fig. 5(c)). It was found that immune and normal cells fought against the tumor cells but failed to stop the growth of the tumor cells.
In Figs. 5 and 6, it can be observed that the proliferation of tumor cells depends on the time taken by the immune system's response to the tumor cells. To balance the effect of delay, it is necessary to increase the growth rate of immune cells and the strength of normal cells in their competition against tumor cells. Hence, in the biological sense, the results have practical significance in terms of determining the amount of drugs required to eliminate the tumor, otherwise, the tumor may grow in a periodic way putting the health of the patient in danger.

Conclusion
In this paper, we studied the dynamical behavior of a nonlinear model which was originally proposed by De Pillis et al. by inducing a delay term in the immune system-tumor interaction term. It was done to make the model more realistic. Mathematical analysis of the model demonstrated that, when r 1 < (m 2 μ/d 1 ) + h 1 , the tumor cells can be eradicated by the combined effort of the immune system and normal cells without the application of any drug. When r 1 > (m 2 μ/d 1 ) + h 1 and the immune system takes a long time to recognize the tumor cells to give an adequate response (i.e., delay term is large), the tumor growth rate is faster and the system loses its stability and moves away from the tumor-free equilibrium point P 1 and, as a result, the immune-normal cells fail to eradicate the tumor load.
To understand more details about the dynamical behavior of the system, we explored the coexisting equilibrium point when r 1 ≥ (m 2 μ/d 1 ) + h 1 . After analysis, it was concluded that the Hopf bifurcation occurs when the time delay increases to τ = τ 0 . When the delay term is τ < τ 0 , the coexisting equilibrium is locally asymptotically stable at the coexisting equilibrium point P 2 and becomes unstable when τ > τ 0 . In such a situation, to balance the effect of delay, it is necessary to increase the growth rate of immune and normal cells in their competition against tumor cells by employing some drugs.
All of our conclusions are based on the theoretical investigation and some numerical simulations supporting the results. Experimental or clinical verifications will give proper feedback whether the assumptions considered in the formulation of the model and the parameter values used are correct or if they need some modifications. Moreover, the numerical simulations involved in the paper were carried out by using the inbuilt software of Mathematica. The order of convergence of the methods employed, the CPU time required, etc., were not addressed as we were interested in studying only the dynamical behaviors shown by the system. Readers who are interested in those aspects and also in different methodologies to carry out the investigation of different problems in epidemiology are referred to the papers [34][35][36][37][38][39][40][41][42][43].