EFFECTS OF DELAYED IMMUNE-ACTIVATION IN THE DYNAMICS OF TUMOR-IMMUNE INTERACTIONS

. This article presents the impact of distributed and discrete delays that emerge in the formulation of a mathematical model of the human immunological system describing the interactions of eﬀector cells (ECs), tumor cells (TCs) and helper T-cells (HTCs). We investigate the stability of equilibria and the commencement of sustained oscillations after Hopf-bifurcation. Moreover, based on the center manifold theorem and normal form theory, the expression for direction and stability of Hopf-bifurcation occurring at tumor presence equilibrium point of the system has been derived explicitly. The eﬀect of distributed delay involved in immune-activation on the system dynamics of the tumor is demonstrated. Numerical simulations are also illustrated for elucidating the change of dynamic behavior by varying system parameters.


P. DAS ET AL.
A summary of significant works of tumor-immune system(T-IS) interaction can be obtained in [2,3].In 1994 Kuznetsov et al. [30] proposed and studied a mathematical model of a cell-mediated tumor population.The main difference from other models is that it takes into account the percolation of tumor cells (TCs) by effector cells (ECs) as well as inactivation of effector cells.Generalized Kuznetsov model was studied by Kirshner and Panetta [29] to describe the dynamics of interactions between the TCs, ECs, and IL-2.They incorporated adoptive cellular immunotherapy to portray short term oscillations in tumor growth level and also long term tumor relapse.Zhange et al. proposed a model which generates periodic and chaotic oscillations in T-IS interaction [7].De Pills and Radunskaya suggested a model regenerating asynchronous tumor-immune resistance through drug therapy known as Jeff's phenomena [16].Forys studied Marchuk's model of the immune system for chronic state [22].Interaction between cancerous cell and micro-environment introducing complex phenomena was discussed in [15,24].Banerjee et al. emphasized on delay-induced malignant tumor growth and control [4].Chaotic cancer model also revealed a new approach of research in tumor growth before and after treatment [14].
Incorporating time-lag in T-IS interaction for molecule production, proliferation, differentiation and transportation of cells etc. is essential in mathematical model [5,33] and it influences the dynamics of the physiological system [13,32].The study on dynamics of T-IS interaction with discrete delay has been noteworthy interest for a long time [6,17,19,21,23,34,36].Ruan et al. [7] investigated tumor-immune interactions with different discrete delay and established asymptotic stability below some threshold value of delay and Hopf bifurcation at its critical value.Zhang et al. studied a model with three delays which results in periodic and chaotic oscillations in tumor-immune system [7].
Recently, distributed delays are considered in cancer-immune system [32].Some of the researchers studied how the distributed delays affect the dynamics of the system differently from discrete delay [20].Signal transmission during cellular phenomena can be described by a sequence of linear equation with distributed delay.Immune response is no longer instantaneous and is continuous in nature.With distributed delay, it is easier to shorten the oscillation than with a discrete delay due to negative feedback on the system [11].A non-negative bounded delay kernel K(.)(say) is defined on [0, ∞) is considered which reflects the influence of past states on current dynamics [12].Dong et al. proposed a model on helper T-cells in immune system [17].Again Dong et al. introduced delay kernel as immune-activation in HTCs without taking into account the stimulation of TCs in the presence of ECs along with discrete delays [38].Helper T-cells do not kill the cancer cell, instead, they continuously help activated immune-effector cell such as cytotoxic T-cells to kill target cell and also help stimulated B-lymphocyte to secrete antibodies and macrophages to kill ingested microbes.Furthermore, effector cells cannot defend the infected target cells without helper-T cells.
Encouraged by the studies briefly outlined above, in this paper, we explore the work of Dong et al. [18] by introducing immune-activation delay kernel as an infinite distributed delay for enhancing the continuous stimulation of HTCs.We are interested not only in investigating the model system describing possible phenomena quite exhaustively at the cell level but also in tissue level.Our aim is to study the effect of interaction rate of ECs with TCs, interaction rate of HTCs with TCs and also effect at mean delay which represents the strength of immune activation delay to boost the immunity of HTCs in the dynamics of the system.
The organization of this paper is as follows: Section 2 is devoted to the preliminary description of the delay differential equation governing the interactions between TCs, ECs, and HTCs.In Section 3, the qualitative and analytical study of the system is illustrated which includes positivity of the system, boundedness of the solutions, stability analysis of meaningful biological equilibrium point.We also investigated the existence of Hopf-bifurcation with respect to key parameters.In Section 4 the direction and stability of Hopf-bifurcation with respect to interaction delay is discussed.In Section 5 we carried out extensive numerical simulations to validate our analytical findings.Finally, Section 6 ends with concluding remarks about the key findings of the problem.

Mathematical model and preliminaries
It is known that instead of directly killing cancer cell, the helper-T cell help to activate immune-effector cell such as cytotoxic T-cells to attack cancer cell and also stimulate B-lymphocyte to secrete antibodies.As the overall process does not happen instantaneously, the time-lag is needed to make more realistic.Here, we consider the continuous delay kernel on helper-T cells while stimulating effector cells.The delay distribution equation with Gamma type kernel is more appropriate from the aspect of mathematical modeling and analytical tractability.For this purpose, we take Gamma distribution delay kernel as follows, G k (t) = ζ k+1 t k e −ζt k , k = 0, 1, 2... where, ζ > 0 indicates the strength of decay of helper-T cells of past states.If k = 0, then G(t) = ζe −ζt describes the weak kernel indicating the maximum weighted response of growth rate in current cell density, i.e., past cell density have decreasing effects.In the contrast, if k = 1,then G(t) = ζ 2 te −ζt describes strong kernel indicating the maximum influence on growth response past time.Here, we consider only the weak kernel in our study.Since delay kernel belongs to a family of probability distribution, we consider cell cycle time of helper-T cell positive and To apply the Gamma distribution delay kernel on helper-T cells, we define the integro-differential equation where, l > 0. le −lv is weight function for H(t − v) and monotonically increasing with respect to v. Whenever l is small, immune activation is considered previous time (t − v) depending on densities of HTCs, similarly when l is large, t takes present time.The present tumor-immune system interaction model is a modification of a model proposed by Dong et al. [18].So, the delay kernel function as immune activation continuous delay is incorporated in the proliferation of effector-cells (ECs) which is stimulated by helper-T cells (HTCs).The following system of DDEs are given where  For the sake of simplicity, we non-dimensionalize model system (2.1) by using the following transformation: We take T 0 = E 0 = H 0 = 10 6 cells/cm 3 [18,30] for the improvement of numerical performance.We also consider τ = τ and replace ν by t.Hence, we get the following dimensionless model equations with initial conditions φ = (φ 1 , φ 2 , φ 3 , φ 4 ), defined in the space where τ = max{τ 1 , τ 2 }φ i (ξ) ≥ 0, i = 1, 2, 3, 4, ξ ∈ [−τ, 0] and C is Banach space of continuous functions.φ : [−τ, 0] → R 4 + , with suitable sub-norm, and R 4 + = {(x, y, z, g) : x ≥ 0, y ≥ 0, z ≥, g ≥ 0}.
3. The qualitative analysis of the model where the mapping F : ), then the system can be written as )e ( t s e r 1 (τ 1 +ξ) dξ) ]ds = M 1 (say).

Stability analysis
In order to study the local dynamics of the system (2.2), the Jacobian matrix J E * of the system (2.2) about positive equilibrium points E * (x * , y * , z * , g * ) is given by and corresponding characteristic equation is, where λ is eigenvalue and Theorem 3.3.The tumor-free equilibrium point E 0 of the model system (2.2) is locally asymptotically stable if η 1 < η < η 0 for any τ 1 , τ 2 > 0.
Proof.The characteristic equation (3.2) at E 0 = (0, y 0 , z 0 , g 0 ) becomes Now we investigate the local stability of tumor-presence equilibrium with influence of discrete time delay τ 1 , τ 2 .Case I: The characteristic equation at E 1 (x 1 , y 1 , z 1 , g 1 ) is given by where, By well-known Routh-Hurwitz criterion, all the roots of equation (3.4) have negative real parts if and only if It is clear that the tumor-presence equilibrium is locally asymptotically stable if the condition (3.5) holds.
Case II: τ 1 > 0, τ 2 = 0.The characteristic equation (3.2) becomes The characteristic equation at E 1 (x 1 , y 1 , z 1 , g 1 ) is given by where, It can be mentioned that the characteristic equation (3.7) containing λ, τ has infinitely many eigenvalues.So, Routh-Hurwitz criterion is unable to investigate the locally asymptotical stability of equation (3.7).The local asymptotical stability of E 1 occurs if all roots of equation (3.7) have negative real parts, otherwise, stability is lost if a purely complex root crosses the imaginary axis from left to right.To investigate the sign of the real part of the roots, we use analytic arguments of Ruan and Wei [35].We substitute λ = iθ into the equation (3.7) to get periodic solutions which are significant in tumor dynamics.Separating the real and imaginary parts, we obtain the following transcendental equations Squaring and adding both the equations of (3.7), we have which implies where, Then simplified version of equation (3.9) is where k = 1, 2, 3, 4; j = 0, 1, 3....., and ±iθ k is pair of pure imaginary roots of equation (3.9) with τ j 1,k given by (3.11).We define Case III: The characteristic equation (3.2) becomes The characteristic equation at E 1 (x 1 , y 1 , z 1 , g 1 ) is given by where, Now, we put λ = iφ in equation (3.14) and proceed in similar way as (3.8)-(3.10).We obtain the expression for time delay τ 2 as where k = 1, 2, 3, 4; j = 0, 1, 3....., and ±iφ k is pair of pure imaginary roots of equation (3.13) with τ j 2,k given by (3.15).We define (3.16)

Analysis of Hopf Bifurcation
As pair of purely imaginary roots arise, we need to investigate the Hopf bifurcation analysis [26] in the system (2.2) and verify the transversality condition d(Reλ) dτ1 | τ1=τ * 1 > 0 which also preserves the conditions for the existence of periodic solution.We put ±iθ k in characteristic equation (3.7) such that |D(iθ k )| = |Q(iθ k )|,which also gives set of values of τ j k .Our main aim is to study the direction of motion of λ whenever τ 1 varies.Now we calculate Differentiating equation (3.7) with respect to τ 1 , we have Therefore, if F (θ * 2 ) = 0, the transversality condition holds and Hopf bifurcation occurs at θ = θ * , τ 1 = τ * 1 .Hence, delay based tumor model has a small amplitude periodic solution bifurcating from the tumor free equilibrium point as bifurcating parameter τ 1 crosses critical value τ 1 = τ * 1 , where τ * 1 is least positive value given by (3.11).Now we summarize the results as a theorem: Theorem 3.4.For τ 1 > 0, τ 2 = 0, assume that the condition (3.7) holds for the model system (2.2) with initial conditions (2.3).If e 0 < 0, i .e., a 2 0 − b 2 0 < 0.Then, i) the tumor presence equilibrium E 1 is locally asymptotically stable, if 2) experiences a Hopf-bifurcation around the tumor presence equilibrium E 1 at τ 1 = τ * 1 , where, (3.17) In similar way, we can formulate a theorem for time delay τ 2 .

Direction and stability analysis of Hopf bifurcation
From the previous section, we know that the model system (2.2) undergoes Hopf-bifurcation at E 1 under certain conditions.In this section we will study the stability and direction of Hopf-bifurcating periodic solution when the bifurcation parameter τ 1 passes through critical value τ 1 = τ * 1 using center manifold reduction of following the ideas of Hassard et al. [26].Now we translate the point E 1 = (x 1 , y 1 , z 1 , g 1 ) to the origin by using translator u 1 = x(t) − x, u 2 = y(t) − y, u 3 = z(t) − z, u 1 = x(t) − x, u 4 = g(t) − g.Let x i = u i (τ t) and τ = τ 0 + t, for simplicity τ 0 = τ j 1,k , where k = 1, 2; j = 0, 1, 2, 3 and µ = R. Then the system (2.2) can be expressed as a form of functional differential equation(FDE where (x Now we may take where δ is Dirac delta function.Let us define then the system(4.1)is equivalent to where Again we define a bilinear inner product R(0) and R * are known as adjoint operators.From the previous section, we have ±iτ 0 θ k which are eigenvalues of R(0) as well as of R * .Now we are interested to evaluate the eigenvectors R(0) and R * .Let p(ν) = (1, α 1 , α 2 , 0) T e iτ0θ k ν is eigenvector of R(0) corresponding to iτ 0 θ k , then R(0)p(ν) = iτ 0 θ k p(ν).By using (4.4),(4.5),(4.14) and the definition of R(0), it follows that, where, a 21 = −µ 1 y 1 e −λτ1 , a 31 = −µ 2 z 1 e −λτ2 , a 22 = −r 1 + ηg 1 + µ 1 yx1e −λτ1 and a 33 = −r 2 + µ 2 x 1 e −λτ2 to find α 1 and α 2 .Similarly, assume that p * (s) = G(1, ρ 1 , ρ 2 , ρ 3 )e iτ0θ k s is eigenvector of R * corresponding to eigenvalue−iτ 0 θ k to find ρ 1 , ρ 2 and ρ 3 .In order to verify p * (s), p(ν) = 1, it is essential to evaluate the expression of G.By using (4.15)G is calculated as follows: Therefore, we choose Now, we need to evaluate the coordinates describing the center manifold C 0 .Since µ = 0, assume that x t be solution of (4.13) at µ = 0.For this, we define, By center manifold C 0 , we get, W (t, ν) = W (u(t), u(t), ν), where u and u are the local coordinates for the center manifold C 0 in the direction of p and p * .W (u, u, ν) will be real if x t is real.Now our aim is to find only real solutions.Since µ = 0, we get which can be expressed as u(t) = iτ 0 θ k u + h(u, u), where Now, it can be obtained from (4.11) and (4.12) that From the definition of h(u, u), we get Comparing with the coefficients of equation (4.13), we have As the expression for 21 contains W 20 and W 11 , we need to classify W 20 and W 11 .Using (4.13) and (4.16), we obtain where On the other side, by the center manifold C neat at origin, it can be written as Ẇ = W u u + W u u Putting the corresponding series into (4.14) and after comparing, we have Again comparing the coefficient with (4.15), we get Using (4.18) and definition of R, it is obtained that Ẇ20 = 2iτ 0 θ k W 20 (ν) + h 20 p(ν) + g 02 p(ν).Therefore, Again in the same way, from (4.16) and (4.18), we get where 2 ) are constant vectors in R 4 .Now, we need to evaluate an suitable form of constant vectors V 2 and V 2 .From definition of R and (4.16), we get Again we use (4.14) to get N 20 and N 11 as follows: 1 can be obtained.Again, substituting (4.20) and (4.21b) into (4.22b),we have Hence, 2 can also be evaluated by Crammer's rule.Thus, we get h 21 using W 20 (ν) and W 11 (ν) from (4. 19) and (4.20) respectively.Therefore, we can evaluate the following values ) Here, κ 1 indicates the direction of Hopf-bifurcation, κ 2 indicates the stability of Hopf-bifurcating periodic solution and T 2 represents periodic solution at critical value τ = τ j k .By the center manifold theorem and results of Hassard et al., the properties of Hopf-bifurcation can be summarised as a theorem.

Numerical simulations and biological implications
In this section, we present the numerical results of the model system to study the effect of immune-activation distributed delay, discrete delays and different parameters.The parameter values given in Table 1 have been used to obtain numerical simulations.For parameter values b = 4, µ 1 = 0.04, µ 2 = 0.015 and η = 0.0375, the tumor free equilibrium and tumor presence equilibrium are obtained as E 0 (0, 1.02509, 6.90909, 6.90909) and E 1 (0.362843, 1.63481, 7.66788, 7.66788) respectively.After a few calculations, we have η 0 0.1743 and η 1 0.0437.From Theorem 3.3, it can be noted that the tumor-free equilibrium is asymptotically stable for 0.0437 < η < 0.1743.Using Theorem 4.1, the values of κ 1 , κ 2 and T 2 can be evaluated.From explicit expressions (4.23) and given parameter values, one can calculate κ 1 = 59.6107(> 0), κ 2 = −2.5379(<0) and T 2 = 2.1523(> 0).As κ 1 > 0 and κ 2 < 0, the Hopf bifurcation is supercritical and the bifurcating periodic solution is stable and its period increases (as T 2 > 0) when τ 1 crosses the critical value τ * 1 from left to right.

Influence of µ 1 : Stimulation rate of ECs in presence of TCs
Figure 2a-f captures the effect of the stimulation rate µ 1 of effector cells (ECs) in presence of tumor cells (TCs) on the dynamics of tumor-immune interaction.For µ 1 = 0.038 and τ 1 = 0.2 days < τ * 1 (bifurcation point) = 0.282 days = 6.768 h, the system shows an unstable equilibrium with a regular periodic behavior to tumor presence equilibrium point E 1 .This implies that expansion of tumor cells in immunizing process is under host-defensive as well as tumor-enhancement actions of immunity [7].As the value of µ 1 increases, the system exhibits a change from unstable to stable regime.Figure 2d-f shows stable dynamics of the system at µ 1 = 0.05 and τ 1 = 0.2 days = 4.8 h.The system admits Hopf bifurcation at µ * 1 (critical point) = 0.045 for τ 1 = 0.2 days (Fig. 7a-c) and sustains periodic oscillation.With increasing µ 1 , the amplitude of periodic oscillation for tumor cells changes to a damped oscillatory solution.Ultimately, the effector cells and helper T cells are capable of lysing tumor cells.At this stage tumor cells are non-invasive [28].

Influence of µ 2 : Stimulation rate of HTCs in presence of TCs
Figure 3a-f illustrates how the dynamics of tumor-immune interaction depends on the stimulation rate µ 2 of helper T cells (HTCs)in presence of tumor cells(TCs).For µ 2 = 0.215 and τ 2 = 0.25 days = 6 h, the system approaches a steady state with damping oscillatory behavior(Fig.3c).With increasing µ 2 , the system undergoes P. DAS ET AL.  (d,e) present the phase portrait of TCs-ECs, TCs-HTCs and (f) presents time evaluation of system for µ 2 = 0.28, τ 2 = 0.25 days = 6 h. a change in stability property and a small amplitude of periodic oscillation is demonstrated.Thus the system exhibits a Hopf bifurcation at µ * 2 (bifurcation point) = 0.223 (Fig. 8a-c).Finally, beyond Hopf threshold, a stable oscillation is observed with very low amplitude governing the persistent oscillatory behavior of tumor cells that may be mentioned as dormancy [23,30].

Influence of τ 1 : Interaction delay between ECs and TCs
The dynamics of tumor-immune interaction due to change in discrete time delay τ 1 is observed in Figure 4a-f.For τ 1 = 0.21 days < τ * 1 (bifurcation point) = 0.282 days and b = 4, the system shows stable equilibrium (Fig. 4c) and approaches the tumor presence equilibrium E 1 with a damping oscillatory behavior.As τ 1 increases, the system experiences a changing dynamics from stability to instability.At τ 1 = 0.282 days = 6.768 h, a Hopfbifurcation is observed (Fig. 9a) in the system.Consequently high amplitude of periodic oscillation is observed (Fig. 4f) as τ 1 increases.This means that the 'patient's situation' is in advance stage due to increase of interaction time delay τ 1 , which corresponds to metastatic state [31].

Influence of τ 2 : Maturation delay between HTCs and TCs
Figure 5a-d show the effect of maturation delay τ 2 on the dynamics of tumor-immune interaction.For τ 2 = 5 days and b = 4, a very small periodic oscillatory behavior at steady-state value E 1 is observed in Figure 5a which is dormant state of tumor.It can be observed in Figure 5b that the immune system can control tumor burden for a short time.As τ 2 increases, a quasi-periodic behavior has flourished as the system shows rare and irregular long periodic oscillations.This means a long delay in immune response may develop the tumor to more incursive state [7].Also further increase of τ 2 , there occurs a stability switch (Fig. 5c), the system shows unstable behavior leading to quasi-periodic phenomena.This indicates that tumor growth is beyond controlled.As a result, tumor cells are highly dependent on healthy tissue cells.behavior of the system of tumor presence equilibrium E 1 is observed.As the value of b increases, the oscillation of tumor cells decreases after some recovery time(Fig.6b).For τ 1 = 0.35 days = 8.4 h, the system exhibits Hopf-bifurcation at b = 0.218.With increasing b, the amplitude of oscillation tumor cells increase for tumor cells (Fig. 6d) and there is a change in dynamics from damped oscillations to periodic oscillation, i.e., effector cells and helper T cells do not have any impact on tumor cells.

Bifurcation analysis
The analysis of Hopf-bifurcation is more important in the dynamics of tumor-immune interaction due to the occurrence of limit cycle around the critical point and resulting stable and unstable periodic solutions.The periodic solutions provide the clinical symptoms of tumor dynamics as it indicates the oscillation of tumor growth levels around the equilibrium point due to lack of any kind of treatments.Now, our aim is to investigate the different situations for which we change the parameter values.Here µ 1 , µ 2 , τ 1 , τ 2 and b are bifurcating parameters and remaining parameter values are given in Table 1.

Bifurcation diagram for µ 1
Figure 8a-c shows the bifurcation of the system with respect to µ 1 , the stimulation rate of effector cells (ECs) in presence of tumor cells (TCs).Below the threshold value µ * 1 (bifurcation point) = 0.045, the tumor cells are in the unstable state.Beyond the Hopf threshold, the population of tumor cells, effector cells and helper T cells gradually reduce the periodicity.Finally, the system reaches a stable steady state.This implies that the stimulation rate does not affect the system in range µ 1 > µ  1.

Bifurcation diagram for µ 2
Figure 9a-c illustrates the Hopf bifurcation of the system with respect to the proliferation rate µ 2 of helper T cells (HTCs) in presence of tumor cells (TCs) and effector cell (ECs).In the range 0 < µ 2 < µ * 2 (bifurcation point) = 0.223, the tumor cells are controlled by effector cells and helper T cells and a stable limit cycles exists.Beyond the critical value µ * 2 = 0.223, the growth of tumor cells highly depends on µ 2 .

Bifurcation diagram for τ 1
Figure 10a shows the Hopf-bifurcation of the system with respect to interaction delay τ 1 .In the range 0 < τ 1 < τ * 1 (bifurcation point) = 0.282 days = 6.768 h, the system shows the stable steady state of tumor presence equilibrium point E 1 and this implies that the density of tumor cells is well controlled by ECs and HTCs.At this stage tumor characteristic is non-metastatic and dormant.Beyond the critical value τ * 1 , the tumor cells exhibit stable oscillator behavior with higher amplitude.This implies that the growth of tumor cells depend on τ 1 .1.

Bifurcation diagram for b
Figure 11a demonstrates Hopf-bifurcation with respect to b, the strength of distributed delay in presence of interaction delay τ 1 .The stability of the system changes around the critical value b * τ1 (bifurcation point in presence of τ 1 ) = 0.218.Below the threshold value, the system shows damped oscillations with stable limit cycles of tumor presence equilibrium E 1 .Beyond the critical value b * τ1 , tumor cells show a regular periodic oscillation.Similarly in Figure 11a, the system admits Hopf-bifurcation with respect to b around b * τ2 (bifurcation point in presence of τ 2 ) = 0.23 in presence of maturation delay τ 2 in Figure 11b.It is observed that the immune activation is highly dependent on the population of helper-T cells (HTCs).At time t, the larger value of b shows weak delay effect while the smaller value of b at past time t − v may demonstrate strong delay effect.Figure 11c illustrates Hopf-bifurcation diagram.The system shows stable steady state of tumor presence equilibrium with increase of b, strength of immune activation delay in presence of discrete time delays.At b * τ1τ2 (bifurcation point in presence of τ 1 and τ 2 ) = 0.361, Hopf-bifurcation has occurred.Furthermore, it can be noticed that stability regime for parameter b is relatively larger than that of any discrete delay.This implies that uniformly activated HTCs help in ECs stimulation and also controls tumor burden simultaneously.

Concluding remarks
Nowadays it is more essential to consider natural phenomena by introducing time delays in the demonstration of tumor-immune interaction.A simple mathematical model can play an important role to interpret the nature of a tumor-immune dynamics and help to reveal better treatment technique.Indeed, it may be challenging to deal with experimental studies to trace real delay patterns in T-IS interplay where our present work can give some important insights.
In the previous study [17], the existing tumor-immune interaction model considered the discrete delays between tumor cell-effector cell and tumor cell-helper T cell.Then the dynamics have also been investigated in the presence of delays.However, the activation delays can even stabilize the immune-control equilibrium.But current study reveals that the activation of a helper-T cell is not an instantaneous process.Helper-T cells become activated by interaction with an antigen-presenting cell such as macrophages.After this, the activated HTCs stimulate the immune-effector cell.Due to these reasons, in this research, we have incorporated a continuous delay kernel in HTCs to illustrate the inner scenarios of tumor-immune interaction.Our analytical findings are well complemented with numerical simulations.Here, we have demonstrated positivity of the solutions, boundedness of the system and investigated local stability by Routh-Hurwitz criterion.The qualitative analysis of Hopf-bifurcation with respect to different parameters has been studied.The direction of Hopf-bifurcation and the stability of bifurcating periodic solution based on center manifold theorem has been derived explicitly.While observing the dynamics, the influence of the immune activation delay along with discrete delays has been highlighted.For understanding complex biological scenarios, we have analyzed numerical simulations to validate the analytical findings.It is observed that beyond the Hopf-threshold, our model generates periodic oscillation which is biologically appropriate and realistic as similar kind of dynamics are observed in cancer patients.Moreover, few simulations exhibit that increase of delays may cause a change from stable to unstable dynamics.Otherwise, local stability allows the tumor to reach near values for their carrying capacity(maximum tumor size).As stability switch has occurred in the system, this phenomena interprets that immune system receives the message of oscillatory tumor growth.As a result, the immune system becomes more effective to annihilate tumor cells.But, for excessive delay in immune interaction, the immune system becomes less effective to control tumor growth.
In the beginning, we have observed that effector cells are able to control tumor growth for stimulation rate µ 1 > 0.045.Similarly, if stimulation rate µ 2 of helper-T cells crosses critical value µ * 2 = 0.045, the tumor becomes invasive.Also, interaction delay τ 1 plays a key role, the tumor is beyond the control of effector cells if τ 1 > τ * 1 ( = 0.282 days = 6.768 h) i.e.. it takes more time to stimulate effector cells to initiate annihilation of tumor cells.Hence, the parameters µ 1 and µ 2 play a key role to boost HTCs and ECs in controlling tumor burden.
We have investigated the influence of immune activation delay strength b.It has been discussed that the immune response is a continuous process.So if we increase the strength of distributed delay, the immune response can be increased equivalently.We have observed that tumor is stable in presence of τ 1 for below critical value b * τ1 = 0.218 days = 5.232 h and similarly in presence of τ 2 , tumor is stable for below b * τ2 = 0.230 days = 5.52 h.But for below of b * τ1τ2 = 0.361 days = 8.664 hours, the system remains stable for larger time.Nevertheless, distributed delay along with discrete delays play a crucial role for the larger stability regime of tumor growth than any individual delay or both delays (Fig. 11a-c).Again we have taken τ 1 > 0 and τ 2 = 0 for tumor angiogenesis.There exists Hopf-threshold value for τ 1 crossing which unstable steady state is manifested due to longer delay.This implies a larger interaction delay couldn't affect the tumor burden.Here "Sneaking through" phenomena are observed [30], which refers low doses of tumor cells fail to mount an effective tumorimmune response and gradually tumor grows, similarly medium doses of tumor cells lead to tumor eradication and large doses develop tumor cells through immune defenses.
However, it has been observed that the stability regime for b, the strength of continuous delay is relatively larger than that of the presence of any delay.From a qualitative point of view, a set of parameter values have been derived which may design new diagnostic and control technique to interrupt tumor growth.This assures that the uniformly activated helper-T cells help in the stimulation of ECs and also control the tumor burden simultaneously which indicates the applicability of our proposed research.We would like to deepen our analysis of the study of CD8 + T cells and tumor cells in the future.
and H(t) represent the densities of tumor-cells (TCs), effector-cells (ECs) and helper-T cells (HTCs) respectively.The first equation represents the rate of change of TCs densities, the first term of which describes the logistic type growth of TCs, a is the intrinsic growth rate and m −1 is maximal cell burden or carrying capacity.The second term corresponds to the inhibition of immune ECs due to the presence of TCs at a rate of n.The second equation corresponds to proliferation enhancement rate of ECs, where the first term represents a constant flow at a rate of s 1 of mature ECs into the region of TCs localization.The second term indicates competition between TCs and ECs, where k 1 is the stimulation rate of ECs lysed and TCs debris.The third term represents ECs natural death with an average half-life 1/d 1 .The fourth term describes the stimulation of ECs with a rate of p by kernel-based HTCs K(t).The third equation represents the rate of change of HTCs, where the first term represents the birth of HTCs, the second term describes stimulation rate of HTCs at k 2 in the localization of tumor-specific antigens and last term represents the HTCs natural death with average half-life 1/d 2 .The fourth equation represents the rate of K(t); where the first term describes the increase of HTCs.The second term represents a loss of K(t) with a strength of l.It is assumed that all parameter value are positive.Figure1illustrates the interactions between Tumor-cells(TCs), Effector-cells(ECs) and Helper T-cells(HTCs) with their cellular environment.

5. 5 .
Figure 6a-d illustrate how the strength of distributed delay, b can change the dynamics in presence of interaction delay τ 1 .For b = 0.18 < b * (bifurcation point) = 0.218 and µ 1 = 0.35, a stable damping oscillatory

Figure 7 .
Figure 7. (a) and (c) show the three dimensional phase portrait diagram for (a) b = 0.18 and (c) b = 0.35 with τ 2 = 4 days.(b) and (d) present the time evaluation curve of system for (b) b = 0.18 and (d) b = 0.35 with τ 2 = 4 days.

Table 1 .
Parameter values for numerical simulations.