{'fraction_non_alphanumeric': 0.057842889618590555, 'fraction_numerical': 0.025157868148522354, 'mean_word_length': 3.859106529209622, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
COMPLICATED DYNAMICS OF TUMOR-IMMUNE SYSTEM INTERACTION MODEL WITH DISTRIBUTED TIME DELAY

. In this paper, we propose a distributed delay model to investigate the dynamics of the interactions between tumor and immune system. And we choose a special form of delay kernel which combines two delay kernels: a monotonic delay kernel representing a fading memory and a nonmonotonic delay kernel describing a peaking memory. Then, we discuss the eﬀect of such delay kernel on system dynamics. The results show that the introduction of nonmonotonic delay kernel does not change the stability of tumor-free equilibrium, but it can induce stability switches of tumor-presence equilibrium and cause a rich pattern of dynamical behaviors including stabilization. Moreover, our numerical simulation results reveal that the nonmonotonic delay kernel has more complicated eﬀects on the stability compared with the monotonic delay kernel.

1. Introduction. A tumor (neoplasm) is an abnormal mass of tissue which can be benign, pre-malignant, or malignant. The immune response begins when tumor cells (TCs) are identified by many different kinds of immune cells: for example, B cells producing and secreting antibodies into the blood, cytotoxic T lymphocytes (CTLs) considered as effector cells (ECs) with the ability to kill TCs, Helper T cells (HTCs) secreting interleukins and hence stimulating both B and T cells, and so on [27]. Although the functions of the immune system during tumor development are very complicated, being nonlinear and evolutive in some extent [10], the mechanisms underlying the interactions between the tumor and immune system are critically important for understanding tumorigenesis and immunotherapy [15,34].
Enormous modeling efforts have been made to explore the dynamics of tumorimmune system interaction since the early 1990s and have evolved to capture much more complex aspects of the immune response [1,3,4,5,6,7,25,9,10,11,12,13,14,15,8,17,18,20,22,24,26,27,28,29,30,33,34,36]. For example, Kuznetsov et al. [22] presented a well-known model describing the CTL response to the growth of an immunogenic tumor, where the sneaking through mechanism was described: the tumor stays in a very small size for a relatively long time period, and subsequently grows dangerously large. Kirschner et al. [20] first introduced immnunotherapy into their models which can explain both short-term tumor oscillations in tumor sizes as well as long-term tumor relapse. Galach [18] introduced the immune reaction time delay into the simplification of Kuznetsov-Taylor model [22] and observed a state of the "returning" tumor. d'Onofrio et al. [9] proposed a general scheme to model tumor-immune system competition and immunotherapy. To study the influence of Natural Killer (NK) cells and CTLs in tumor surveillance, de Pillis et al. [7] presented a mathematical model validating by both experimental and human data. Mallet et al. [25] considered a hybrid cellular automata model to investigate the interactions of a tumor and the host immune system that involved important signalling molecules. Dong et al. [12] developed a three dimensional ordinary differential equations (ODEs) model to describe the interactions among TCs, ECs and HTCs, and observed oscillated behaviours. Furthermore, Dong et al. [13] incorporated two immune activation delays into their previous model [12] . It is quite common that oscillations occur in the immune system [31], however, there is debate whether the immune response delay can cause the oscillations [5,6,11,34]. Piotrowska et al. [28] took into account distributed time delays in a generalisation of tumor-immune system interaction model, and numerically solved the system with Erlang probability densities by applying the classic linear chain trick [23,21]. Thurley et al. investigated simple cell-to-cell communication circuit motifs by using response-time distributions [32]. For good reviews in the modeling of tumor-immune system interaction, we refer to [1,14] and references therein.
Recently, Yu et al. [36] formulated a three-dimensional tumor-immune system to study the interactions among TCs, ECs and HTCs as follows: Following the method in [22], system (1) can be rescaled by introducing the new variables and parameters: with initial conditions x(0) > 0, y(0) ≥ 0, z(0) ≥ 0, where x(t), y(t) and z(t) denote the dimensionless density of TCs, ECs and HTCs population, respectively.
In reality, the transmission of signals and biological processes does not occur instantaneously, but rather depends on the influence of the past history of the system. Therefore, it is often necessary to incorporate time delays into model foundation to reflect the real system. As we know, the activation of immune response against the development of tumors is very complicated and often delayed, such as, the immune activation delays of ECs by TCs and HTCs by TCs [13], and the immune activation delay of ECs by HTCs [36] and so on. Since the immune activation delay of ECs by HTCs is relatively longer than the others, we only consider such delay in this study.
The time lags can be generally classified into two types: discrete delay and distributed delay [19]. Here we use the distributed immune activation delay t −∞ F (t − s)z(s)ds to replace z(t) in the second equation of (2), then we obtain the following system: where F (t) is the distributed delay kernel. And we choose a special case of F (t) as the convex combination of a monotonic delay kernel and a nonmonotonic delay kernel, that is F (t) = c 1 F 1 (t) + c 2 F 2 (t), where F 1 (t) = be −bt , F 2 (t) = b 2 te −bt , and b, c 1 , c 2 ≥ 0, c 1 + c 2 = 1. F 1 (t) is a monotonic decreasing function describing a fading memory (see Figure  1(a)). It has a maximum value at the present time t = 0, that is to say, the present time has the greatest effect on the system. F 2 (t) is a nonmonotonic function describing a peaking memory (see Figure 1(b)). It has a maximum value at the past time t = 1 b , that is to say, the past time 1 b has a greater effect on the system. From Figure 1(b), we can see that the peak can shift to the left as the increase of b. F (t) is the linear combination of F 1 (t) and F 2 (t), and its properties are presented in Appendix A. Let By using the linear chain trick [23,21], system (3) can be transformed into the following system: where b = 1/τ is a positive constant, and τ represents average time delay. A relatively small b indicates a strong delay effect, which means that the system is more dependent on the past time, while a relatively large b indicates a weak delay effect, which means that the system is more dependent on the present time (see Figure 1).
When c 2 = 0, system (4) can be reduced to a four-dimensional system (with only a fading memory) which was well studied in [36] taking the following form: The analytical results on the stability of system (5) are summarized in Appendix B.
As far as we know, this is the first attempt to consider such nonmonotonic delay kernel in the modeling of tumor-immune system interaction. And we are aimed to investigate the impact of this nonmonotonic delay kernel on the dynamical behaviors of system (4). The rest of the paper is organized as follows. In Section 2, we discuss the local stability of tumor-free and tumor-presence equilibria of system (4). In Section 3, we perform the numerical simulations to draw the stability boundary curve for the tumor-presence equilibrium. In Section 4, we give a brief summary and a discussion of our findings.
2. Stability analysis. In order to investigate the steady states of system (4), we set dx/dt = 0, dy/dt = 0, dz/dt = 0, dG 1 /dt = 0 and dG 2 /dt = 0. Substituting x = 0 leads to the tumor-free equilibrium, E − exists, when ρ < δ1δ2 σ2 ≡ ρ 0 . Following the similar analysis in [36], we can also obtain the unique tumorpresence equilibrium E + = (x + , y + , z + , . Now we discuss the stability of two equilibria E − and E + of system (4). First, 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.
Second, we compute the Jacobian matrix of system (4) at E + , denoted by The characteristic equation of J(E + ) takes the form where and Applying the Routh-Hurwitz criterion, we know that all roots of (6) have negative real parts if and only if Following we choose the parameters which locate in the stability regions of E * = (x + , y + , z + ) of system (2) and E * = (x + , y + , z + , G 1 + ) of system (5). That is, we have . Therefore, the sufficient condition for the stability of E + is given by where It is easy to see that equation H(D 4 ) = 0 has two positive roots Note that D 4 is a decreasing function with respect to c 2 (see Figure 2(a)), and H(D 4 ) is a quadratic function with respect to D 4 (see Figure 2(b)). We define d = bB 3 + B 4 and d = bB 3 Figure 2 Figure 2(b)). The remaining matter is to determine the sign of H(d). In summary, we have two cases: There must exist a root c * 2 which satisfies D 4 (c * 2 ) = D * 4 . Thus we have the following result.
Since the condition c 1 + c 2 = 1 is satisfied, we choose c 2 as a bifurcation parameter. The stability boundary curves in Figure 3 are derived from H(c 2 , b) = 0 for each fixed ρ and ω. Note that function H given by (9) is now a function of c 2 and b, since E + does not depend on c 2 and b. In Figure 3(a), two stability boundary curves are drawn for system (4) in the c 2 − b parameter plane when ρ = 0.01 and ω = 0.00035 are fixed. E + is stable in the area outside these two stability boundary curves, while E + is unstable in the area between these two stability boundary curves. Particularly, when we fix c 2 = 0 (that is, c 1 = 1: case for a fading memory of system (5)), we can see that E + is stable, and becomes unstable, and then becomes stable again as the increase of b (decrease of delay effect), which is in accordance with the previous observations (dual role of delay effect) in [36]. Actually, system (4) can be reduced to system (5) when c 2 = 0. In this case, we observe that the instability region increases as the increase of c 2 . In Figure 3(b), one stability boundary curve is drawn for system (4) in the c 2 − b parameter plane when ρ = 0.03 and ω = 0.00035 are fixed. E + is stable in the area below the stability boundary curve, while E + is unstable in the area above the stability boundary curve. Similarly, when we fix c 2 = 0, we can see that E + is stable, and becomes unstable with increasing b, which is in accordance with the previous observations in [36]. Note that we observe the dual role of delay effect when an activation rate ρ of ECs by HTCs is relatively small even under a peaking memory (see Figure 3(a)).  (a) E + is stable in the area outside these two stability boundary curves, while E + is unstable in the area between these two stability boundary curves. (b) E + is stable in the area below the stability boundary curve, while E + is unstable in the area above the stability boundary curve.  In order to check how the parameter c 2 affects the stability of E + , we choose two representative values c 2 = 0.2 (the distributed delay kernel dominated by a fading memory) and c 2 = 0.8 (the distributed delay kernel dominated by a peaking memory) to draw stability boundary curves of E + in ρ − ω parameter plane for different b (see Figure 4). E + is stable in the area below the stability boundary curve, while E + is unstable in the area above the stability boundary curve. These curves satisfy H(ρ, ω) = 0 for each fixed c 2 and b. For the both cases c 2 = 0.2 and 0.8, we find that the stability region decreases first, and then increases with increasing b, which is similar to the system (5) with only a fading memory. This means that the introduction of c 2 does not change the variation tendency of stability, but it causes the change of stability region inducing the possibility of stability switch (see Figure 4(a) and (b)). In general, the increasing delay gives the tendency to promote the instability of the system. But for both cases in Figure 4, when b is relatively small (b < 0.1, indicating that delay effect is stronger), the stability region decreases with the increase of b, which is counterintuitive (against the general tendency). On the other hand, when b is relatively large (b > 0.1, indicating that delay effect is weaker), the stability region increases with the increase of b, which agrees to intuitive sense. Now let us compare two cases in Figure 4. When we fix b for the same value in Figure 4(a) and (b), we find the stability region for relatively small c 2 (the case dominated by a fading memory) is smaller than one for relatively large c 2 (the case dominated by a peaking memory). This implies that a peaking memory will make the system more stable. We now observe this is true or not by Figure 5.
We draw the stability boundary curves in ρ − ω parameter plane as c 2 varies for fixed b = 0.001 and b = 0.1. The stability boundary curves satisfy H(ρ, ω) = 0 for each fixed c 2 and b. From Figure 5(b), when b is relatively large (weak delay effect: b = 0.1), the stability region decreases monotonically as c 2 increases. This agrees to an intuitive sense of delay effects. That is, a peaking memory promotes the instability of the system. On the other hand, from Figure 5(a), when b is relatively small (strong delay effect: b = 0.001), the stability region becomes larger when c 2 increases as c 2 = 0(cyan curve), 0.2(magenta curve), 0.5(green curve), 0.8(yellow curve) for 0 < ρ < 0.04 and the stability region turns to be much smaller for c 2 = 1(red curve). The former phenomenon is also true for any ρ < ρ 1 , 0 < c 2 < 0.5. This implies that a peaking memory promotes the stability of the system when a fading memory is dominated (c 2 < 0.5) under any activation rate ρ of ECs by HTCs. After the peaking memory becomes dominant (c 2 > 0.5), still it promotes the stability under small ρ(< 0.04), but hinders the stability for large ρ. It is shown that the delay distribution function F (t) has a maximum peak F max only for c 2 > 0.5. Hence, dominated peaking memory brings instability under large activation rate ρ of ECs by HTCs. It is interesting to observe in Figure 5(a) that the stability region becomes smallest when c 2 reaches at its maximum (c 2 = 1). This suggests that there may exist an optimal valueĉ 2 which gives the maximum stability region at c 2 =ĉ 2 satisfying 0.5 <ĉ 2 < 1. This also shows that a peaking memory can promote the stability of the system. 4. Conclusion and discussion. In this study, we employed a special form of distributed delay with the kernel, F (t) = c 1 F 1 (t) + c 2 F 2 (t), (c 1 , c 2 ≥ 0, c 1 + c 2 = 1), consisting of a monotonic delay kernel F 1 (t) = be −bt representing a fading memory and a nonmonotonic delay kernel F 2 (t) = b 2 te −bt illustrating a peaking memory to model the time lag of the immune activation of ECs by HTCs. A relatively small b indicates a strong delay effect, while a relatively large b indicates a weak delay effect. c 2 = 0.5 is a threshold determining which delay kernel is dominated. When c 2 is relatively small (< 0.5), F (t) is dominated by F 1 (t). In this case, F (t) has the maximum value at t = 0 denoted by F max = F (0) = (1 − c 2 )b. The immune activation of ECs mainly depends on the density of HTCs at the present time t = 0. On the other hand, when c 2 is relatively large (> 0.5), F (t) is dominated by F 2 (t). In this case, F (t) has the maximum value at t = t * denoted by F max = F (t * ) = bc 2 e −2+ 1 c 2 . The immune activation of ECs mainly depends on the density of HTCs at the past time t = t * = 1 b (2 − 1 c2 ). We found that c 2 and b play an opposite role in the effect of F max .
The introduction of new parameter c 2 > 0 in system (4) does not change the values of tumor-free equilibrium E − and tumor-presence equilibrium E + . The analytical results showed that varying c 2 does not change the stability of E − . However, the analytical results as well as numerical simulation results showed that the stability of E + can be changed. We presented some numerical examples to show the dynamical changes at the branch of E + . The dual role of delay effect is not changed by introduction of a peaking memory (see Figure 3). We showed the effect of c 2 on the stability of E + (see Figures 4 and 5) and suggest the existence of an optimal c 2 to make the stability region maximum. From Figure 7, we know that t * (the time giving maximum value F max of delay kernel F (t)) is increasing with respect to increasing of c 2 (> 0.5), but F max is decreasing with respect to c 2 (> 0.5). This will give an optimalĉ 2 (> 0.5) to ensure the stability region maximum.
Note that we only considered the model with the special delay kernel F (t) = c 1 F 1 (t)+c 2 F 2 (t), (c 1 , c 2 ≥ 0, c 1 +c 2 = 1). In general, the delay kernel can be defined as F (t) = n j=1 c j F j (t), where c j ≥ 0, F j (t) = b j (j−1)! t j−1 e −bt , b > 0, j = 1, 2, .., n and n j=1 c j = 1. Finally, we have to point out that how to analyze the original system (3) is still an open problem. Besides, in this study, we only consider the anti-tumor immunity factors, such as ECs and HTCs. However, the pro-tumor immunity factors, such as Th17 cells and regulatory T lymphocytes are also present in the tumor microenvironment [35]. So it is important to consider the dual role of the immune response: tumor-inhibiting and tumor-promoting in the modeling of tumor immune system interaction [2].
Appendix A. Properties of delay kernel F (t). By calculating the derivative of F (t) with respect to t, we have dF (t) Then we have dF (t) dt | t=t * = 0. From (10), we can see that t * is a monotonic increasing function of c 2 , while it is a monotonic decreasing function of b. There exists a critical threshold c 2 = 0.5 that can determine the sign of t * . If c 2 ≤ 0.5, we have t * ≤ 0, further if t ≥ 0, then we have t ≥ t * and dF (t) dt ≤ 0. So F (t) is a monotonic decreasing function for all t ≥ 0. If c 2 > 0.5, we have t * > 0, and when t < t * we have dF (t) dt > 0, and when t > t * we have dF (t) dt < 0. So F (t) is a nonmonotonic function with a peak at t = t * . For relatively small c 2 (< 0.5), the curve shape of function F (t) is very similar to function F 1 (t), which means that F (t) is dominated  by F 1 (t) (see Figure 6(a)), while for relatively large c 2 (> 0.5), the curve shape of function F (t) is very similar to function F 2 (t), which means that F (t) is dominated by F 2 (t) (see Figure 6(b)). For any 0 < c 2 ≤ 1, F (t) always has the maximum value at t = t * . Denote F max = F (t * ) = c 2 be −2+ 1 c 2 .
(11) By calculating the derivative of F max with respect to c 2 , we have dFmax c 2 ≤ 0. So F max is a monotonic decreasing function of c 2 . From (11), we can see that F max is a monotonic increasing function of b. We can conclude that c 2 has an opposite role on t * and F max compared to the role of b. Large t * or small F max indicates that the system depends on the past time more than the present time. Note that F = c 1 F 1 + c 2 F 2 = (bt − 1)be −bt c 2 + be −bt . We definet = 1 b , then we have F (t) = b e . From (10), we can see that the limit of t * ist as c 2 approaches 1. When t <t, F is a monotonic decreasing function of c 2 . And when t >t, F is a monotonic increasing function of c 2 (see Figure 7). Appendix B. Summary of previous stability results of systems (2) and (5). We summarize the stability results from [36] which are used in this study. Three-dimensional system (2) has a tumor-free equilibrium E 0 = (x − , y − , z − ) = (0, σ 1 δ 2 δ 1 δ 2 − ρσ 2 , σ 2 δ 2 ).
Then we have the following result.