Mathematical analysis of a cancer model with time-delay in tumor-immune interaction and stimulation processes

In this study, we discuss a cancer model considering discrete time-delay in tumor-immune interaction and stimulation processes. This study aims to analyze and observe the dynamics of the model along with variation of vital parameters and the delay effect on anti-tumor immune responses. We obtain sufficient conditions for the existence of equilibrium points and their stability. Existence of Hopf bifurcation at co-axial equilibrium is investigated. The stability of bifurcating periodic solutions is discussed, and the time length for which the solutions preserve the stability is estimated. Furthermore, we have derived the conditions for the direction of bifurcating periodic solutions. Theoretically, it was observed that the system undergoes different states if we vary the system’s parameters. Some numerical simulations are presented to verify the obtained mathematical results.

Caputo-Fabrizio fractional-order derivative. They calculated the model results for different fractional order and compared those with the results of the integer-order model. Aydogan et al. [3] examined a Caputo-Fabrizio fractional-order mathematical model of Rabies disease with the use of the Laplace Adomian decomposition method. A box model of mumps-induced hearing loss in children using the Caputo-Fabrizio fractional-order derivative that preserves the system's historical memory was investigated in [4]. The authors also determined the optimal control problem for the proposed model considering treatment as a control parameter to reduce the infected population. By examining the sensitivity of basic reproduction numbers to each of the model parameters, they showed that the basic reproduction number increases with the increase of disease transmission rate and daily birth rate. Also, it reduces with an increase in the recovery rate, normal mortality rate.
The tool of mathematical modeling is widely used by researchers in the field of cancer modeling too. Several significant works [5][6][7][8][9] through mathematical modeling have been done to understand the response of the immune system with tumor. In [5], the authors analyzed a cancer model by considering the interactions between cancer cells, tumor angiogenesis and endothelial cells. Lopez et al. [6] estimated the decay rate of tumor cells in the presence of immune response and set a threshold value for which the immune system can eradicate the tumor. Dong et al. [7] explored the effect of CD4 + T cells in a tumorimmune system incorporated with adoptive cellular immunotherapy (ACI) and suggested that CD4 + T cells play a crucial role in the tumor eradication process under ACI therapy. Arlotti et al. [8] proposed a bilinear model of integro-differential equations to describe the dynamics of cellular interaction between tumor and immune cells. By introducing two phases, namely, the Inter-phase and M-phase at which cells are generally divided, Awang et al. [9] newly described the interactions between tumor cells and immune system. Zeng and Ma [10] analyzed a deterministic tumor-immune model under the Allee effect. This Allee effect disturbed the growth and reproduction of tumor-immune cells. They found a range of Allee threshold values for which the system can be stabilized.
The cytotoxic T cells are primarily responsible for tumor suppression and are found in all tissues in the body. Kuznetzov et al. [11] described the response of CTL cells to the growth of an immunogenic tumor. Quinonez et al. [12] examined a mathematical model to show the immune response in the presence of synthetic tumor vaccines to mitigate developing cancer. Li and Li [13] modified the Kuznetsov et al. [11] and Galach model [14] to a stochastic one by perturbing environmental noise. They found that the reduced rate of tumors in their model is faster than in the earlier models. Also, their result reveals that environmental noise is favorable for the extinction of tumor cells under immune surveillance, and noise is ineffective when the immune system's ability is strong enough. De Pillis et al. [15] explored the role of natural killer (NK) and CD8 + T cells in the immune-tumor interaction mechanism and tumor surveillance through a mathematical model. Another mathematical model was developed by Dritschel et al. [16] to investigate the anti-tumor immune response of helper and cytotoxic T cells. They concluded that the tumor growth could be reduced if we use treatments like IL-2 therapy, adoptive T cell therapy which boosts the immune system, and antibody therapy which blocks tumor-induced immunosuppression. Recently, Pang et al. [17] proposed a modified model of [18,19] to reflect the clinical phenomena of anti-tumor immune responses and observed that the tumor cells show initial exponential growth to the final stable position at zero depending upon the flow rate of mature immune cells.
In many biological complex systems, time delay plays a vital role. The introduction of time delay forces a system to depend not only on the present state but also on the past state. In cancer modeling, the time delay can describe the required time for cell differentiation, cell proliferation, the response of one cell to other cells, etc. Banerjee and Sarkar [20] studied a delay-induced tumor-immune model to control the growth of malignant cells. They varied parameters, analyzed the model, and observed that Hopf bifurcation occurs for delay term as the bifurcation parameter. Rihan et al. [21] considered a family of differential models to explore the effects of ACI therapy to control tumor burden. In a study of Bi and Xiao [22], they illustrated the bifurcation analysis of a modified version of the Kuznetsov model [11] with the introduction of two-time delays for immune response into the model. Khajanchi [23] and Khajanchi et al. [24] analyzed in detail the influence of time delay on the chaotic dynamics of tumor-immune interaction model [25]. In [23,24], the authors presented that in the presence of delay, the model shows long-term chaotic behavior with Hopf bifurcations. They also estimated the length of the time delay to preserve stability and direction of Hopf bifurcation. Ghosh et al. [26] described the interaction between tumor cells and micro-environment immune and host cells with the use of two time delays, one for immune interaction with tumor and the other for immune action on the tumor. Dong et al. [27] converted their previously proposed model [7] to a delayed one with the use of two delays, namely the immune activation delay for ECs and immune activation delay for HTCs. Their results showed that the unstable equilibrium goes to a stable position for the activation delay of HTCs. In the work [28], the authors modified the model proposed by Dong et al. [7] with the use of one delay term for the activation of ECs by HTCs. This delay term induced two effects on the model. The first effect is that stability switches to instability, and the second stabilizes tumor-presence equilibrium. Further modification of the model [7] was carried out by Das et al. [29] with the use of distributed and discrete-time delay. They observed that uniform activation of helper T-cells can help in ECs stimulation and tumor control. Considering a reaction-diffusion system, including time delay under the Neumann boundary conditions, Kayan et al. [30] modified the model [11] which described tumor-immune competitions. They studied the Hopf bifurcation analysis and found that the effect of diffusion of tumor-immune interaction can significantly change the dynamics of the model.
In this article, we investigate the model proposed by Pang et al. [17] introducing delay term for anti-tumor immune responses of matured T lymphocytes to destroy tumor cells. In Sect. 2, we formulate our proposed model. A qualitative analysis is done in Sect. 3. In Sect. 4, conditions of Hopf bifurcation of the system are analyzed. Section 5 and Sect. 6 deal with the stability of the limit cycle and direction and stability of Hopf bifurcation. Numerical examples, discussion, and conclusions are carried out in respective Sects. 7 and 8.

The model
Mathematical models can give a better insight into tumor-immune interaction. In the above literature, we have found that each model is proposed to understand the tumor cells' mechanism and reduce the tumor burden in the body. T lymphocytes are the most important cells in the immune system, capable of killing the tumor cells through kinetic processes. However, it is important to note that the tumor cells can also compete with T lymphocytes and make them functionally inactive. Also, tumor cells secrete cytokines, which are responsible for tumor cells proliferation. As a result, eradicating all tumor cells becomes difficult for T lymphocytes, and a time delay occurs for the deactivation of tumor cells by T lymphocytes. This time delay is regarded as an interaction and stimulation delay of the tumor-immune system. We account for this time delay, the model of Pang et al. [17] is thus modified into a delay system as follows: where L 1 (τ ), L 2 (τ ), and T(τ ) are densities of immature T lymphocytes, mature T lymphocytes, and tumor cells at any time τ respectively. In the first equation of (2.1), the first term μ describes the fixed production rate of immature T lymphocytes by the body in the absence of tumor cells. The second term λ 1 L 1 (τ ) is used for describing the transformation rate of immature T lymphocytes to mature T lymphocytes. The third term α 1 describes the recruitment term of anti-tumor immune response, where α 1 is the maximum recruitment rate and η is the half-saturation constant. Here, represents the discrete-time delay factor added for interaction and stimulation delay of the tumor-immune system. The second equation describes the dynamics of mature T lymphocytes where α 3 is the inactivation rate of T lymphocytes. The third equation designates the rate of change of tumor cells in which tumor cells can grow exponentially in the absence of immune response, where λ 2 is the exponential tumor growth rate. The term α 2 T(τ -)L 2 (τ -) describes the interaction between tumor and mature T lymphocytes, where α 2 is the rate of tumor cells killed by the mature T lymphocytes.

Basic properties
The following proposition establishes the nonnegativity of the solutions of (2.2) with (2.3). Proof System (2.2) can be written aṡ is locally Lipschitz, and hence the derivatives are bounded, satisfy the conditions According to the lemma by Yang et al. [31], every solution of system (2.2) with initial values (2.3), φ i (t) ∈ R 3 + , say, X(t) = X[t; X(0)], ∀t > 0, that is, it remains positive throughout the region R 3 + , ∀t > 0.
System (2.2) without time delay was discussed in [17]. Since inducing time delay does not affect the existence conditions for equilibria of the system, the system has three equilibria with E 0 (0, 0, 0) as a trivial equilibrium, E 1 (0, a 3 a 2 , 0) as a tumor-free one, which always exist, and E 2 ( a 2 a 4 -a 3 a 1 , a 4 , a 2 a 4 -a 3 a 3 +a 1 a 4 -a 2 a 4 ) as a co-axial equilibrium, which exists for max{0, (a 2a 1 a 4 )} < a 3 < a 2 a 4 .
III. Now, we will investigate the dynamical behavior of the system (2.2) around the co-axial fixed point E 2 (x = a 2 a 4 -a 3 a 1 ,ŷ = a 4 ,ẑ = a 2 a 4 -a 3 a 3 +a 1 a 4 -a 2 a 4 ) under the influence of discrete time lag for both the cases of existence of E 2 . In this case, the eigenvalues of matrix (3.3) can be found from the following equation: For the case of no time lag ( = 0), equation (3.4) becomes By the use of Routh-Hurwitz criterion to (3.5), E 2 is asymptotically stable if the following conditions are satisfied: Now, we shall analyze the dynamical behavior of the system (2.2) with time lag = 0. For this, we assume that there exists a purely imaginary root for (3.4). Hence, we substitute λ = im(m > 0) into (3.4) and, separating the real and imaginary parts, Using the method of cross-multiplication, we solve both the equations of (3.7) and get By squaring and adding both sides of both the equations of (3.7), we get where The equation (3.9) will have a positive root if Suppose that m 0 is a unique non-negative root of the equation (3.9) such that (3.4) has a pair of purely imaginary roots of the form ±im 0 . Then from the equation (3.8) the time lag k corresponding to m 0 is Therefore, the co-axial equilibrium E 2 is locally asymptotically stable under the conditions (3.6) for k = 0. So, this point will also remain stable for k < 0 , where k = * at k = 0 by Butler's lemma [24]. This suggest that for k > 0 the co-axial equilibrium E 2 is unstable. This implies that the tumor cells can proliferate faster if the interaction time delay crosses a given critical value and the system loses its stability at E 2 . The rest of the work which is discussed in Sect. 4, Sect. 5, and Sect. 6 is inspired and followed by the previous works [20,23,24,32,33].

Analysis of Hopf bifurcation
As the equation (3.4) has a complex roots of the form λ = im 0 , it implies that the system (2.2) may undergo a Hopf bifurcation at = k and around the equilibrium E 2 . Here, we establish a condition for which the system (2.2) undergoes a Hopf bifurcation by using Lemma 1 [17]. For this, first we need to verify the transversality condition d(Re λ) Differentiating (3.4) with respect to gives (4.1) Simplifying equation (4.1), we have Therefore, with the increase of the value of k the direction of motion of λ is given by

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 period-1 limit cycle. ,ŷ = a 4 ,ẑ = a 2 a 4 -a 3 a 3 +a 1 a 4 -a 2 a 4 ), which gives uṡ .

Direction and stability of Hopf bifurcation
In the previous sections, we have discussed the conditions for which the periodic solutions of the system (2.2) bifurcate from co-axial equilibrium E 2 at the critical values of k via the Hopf bifurcation. In this section, we will analyze the direction, stability, and periods of the periodic solutions of the system (2.2) using normal theory and the center manifold theorem developed by Hassard et al. [35]. Let , z(t) = u 3 ( t), and = 0 + μ, where 0 is defined by (3.11), and μ ∈ R. The system (2.2) can be written as a functional differential equation in C = C([-1, 0], R 3 ) as where X(t) = (x(t), y(t), z(t)) T ∈ R 3 , and L μ : C → R 3 , f : R × C → R 3 are given, respectively, as follows: for (t) = ( 1 (t), 2 (t), 3  where HOT → higher order terms. By the Riesz representation theorem, there exists a matrix function η(θ , μ) of bounded variation for θ ∈ [-1, 0] such that For the Dirac delta function δ, choose and for ∈ C 1 ([-1, 0], R 3 ), define and (6.9) Hence, system (6.1) is equivalent to the operator equation (ξθ ) dη(θ ) (ξ ) dξ , (6.12) where η(θ ) = η(θ , 0). Then A(0) and A * are adjoint operators. We already assume that ±ιω 0 k are the eigenvalues of A(0). Hence, the eigenvalues of A * are ∓ιω 0 k . We need to compute the eigenvectors of A(0) and A * corresponding to the eigenvalues ιω 0 k and -ιω 0 k respectively. Suppose that v(θ ) = (1, v 1 , v 2 ) T e ιω 0 k θ is the eigenvector of A(0) corresponding to ιω 0 k . Then A(0)v(0) = ιω 0 k v(0). It follows from the definition of A(0) and (6.3), (6.6) and (6.7) that In a similar manner, we can obtain the eigenvector v * (s In order to guarantee v * (s), v(θ ) = 1, we need to determine the expression of D v * (s), v(θ ) Therefore we can choose D as On the other hand, due to adjoint property, we can write , A = A * , .
We have Therefore, v * ,v = 0. Now, we will compute the coordinates describing the center manifold C 0 at μ = 0. Let X t be the solution of (6.10) when μ = 0. We define where z andz are the local coordinates for the center manifold C 0 in the direction ofv * and v * . Note that W is real if x t is real. Here, we are only interested in real solutions. From (6.13), we have For the solution x t ∈ C 0 in (6.10), since μ = 0, hencė (6.15) From (6.13) and (6.14), it follows that Thus, where X 1t (θ ) = ze ιω 0 k θ +ze -ιω 0 k θ + W 1 20 (θ ) Using these values and from (6.3) it follows that Putting the values of X 2t (-1), X 3t (-1), and X 2 3t (-1), computing the above expressions (6.16), and comparing the coefficients of z 2 , zz,z 2 , and z 2z with (6.15), we have and W 20 (θ ) = ιg 20 ω 0 k q(0)e ιω 0 k θ + ιḡ 02 3ω 0 kq (0)e -ιω 0 k θ + E 1 e 2ιω 0 k θ , 2 ) are a constant vector in R 3 satisfying the following equations: Furthermore, we can compute g 21 with respect to parameters and delay. Hence, from the above analysis we can conclude that in order to find each g ij we have to use the parameters and delay in system (2.2). Thus we can compute the following values: Based on our analysis, by the result of Hassard et al. [35], we have the following result.

Numerical simulations
Through our above analysis, we have gained an analytical understanding of the possible dynamics of our proposed nonlinear delay differential equation model (2.2). In this section, bifurcation analysis and parameter sensitivity will be discussed. The delay model (2.2) showed that the system exhibited oscillations. Here, we shall discuss the effects of varying the parameter a 1 on generating the oscillatory behavior. In Fig. 1, X is written as a function of a 1 . The curves (equilibrium branch) having red and black colors represent the stable and unstable steady state branches, respectively. Further, the green curve denotes the maximum and minimum of the limit cycle. The equilibrium branch loses its stability due to the appearance of a Hopf bifurcation point. In the unstable branch, the system (2.2) exhibits oscillations, and from the figure, it can be observed that the amplitude of these oscillations decreases as a 1 increases. Similar effects of a 2 on the steady states of X are discussed in Fig. 2. The red, black curves represent the stable and unstable steady state branches, whereas the green circle denotes the stable limit cycles. For low and high values of a 2 , the system (2.2) has a stable equilibrium point. The system (2.2), however, does not converge to a steady-state when a 1 = 0 or a 2 = 0. These are shown by arrows in Figs. 1 and 2. The reason for this is because the equilibrium expression E 1 is divided by a 2 and E 2 is divided by a 1 . Hence, E 1 , E 2 → ∞.
After discussing the effects of a 1 and a 2 on generating limit cycle oscillations, we need to investigate how the interplay between the two parameters a 1 and a 2 can alter these os-  cillations. To address this question, we shall track the two parameters on a two-parameter plot. Figure 3 shows a two-parameter plot. Within the cusp, the system (2.2) exhibits oscillations. Outside the cusp, the model (2.2) only has a single stable steady-state. In the absence of a 1 or a 2 , the system cannot generate limit cycle oscillations. For the model (2.2) to generate an oscillatory response, it must have a "well" balance between a 1 and a 2 . If this balance is biased, the oscillatory response cannot be achieved.
When selecting values of a 1 and a 2 from the region of oscillations of the figure, the model (2.2) shows oscillation. We shall focus on the effects of a 3 on these oscillations. In other words, how does it alter this response? Figure 4 shows a two-parameter bifurcation   Fig. 1 diagram, and it shows that if either of these two parameters increases beyond a critical level, the oscillations will be destroyed. Additionally, in the absence of a 3 , the model (2.2) can generate an oscillatory behavior.
Here, we shall focus numerically on the effects of the delay term on the limit cycle oscillations. Figure 5 illustrates the relationship between a 1 and on generating limit cycle oscillations. In this Figure, in the absence of the delay term, the system (2.2) can generate oscillations if a 1 is increased beyond the blue vertical line (the Hopf locus).
We shall consider the role of a 2 and the delay term on generating oscillations and this can observe in Fig. 6. The system (2.2) can generate an oscillatory behavior even without the   Fig. 1 delay term as shown in Pang et al. model [17]. Figure 5 and Fig. 6 suggests that the system (2.2) requires greater a 2 value than a 1 value in order to be able to generate oscillations. Figure 7 highlights the contribution of the parameters a 3 and on destabilizing the stable steady state branch. For the model (2.2) to generate oscillations, it requires greater a 3 values than a 1 and a 2 . Thus, the system (2.2) can cross the Hopf locus and generate an oscillatory behavior. Figure 8 shows that the model (2.2) requires greater a 4 values in order to exhibit oscillations that a 1 , a 2 , and a 3 . However, suitable parameter values often give a meaningful biological scenario of the system (2.2). Therefore, we perform some simulation works for a better understanding of our analytical treatment. We consider different values of the parameters and the delay factor ( ) to observe biologically plausible different dynamical scenarios of the model (2.2), enough to merit the mathematical study.
Choosing a 1 = 0.3, a 2 = 0.6, a 3 = 4.0, a 4 = 10.0, then a 2 > a 1 and (a 2a 1 )a 4 < a 3 < a 2 a 4 , which implies the existence of equilibria E 1 and E 2 . Using the conditions of local stability, we get that both the equilibria E 1 and E 2 are unstable. Also, there occurs a periodic solution at E 2 . Figures 9 and 10 show the oscillating behavior as well as the periodic solutions for the system (2.2). Existence of periodic solutions is relevant in cancer models. It implies that the tumor levels may oscillate around a fixed point even in absence of any treatment. Such a phenomenon, which is known as Jeff 's phenomenon, has been observed clinically. We observe that is beneficial for tumor cells. We observe no stability switch in the system (2.2) as the delay factor increases.
Choosing a 1 = 0.4, a 2 = 0.6, a 3 = 3.5, a 4 = 7.0, then a 2 > a 1 and (a 2a 1 )a 4 < a 3 < a 2 a 4 , which implies the existence of equilibria E 1 and E 2 . Using the conditions of local stability,  we get tumor-free equilibrium E 1 as unstable and co-axial equilibrium E 2 as stable in nature. From Figs. 11 and 12, we observe a stability switch in the system (2.2) as the delay factor crosses a threshold.  Choosing a 1 = 0.4, a 2 = 0.6, a 3 = 3.5, a 4 = 7.0, then a 2 > a 1 and 0 < a 3 < (a 2a 1 )a 4 , which indicates the existence of tumor-free equilibrium E 1 . At this equilibrium the system (2.2) shows unstable behavior. Figures 13 and 14 indicate that the number of tumor cells increases without restriction, which is in accordance with the immune escape phenomena of tumor which is observed clinically.
Choosing a 1 = 0.3, a 2 = 0.6, a 3 = 0.5, a 4 = 0.7, then a 3 > a 2 a 4 , which suggests the existence of the stable tumor-free equilibrium E 1 . Figures 15 and 16 clarify that the tumor-free equilibrium is asymptotically stable, which is in accordance with the spontaneous tumor regression phenomena observed clinically.

Conclusion
In this paper, we aimed to investigate the dynamical behavior of the modified Pang et al. model [17]. In their model [17], it was shown that with the increase of normal flow rate of mature immune cells, the system exhibits different states such as tumor dormant, periodic tumor oscillation, immune escape of tumor, and so on. However, the effects of the delay term on the oscillatory behavior were not considered in their model. Therefore, a delay term was included in our model, and we investigated the system behavior with varying system parameters. As a result, the modified model showed that the system (2.2) could generate an oscillatory response even with a delay term. Moreover, it illustrated that these oscillations were persistent and could not be destroyed by the additional delay term. Our bifurcation analysis and numerical simulations revealed that a "careful" selection of the model's parameters must be obtained so that the stable steady-state loses its stability. We showed that the delay term was not necessary to generate oscillations because our model can generate these oscillations even without the delay term. From the bifurcation analysis, in addition to Figs. 9 and 10, it can be shown that neither affects generating of oscillations nor the amplitude of these oscillations. However, varying other parameters such as but not limited to a 1 , a 2 leads to the stabilization of the unstable equilibrium point. A set of realistic parameter values gives us a better insight into the model, which we leave as our future work.
In our entire discussion, the major goal was to have a dynamical analysis of the considered model with the incorporation of a delay term. The numerical results were obtained by applying standard MATLAB software. The aspects of numerical stability, CPU time, minimum error, etc. of the adopted numerical techniques were not investigated. The interested readers and researchers are referred to [36,37] for such kinds of investigations.