Dual role of delay effects in a tumour–immune system

ABSTRACT In this paper, a previous tumour–immune interaction model is simplified by neglecting a relatively weak direct immune activation by the tumour cells, which can still keep the essential dynamics properties of the original model. As the immune activation process is not instantaneous, we now incorporate one delay for the activation of the effector cells (ECs) by helper T cells (HTCs) into the model. Furthermore, we investigate the stability and instability regions of the tumour-presence equilibrium state of the delay-induced system with respect to two parameters, the activation rate of ECs by HTCs and the HTCs stimulation rate by the presence of identified tumour antigens. We show the dual role of this delay that can induce stability switches exhibiting destabilization as well as stabilization of the tumour-presence equilibrium. Besides, our results reveal that an appropriate immune activation time delay plays a significant role in control of tumour growth.


Introduction
The word 'tumour' in reality expresses an entire and wide family of high-mortality pathologies of cellular proliferation which may cause diseases eventually [9]. The immune system is a complicated network of cells and signals that has evolved to protect the host from pathogens and body from tumour cells (TCs). Any new substance in the body that the immune system does not recognize can raise an alarm and cause the immune system to attack it. Sometimes the immune system can recognize TCs, but the response might not be strong and prompt enough to destroy TCs [5]. Effector cells (ECs) are the activated immune system cells, such as cytotoxic T lymphocytes (CTLs) [15]. TCs are mainly the prey of the ECs, on the other hand, proliferation of ECs is stimulated by the existence of TCs. However, TCs also cause a loss of ECs, and intensity of an influx of ECs may rely on the size of the tumour [8]. Helper T cells (HTCs) are stimulated through antigen presentation by macrophages or dendritic cells and release the cytokines to stimulate proliferation of ECs to kill more TCs. The tumour-immune competitive system is complicated, being nonlinear and evolutive in some extent [10].
To capture the dynamics of tumour-immune system interaction, a lot of mathematical efforts have been made [1-4, 6-11, 13-20]. For example, a classical model describing the CTL response to the growth of an immunogenic tumour was presented in [16]. Kirschner et al. first introduced immunotherapy into their models which can explain both short-term tumour oscillations in tumour sizes and long-term tumour relapse [15]. d'Onofrio et al. [8] introduced a general framework for modelling tumour-immune system competition and immunotherapy using the qualitative theory of differential equations. de Pillis et al. [3] proposed a mathematical model of tumour-immune interaction with chemotherapy and strategies for optimally administering treatment. Recently, Dong et al. [7] studied the role of HTCs in the tumour-immune system: where T(t), E(t) and H(t) are the densities of TCs, ECs and HTCs , respectively. The logistic term aT(1 − bT) describes the growth of TCs, where a denotes the maximal growth rate of the TCs population and b −1 denotes the maximal carrying capacity of the biological environment for TCs. n represents the loss rate of TCs by ECs interaction. s 1 is referred to as the constant flow rate of mature ECs into the region of TCs localization. ECs have a natural lifespan of an average 1/d 1 days. k 1 is ECs stimulation rate by ECs-lysed TCs debris. p is an activation rate of ECs by HTCs. s 2 is referred to as birth rate of HTCs produced in the bone marrow. HTCs have a natural lifespan of an average 1/d 2 days. k 2 is the HTCs stimulation rate by the presence of identified tumour antigens. To consider a weak proliferation of ECs stimulating by TCs, we can simplify this system by ignoring the ECs stimulation from TCs (that is to omit the term +k 1 T(t)E(t) in system (1)), while keeping the essential dynamics properties of the original system. Then, we obtain the simplified system: with initial conditions As suggested in [16], we choose the scaling E 0 = T 0 = H 0 to improve the performance of numerical simulations. Then, system (2) can be rescaled by introducing the new variables: x = T/T 0 , y = E/T 0 , z = H/T 0 , τ = nT 0 t, α = a/nT 0 , β = bT 0 , σ 1 = s 1 /nT 0 2 , δ 1 = d 1 /nT 0 , ρ = p/n, σ 2 = s 2 /nT 0 2 , ω = k 2 /n, δ 2 = d 2 /nT 0 . By replacing τ by t, system (2) becomes with initial conditions x(0) > 0, y(0) ≥ 0 and z(0) ≥ 0, where x(t), y(t) and z(t) denote the dimensionless density of the TCs, ECs and HTCs population, respectively. In reality, the immune system will take more time to give a suitable response after recognizing the non-self cells. Thus, we need to take into account the effect of time delay to describe the proliferation of ECs stimulating by HTCs. There are mainly two ways to incorporate delay effects into the ordinary differential equations (ODEs) system. One way is to directly introduce a delay term (e.g. +ρy(t)z(t − θ)) with θ > 0. Following this, Dong et al. [6] constructed the dynamical behaviour of a tumour-immune system interaction model with two discrete delays. And the other way is to introduce a kernel function as an additional compartment. In this study, we use a function G(t) = t −∞ be −b(t−s) z(s) ds to represent the immune activation delay taking the following form: where b is a positive constant expressing the strength of delay. In fact, G(t) can be rewritten as G(t) = ∞ 0 be −bu z(t − u) du. Note that be −bu is the weight function for z(t − u) and monotonically decreasing with respect to u. The mean delay of the weight function be −bu is 1/b. Hence, the immune activation depends on the density of HTCs at the present time t when b is large. On the other hand, when b is small, the activation depends also on the density of HTCs at the past time t−u. This study is aimed to investigate the role of delay effects in the tumour-immune system.
The rest of the paper is organized as follows. The stability of tumour-free and tumourpresence equilibria of systems (3) and (4) was discussed, respectively, in Section 2. The numerical simulations were performed in Section 3. The conclusion and discussion were given in Section 4.
The Jacobian matrix of Equation (3) at E * is The characteristic equation of J(E * ) takes the form where According to the Routh-Hurwitz criterion, we know that all roots of Equation (6) have negative real parts if and only if b Thus, the sufficient condition for stability of E * is given by where y * = α(1 − βx * ) and z * = σ 2 /(δ 2 − ωx * ). Thus, we have the following result.

Dynamics of system (4)
Following the similar analysis of the steady states in system (3), for system (4) we obtain the tumour-free equilibrium, And we can also obtain the unique tumour-presence equilibrium E + = (x + , y + , z + , . Now, we discuss the stability of two equilibria E − and E + of system (4).
We compute the Jacobian matrix of system (4) at E − , denoted by It can be seen that all the eigenvalues are negative if and only if ρ 1 < ρ < ρ 0 . Thus, we have the following result.
The Jacobian matrix of system (4) at E + is The characteristic equation of J(E + ) takes the form where According to the Routh-Hurwitz criterion, we know that all roots of Equation (8) have negative real parts if and only if B 1 > 0, (6). Thus, the sufficient condition for stability of E + is given by where Thus, we have the following result.

Theorem 2.4: System (4) has one tumour-presence equilibrium E
When ω tends to 0, all coefficients of Q(b) are positive. Notice that b 4 → 0 as ω → 0. We have Q(b) > 0 for all b > 0 and we have no stability change.
First, let us suppose that b 1 b 2 − b 3 > 0, that is, a positive equilibrium point E * of Equation (3) is locally asymptotically stable. We are interested in the possibility of stability change by the delay effect b. Note that b 2 1 > 2b 2 > 0 by the definition (6). Further P 2 > 0 when P 1 > 0. Hence, under the assumption b 1 b 2 − b 3 > 0, we have the following three cases: Hence, for both cases, we have that Note that for Case (III) (P 0 > 0, P 1 < 0, P 2 < 0, P 3 > 0), Equation (9) is always satisfied if 9P 0 P 3 < P 1 P 2 . Thus, we have the following result. The equilibrium E + is locally asymptotically stable for 0 < b < b or b > b and unstable for b < b < b. If Equation (9) is not satisfied, the equilibrium E + is always locally asymptotically stable.
Second, let us consider the opposite condition b 1 b 2 − b 3 < 0. Under this condition, a positive equilibrium E * of Equation (3) is unstable and we have a periodic solution. It is easy to check that we have P 0 > 0 and P 1 < 0, P 2 < 0, P 3 < 0. Note that This implies that small b > 0 (i.e. strong delay effect) can stabilize the equilibrium point. For small b, regardless of all Cases (I)-(III), we have Q(b) > 0 since P 0 > 0. Therefore, we have the following result.

Proposition 2.6:
The equilibrium E + is always locally asymptotically stable if b is very small.
In Figure 1, two bifurcation curves are drawn for system (3) (b = +∞) and system (4) (b = 0.004) in the same ρ − ω parameter plane. Both E * and E + are stable below the corresponding bifurcation curve and are unstable above the bifurcation curve. The bifurcation curves in Figure 1 are derived from equations b 1 (ρ, ω) ω) = 0, respectively. The shape of the bifurcation curve for system (3) is similar to the one in [7], which means that the simplified system (3) keeps the dynamics properties of the original model (1). When we consider a small b = 0.004 for system (4), we found that the shape of the bifurcation curve has been changed. Correspondingly, the stable and unstable regions are also changed, which offers the possibility of stability switch due to the delay effect.

Conclusion and discussion
In this paper, first, we reduced the previous system in [7] by ignoring a weak direct immune activation by the TCs. We found that the simplified three-dimensional system (3) can keep the intrinsic dynamics properties of the original model (1) due to the similar expression (7) and the shape of the bifurcation curve (see Figure 1). Second, we incorporated the time delay into the ODE model to modelling the lag of the activation of ECs from HTCs. If we directly incorporate the delay term, it is difficult to calculate the transcendent characteristic equation. Alternatively, we used a function G(t) = ∞ 0 be −bu z(t − u) du to model this immune activation delay effect at the cost of increasing one dimension. The immune activation depends on the density of HTCs at the present time t when b is large representing a weak delay effect, on the other hand, the activation depends on the density of HTCs at the past time t−u when b is small indicating a strong delay effect.
Compared with system (3), the introduction of new parameter b > 0 in system (4) did not change the values of tumour-free and tumour-presence equilibria; however, it can substantially change the dynamics properties. Specifically, the delay can lead to stability changes exhibiting destabilization as well as stabilization of the tumourpresence equilibrium, which plays a dual role. We provided some numerical simulations to show that dynamics can be changed at the branch of the tumour-presence equilibria E * and E + , resulting in Hopf bifurcation. In Figure 1, both E * and E + are stable below the corresponding bifurcation curve and are unstable above the bifurcation curve. To show the possible stability change of E * and E + , we chose four points (0.01, 0.00035), (0.03, 0.00035), (0.0425, 0.00035) and (0.043, 0.00035) in the ρ − ω plane. By comparing the left panel for system (3) and the right panel for system (4) in Figure 2, we found that the introduction of time delay indeed caused the dynamics change.
The bifurcation curves (b = 0) are satisfied with Q(ρ, ω) = 0 for each fixed b. When we consider a very small b (i.e., a strong delay effect, that is, the activation of ECs by HTCs is very slow), the stable region can be very large. And the stable region will tend to infinite as b tends to 0 (see the grey curve b = 0.0001, cyan curve b = 0.001, blue curve b = 0.004 and magenta curve b = 0.01 in Figure 4). As a consequence, the tumour-presence equilibrium is always stable and dominant under the condition b 1 b 2 − b 3 > 0 ensured by Proposition 2.6, which means that TCs cannot be eradicated from the body. When we further increase b, the bifurcation curve will go down and the stability region in the ρ − ω plane will become smaller. But the bifurcation curve can never reach the ρ axis, because when ω is small, the stability region always exists regardless of any changes of ρ and b. There must exist a critical valueb, which makes the stability region to be minimal. When we increase b to exceed thiŝ b, the bifurcation curve of system (4) will go up and stability regions will become larger (see the green dashed curve b = 0.1, yellow dashed curve b = 1 and red dashed curve b = 10 in Figure 4). Finally, the bifurcation curve of system (4) will tend to the bifurcation curve of system (3) as b tends to infinity (see the red dashed curve b = 10 and the black curve b = +∞ of system (3)). In other words, the decrease of delay will decrease first and then increase the stability region. In a word, a critical time delay exists and makes the stable region of tumour-presence equilibrium consisting of Q(b) = 0 and axes to be minimal. This means that an appropriate immune activation time delayb can minimally control the tumour growth.
In system (4), we have introduced a delay effect by a weak kernel, i.e., F(u) = be −bu (b > 0). For the further study, we can introduce a delay effect by a non-monotone kernel, i.e., F(u) = b 2 ue −bu (b > 0).

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This research was supported by Grants-in-Aid Scientific Research (C) No. 26400211, Japan Society for Promotion of Science.