key: cord-0932802-gmjmv9dx
authors: Dehingia, Kaushik; Sarmah, Hemanta Kumar; Alharbi, Yamen; Hosseini, Kamyar
title: Mathematical analysis of a cancer model with time-delay in tumor-immune interaction and stimulation processes
date: 2021-10-26
journal: Adv Differ Equ
DOI: 10.1186/s13662-021-03621-4
sha: 2630bd49d4a1e0cf12756ec9ccda676cd7e07245
doc_id: 932802
cord_uid: gmjmv9dx

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.

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.

For the sake of discussion of model (2.1), we substitute μ = λ 0 L 0 and use nondimensional variables and parameters as

The normalized model of (2.1) is

with initial history:

with φ i ≥ 0, ∀i = 1, 2, 3 for θ ∈ [-, 0], where φ i (θ ) ∈ R 3 + are the continuous functions on [ , 0) that may display jumps at θ = 0.

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 .

• For a 3 > a 2 a 4 , the tumor-free equilibrium E 1 (0, a 3 a 2 , 0) exists. • There exists only one equilibrium E 1 for a 2 > a 1 and 0 < a 3 < (a 2a 1 )a 4 . • For a 2 < a 1 and 0 < a 3 < a 2 a 4 , two equilibria E 1 and E 2 exist. • For a 2 > a 1 and (a 2a 1 )a 4 < a 3 < a 2 a 4 , two equilibria E 1 and E 2 exist.

In order to check the local stability at each of the equilibrium points, we calculate the following Jacobian matrix for the system (2.2):

I. The "no living cell" fixed point E 0 (0, 0, 0) is stable if a 4 < 0. As corresponding to E 0 , the eigenvalues of the matrix (3.3) are λ 0,1 = -1 (< 0), λ 0,2 = -a 2 (< 0), and λ 0,3 = a 4 .

II. The eigenvalues of the matrix (3.3) are λ 1,1 = -1, λ 1,2 = -a 2 , and λ 1,3 = -a 3 -a 2 a 4 a 2

corresponding to the tumor-free fixed point E 1 (0, a 3 a 2 , 0). So, E 1 is stable if a 3 > a 2 a 4 , otherwise unstable. For the case of a 2 > a 1 and 0 < a 3 < (a 2a 1 )a 4 , the fixed point E 1 is a saddle-focus point with two negative eigenvalues and one positive eigenvalue.

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] .

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 

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ṡ

.

Using Laplace transformation into (5.1), we have

and L x (ω), L y (ω), and L z (ω) are the Laplace transformations of x(t), y(t), and z(t) respectively. Now, by using the theory provided by Freedman et al. [34] and the classical Nyquist criteria, the equilibrium point E 2 is asymptotically stable if, for the equation 

Using (5.6), equation (5.9) becomes

Using the inequality ζ 2 0 < (a 2a 2 a 4a 4 + a 2ŷ +ŷ -a 1ẑ 1+ẑ ) for = 0, the above equation has the form a 1 a 4ẑ 1 +ẑ + a 2ŷ -ŷζ 2 0 -(1 + a 2a 4 ) a 2ŷ +ŷ -a 1ẑ 1 +ẑ cos(ζ 0 ) -1 + a 2ŷ +ŷ -a 1ẑ 1 +ẑ (1 + a 2a 4 )ŷ ζ 0 + (1 + a 2a 4 )( a 1 a 4ẑ 1+ẑ + a 2ŷ ) ζ 0 sin(ζ 0 ) < (1 + a 2a 4 +ŷ) a 2a 2 a 4a 4 + a 2ŷ +ŷ -a 1ẑ 1 +ẑ --a 2 a 4 + a 1 a 4ẑ 1 +ẑ + a 2ŷ .

Now the first term and the second term of the L.H.S. of (5.11) can be written respectively as a 1 a 4ẑ 1 +ẑ + a 2ŷ -ŷζ 2 0 -(1 + a 2a 4 ) a 2ŷ +ŷ - 

where

a 1 a 4ẑ 1 +ẑ + a 2ŷ ,

Now, equation (5.12) gives + = 1 2χ 1 -χ 2 + χ 2 2 + 4χ 1 χ 3 for 0 ≤ ≤ + . (5.13) Therefore, the system (2.2) preserves the stability around the equilibrium E 2 for the maximum length of the time lag + with period-1 limit cycle.

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.

The sign of μ 2 , β 2 , and T 2 determines the directions of Hopf bifurcations, stability of the bifurcating periodic solutions, and the period of bifurcating periodic solutions respectively for = k . In view of (6.22), the following results hold for system (2.2): (a) If μ 2 < 0 (μ 2 > 0), the Hopf bifurcation is subcritical (supercritical). (b) If β 2 > 0 (β 2 < 0), the bifurcation periodic solutions are unstable (stable). (c) If T 2 < 0 (T 2 > 0), the period of the bifurcated periodic solution decreases (increases).

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.

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.

SEIR epidemic model for COVID-19 transmission by Caputo derivative of fractional order

A study on the AH1N1/09 influenza transmission model with the fractional Caputo-Fabrizio derivative

On the mathematical model of Rabies by using the fractional Caputo-Fabrizio derivative

A theoretical study of the Caputo-Fabrizio fractional modeling for hearing loss due to Mumps virus with optimal control

Stability of a mathematical model of tumour-induced angiogenesis

Bifurcation analysis and nonlinear decay of a tumor in the presence of an immune response

Mathematical modeling on helper T cells in a tumor immune system

A kinetic model of tumor/immune system cellular interactions

Analysis of tumor populations and immune system interaction model

Dynamic analysis of a tumor-immune system under Allee effect

Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis

A mathematical investigation on tumor-immune dynamics: the impact of vaccines on the immune response

Long time behavior of a tumor-immune system competition model perturbed by environmental noise

Dynamics of the tumor-immune system competition-the effect of time delay

A validated mathematical model of cell-mediated immune response to tumor growth

A mathematical model of cytotoxic and helper T cell interactions in a tumour microenvironment

Mathematical modelling and dynamic analysis of anti-tumor immune response

Immune surveillance and neoplasia-I: a minimal mathematical model

Stable periodic oscillations in a two-stage cancer model of tumor and immune system interactions

Delay-induced model for tumor-immune interaction and control of malignant tumor growth

Qualitative and computational analysis of a mathematical model for tumor-immune interactions

Bifurcations of tumor-immune competition systems with delay

Chaotic dynamics of a delayed tumor-immune interaction model

The influence of time delay in a chaotic cancer model

The dynamics of an optimally controlled tumor model: a case study

How tumor growth can be influenced by delayed interactions between cancer cells and the microenvironment?

Dynamics in a tumor immune system with time delays

Dual role of delay effects in a tumour-immune system

Effects of delayed immune-activation in the dynamics of tumor-immune interactions

Bifurcation analysis of a modified tumor-immune system interaction model involving time delay

Permanence and positive periodic solution for the single-species nonautonomous delay diffusive model

On the stability and Hopf bifurcation of a predator-prey model

Hopf bifurcation analysis of a delayed SEIR epidemic model with infectious force in latent and infected period

The trade-off between mutual interference and time lags in predator-prey systems

Theory of Hopf Bifurcation

Nonlinear waves propagation and stability analysis for planar waves at far field using quintic B-spline collocation method

A computational approach for solving time fractional differential equation via spline functions

The authors of this paper thank both the anonymous reviewers on the basis of whose reports and suggestions the improvement of the paper was possible.

There is no funding for this work.

Not applicable.

The authors declare that they have no competing interests.

All authors jointly worked on the results, and they read and approved the final manuscript. 

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.