THRESHOLD BEHAVIOR OF AN AGE-STRUCTURED TUMOR IMMUNE MODEL ∗

. In this paper, we present and analyze an age-structured tumor immune model. Based on the fact that tumor cells of diﬀerent ages tend to exhibit diﬀerent physiological behaviors, we consider the age structure of tumor cells, age-based proliferation function and age-dependent death function in the model. The threshold R 0 for the existence of tumor-free steady state is derived. It is found that if R 0 < 1, the tumor-free steady state is not only locally stable but also globally stable. Moreover, numerical simulation shows that the threshold R 0 may be regarded as an index to reﬂect the ability of “tumor immune surveillance”, that is, the smaller the R 0 , the better the ability of tumor immune surveillance. If R 0 > 1, it is proved that the tumor steady state is existent and uniformly persistent. The local stability of the tumor steady state is investigated under some further conditions besides R 0 > 1. In the end, we estimate the system parameters, verify the theoretical results and analyze some system parameters’ sensitivities.


Introduction
As one of major diseases that endanger human health, there are many kinds of treatments for cancers (such as surgery, chemotherapy, radiotherapy, etc.), but these treatments are usually used only when tumors are detected and cause some damage to the body. The immune system is one of the most important lines of defense against tumors until humans intervene to treat them. It is particularly important to study the interrelationship between tumors and the immune system. Cells of the immune system are constantly monitoring our body tissues. When natural killer cells (NK cells) find damaged cancerous cells, they release stimulus-related molecules to activate dendritic cells (DC cells). Cytotoxic T cells release tumor-associated antigens through T-cell receptors and other co-receptors after DC cells activate them. Helper T cells help DC cells to activate cytotoxic T cells and secrete the cytokine interferon, which facilitates the recruitment and activation of more NK cells. Activated NK cells and cytotoxic T cells release relevant cytotoxic substances such as perforin, which cause gradually apoptosis of cancerous cells. The immune system monitors and removes cancerous cells through a complex mechanism. This process is called "tumor immune surveillance" [21,23].
Many scholars have proposed tumor immune related mathematical models to describe and explain the interactions between tumors and the immune system. In 1994, Vladimir et al. [16] proposed a mathematical model of cytotoxic T cells and immunogenic tumors and analyzed the global bifurcation of the model. They derived estimations of the model parameters by comparing the model with experimental data from mice. In 1998, Denise et al. [6] established a mathematical model to study the relationship between the interleukin-2 (IL-2) and the tumor immune and partially explained the effect of which to the short term oscillations in tumor sizes and the long term recurrence of tumors. In 2005, Lisette et al. [25] proposed a tumor immunity model containing NK cells and CD8 + T cells and estimated model parameters using a least squares method combined with data from mouse and human studies. In the following years, many scholars considered tumor models with treatments such as immune-modifying therapy [7], pulsed chemotherapy [30], comprehensive therapy [27,29], oncolytic virus therapy [15] and so on and studied the models dynamical behaviors. In 2021, Ruan [26] reviewed and described some findings on the dynamics of delayed differential equation (DDE) models of immune and tumor cells interactions. These works have greatly enriched the tumor immune models and provided some suggestions for tumor treatment.
Individuals with different ages in the population tend to exhibit different physiological behaviors [22]. For example, tumor cells of different ages tend to have different mortality and proliferation rates. A large number of authors have taken into account cell physiological age in tumor models. In 1987, Gyllenberg and Webb [12] proposed and analyzed a framework for the continuous age-size distribution of a population that has both normal growth and quiescent individuals. They demonstrated that the cell population's age-size distribution tends to converge to a stable state. In 1997, Arion et al. [1] investigated a linear model for age-structured cell populations and obtained sufficient and necessary conditions for asynchronous exponential growth of the model (exponential growth in population size and stabilization of population age distribution).
In 2002, Dyson et al. [8] proposed a tumor model that divides age-structured cell population into proliferative and quiescent states, allowing cells to transition between the two states and newly cells enter the quiescent state with age 0. They also obtained that the population eventually shows asynchronous exponential growth and demonstrated that the final age distribution of the population is independent of the initial age distribution. In 2008, Brikci et al. [4] proposed a nonlinear model of age-structured cell populations and considered the effect of cyclins in the model. Besides, they analyzed the asymptotic behavior of this nonlinear model solutions. In 2012, Gabriel et al. [10] analyzed an age-structured tumor model for the effect of anticancer compounds and fitted age-structured proliferation functions based on experimental data. In 2018, Liu et al. [19] considered a nonlinear age-structured tumor model and analyzed the stability of the model, after which a series of works were done to refine and improve this model [17,18]. The consideration of cellular physiological age in tumor growth models has been studied extensively and it is natural to think of taking age structure into account in tumor immune models.
To our knowledge, little work has been done to consider the effect of age structure on tumor immune models. Hence, this paper considers the age structure of tumor cells in the tumor immune model and analyzes the threshold conditions for the presence of tumors. Compared to the description of tumor population growth in most tumor immune models, our age-structured tumor immune model has the following improvements.
(i) Compared with traditional tumor immune models that apply Logistic form [6,16] describing tumor populations' growth, this paper considers the age structure of tumor cell population in the tumor immune model. (ii) Age-based proliferation and death functions are considered in the model that can reflect different physiological behaviors of tumor cells at different ages [10,19]. (iii) The tumor-free threshold R 0 is proposed. Moreover, the correlation between R 0 and system parameters is appropriately illustrated by parameter sensitivity analysis.
The article following this paper is structured as follows. In Section 2, we present the basic age-structured tumor immune model. In Section 3, some basic properties and dynamics behavior of the model are proved. In Section 4, the parameters are estimated, the theoretical results are verified by numerical simulations and the sensitivity of parameters is analyzed. In Section 5, we give a short discussion.

Age-structured tumor immune model
We construct partial and ordinary differential equations to characterize the interaction of an age-structured tumor population with immune cells in vivo. The parameters associated with tumor population depend on cellular physiological ages. Newly dividing tumor cells enter directly into the tumor population with age 0. Tumor immune surveillance mechanisms are complex and often involve the interaction of multiple immune cells with tumor cells. For simplicity, we only consider effector cells that can recognize and clear tumor cells. Assuming that immune cells entering the tumor area are mature and insensitive to cell age, we only consider the age structure of tumor cells.
Let n(t, a) represent the number of tumor cells at time t with age a. Let E(t) represent the number of effector cells at time t. We consider the following age-structured tumor immune model where N (t) = a + 0 n(t, a)da represents the total number of tumor cells at time t. The function β(·) and µ(·) represent the age-based proliferation and mortality rates of tumor cells, respectively. The function σ(·) describes the clearance of tumor cells at different ages by effector cells. The parameter a + is the maximum survival age of tumor cells and it can be simply assumed from the biological significance that n(t, a) = 0 when a > a + . Mature effector cells enter into the tumor area at a constant rate s (assuming that s is not affected by the existence of tumor cells). The positive constant d E represents the mortality rate of effector cells. Q(N, E) represents the promotion effect of the existence of tumor cells on the entry of effector cells into the tumor area. There are two usual forms: pN (t)E(t) [9] and pN (t)E(t) g+N (t) [16] in which g and p are both positive constants. The function κ(·) represents the inactivation rate of immune cells in response to interaction with tumor cells of different ages.
In order to better investigate the dynamical properties of system (2.1)-(2.2), we make some necessary assumptions and simplifications to the model. Assumption 2.1. Throughout this paper, we always let the following assumptions hold.
(A2) Assuming that the promotion rate of tumor cells to effector cells into the tumor area is greater than the inactivation rate of tumor cells to effector cells, the last two terms of the second equation of (2.1) are simplified to τ N (t)E(t) in which τ is a positive constant.
Under the above assumptions, we simplify system (2.1)−(2.2) to the following model

Basic model properties and model dynamics behavior
In this section, we study the basic properties and dynamics behavior of the model. We first discuss the boundedness of system (2.3)−(2.4) and then prove the existence and uniqueness of system solutions by rewriting it to an abstract Cauchy problem.

Boundedness, existence and uniqueness of system solutions
Before proving the boundedness, existence and uniqueness of solutions of system (2.3)−(2.4), we firstly show that the solutions are non-negative.
From the second equation of system (2.3), we know that holds for all E(0) ≥ 0 and t ≥ 0. Applying the characteristic line method to the second equation of system (2.3), we easily obtain Therefore, n(t, a) ≥ 0 holds for all n(0, a) ≥ 0, n(t, 0) ≥ 0 and t ≥ 0.  where the function space Y = L 1 + (0, a + ) × R + . Integrating the first equation of system (2.3), and noting the fact that n(t, a + ) = 0, we get (3.1) For system (3.1), we have the following theorem.

Local and global stabilities of the tumor-free steady state
In this subsection, we study the equilibrium state equations of the system and then study of the local and global stabilities of tumor-free steady state. System (2.3)−(2.4) always admits a tumor-free steady state I 0 = (n 0 (a), E 0 ) with n 0 (a) = 0 and E 0 = s d E . However, the tumor steady state I * = (n * (a), E * ) must satisfy the following time-independent system where N * = a + 0 n * (a)da. Solving the second equation of (3.2), we have where d E − τ N * > 0 has to be satisfied, whereas if not, it contradicts with s > 0. Calculating the first equation of (3.2) yields

Substituting (3.4) into the third equation of (3.2) yields
Because of n * (0) > 0 (n * (0) = 0 leads to the tumor-free steady state), from (3.5) we get needs to be satisfied when the tumor steady state I * exists. Similar to the basic reproduction number in the epidemic model, we define the threshold By a simple calculation we can also obtain that tumor-free steady state satisfies the following equation If R 0 < 1, we can obtain that (3.8) holds if and only if n 0 (0) = 0, that is, there exists a tumor-free steady state. For the special case R 0 = 1 we are able to get If n 0 (0) > 0 we can obtain where N 0 = a 0 n 0 (s)ds > 0. Equation (3.9) contradicts with R 0 = 1. Therefore, n 0 (0) = 0 when R 0 = 1. We obtain the following theorem.
Theorem 3.4. When the time is long enough, the tumor-free stable state must exist if R 0 ≤ 1 and the tumor steady state I * presences must lead to R 0 > 1.
For the local stability of the tumor-free steady state, we give the following theorem.
The proof of Theorem 3.5 will be developed in Appendix C. From Theorem 3.5, one can see that I 0 is locally asymptotically stable if R 0 < 1. Next we will show that I 0 is also globally asymptotically stable if R 0 < 1.
From the second equation of (2.3) we know that It follows that Therefore, for any small enough positive constant ε > 0, there with boundary and initial conditions As in [19], based on the principle of comparison, we get for all t > T 0 , where n(t, ·) and E(t) are integral solutions of (2.3)−(2.4), n 2 (t, ·) and E 2 (t) are integral solutions of (3.11)−(3.13). Using the same approach as the proof in Theorem 3.5, we can obtain that when R 0 < 1. Solving the second equation for (3.11) we get Combining (3.10), (3.14) and (3.16) we have Same as equation (3.7), for N 2 we define Noticing that R ε 0 < R 0 < 1 is satisfied when ε is small enough. Thus for the global stability of the tumor-free steady state I 0 we have the following theorem.

Uniform persistence
In this section we will show that system (2.3)−(2.4) is uniformly persistent and that the tumor is not cured We know that the total number of tumor cells is expressed as N = Similarly, for any ((0, n), E) ∈ ∂M 0 , we have  [13,20], the following theorem can be obtained.
Proof. We prove the uniform persistence of the system by applying the Theorem 4.1 1 in [13]. Conditions (i) − (iii) of Theorem 4.1 in [13] are clearly true. In order to prove that S(t) t≥0 is uniformly persistent, it is necessary to show that it is satisfied by Since R 0 > 1 holds, there exists a positive constant ε * small enough such that Make the contradictory assumption that there exists (3.20) Then ς(a) ≥ 1 and dς(a) da = µ(a) + β(a) + ( is a non-decreasing function for t large enough. Therefore, J(t) ≥ 0 for all t ≥ max{T 0 , T 1 }, which prevents n(t, a) converging to 0 L 1 as t → ∞. A contradiction with ζ ∈ W s ({X 0 }) occurs. This completes the proof.

Existence, uniqueness and local stability of the tumor steady state
In this subsection, we discuss the existence, uniqueness and local stability of the tumor steady state I * with R 0 > 1.
Theorem 3.7 tells us that the presence of tumor steady state I * must lead to R 0 > 1. Now we show that if R 0 > 1, the tumor steady state I * is existent and unique.
It is obvious from (3.2) that n * (a) ∈ C 1 [0, a + ]. Also, it is easy to see that n * (a) is a solution of system (3.2) if and only if it is a solution of equation (3.4). We use the Schauder's fixed point theorem to study the existence of solutions of integral equation (3.4). The lemma is as follows. Next we show that T : S → T (S) is continuous. For any small positive constant ε > 0 and arbitrary ν 1 , ν 2 ∈ S, ν 1 − ν 2 < ε/(2B N + 1) we have   The proof of Theorem 3.10 will be developed in Appendix D.
Remark 3.11. It is noticed that the characteristic equation (D.9) has no monotonicity on λ. Thus, it is difficult to derive explicit sufficient conditions mathematically on the locally asymptotical stability for the tumor steady state. We illustrate the existence (Fig. 4) and stability (Figs. 5 and 6) of the tumor steady state numerically in Section 4, which show that the system really has locally asymptotically stable tumor steady state if coefficients of the system have taken appropriate values.

Numerical simulations and parameter sensitivity analysis
In this section, we estimate the system parameters, verify the obtained theorems and analyze numerically the sensitivities of some system parameters.

Parameters estimation
The choice of appropriate system parameters is important because they determine the dynamical behaviours of the system (types of steady states and stabilities). Different tumor cells have different maximum cell survival ages [2,3,10]. For convenience we selected the maximum cell survival age as a + = 72 hours for all examples. For tumour cells, there are three important age-based coefficients: the proliferation function β(a), the mortality function µ(a) and the tumour clearance function σ(a) by immune cells.
Gabriel et al. [10] fitted several tumor cell proliferation functions β(a) with mouse PC-9 lung cancer cell data. The exponentially modified Gaussian (EMG) can well reflect the intermitotic time distribution of tumor cell mitosis [11]. There is no significant difference by the numerical comparison between the proliferation function fitted by EMG and the proliferation function fitted directly by the error equation (see Fig. 4 in [10]), which shows that the proliferation function fitted with the error equation also reflects the real tumor cells proliferation behavior well. We describe the following proliferation function in terms of an error function (as Eq. (20) in [10]) for all a ∈ [0, a + ] (see [19]). The clearance function is small in order of magnitude [16]. We ignore the effect of age on the clearance function in numerical simulations and assume that the age-based clearance function σ(a) is a constant, that is, for all a ∈ [0, a + ]. β 0 and σ 0 are positive numbers. The forms of β, µ and σ satisfy the biological significance and the assumptions of this paper. For immune cells, there are three important parameters: the constant rate s of effector cells entry into tumor area, the effector cells mortality rate d E and the coefficient τ of the tumor cells to effector cells promotion term. The references and the baseline values of the estimated system parameters are shown in Table 1.
In the following work we will study the effect of the parameters in Table 1 on the system.

Numerical simulations
In this subsection, we perform numerical simulations to study the effect of parameters on R 0 and to verify the asymptotic stability of the tumor-free steady state I 0 as well as the existence of tumor steady state I * . Thereafter, we analyze the sensitivity of the system parameters.
From Theorems 3.4 and 3.9, we know that the kind of steady state depends on values of R 0 . Equations (3.7) and (4.1)−(4.3) tell us that R 0 is related to the values of parameters µ 0 , β 0 , s, d E and σ 0 . We first investigate the effect on R 0 from changes of parameters µ 0 and β 0 (see Fig. 1). Figure 1a shows the effect of the parameter µ 0 on R 0 when β 0 takes values 0.15, 0.20 and 0.25, respectively. Figure 1b shows the effect of µ 0 and β 0 on R 0 . From Figure 1, we know that the value of R 0 increases as β 0 increases and decreases as µ 0 increases. A smaller β 0 and a larger µ 0 are needed to guarantee the condition R 0 < 1 of tumor-free steady state I 0 . Fixed µ 0 , we know that the larger the β 0 the larger the R 0 , which indicates that a larger β 0 leads to a higher the malignancy of the tumor.
It follows from the expression of E 0 and equation (3.7) that the effect of parameters s, d E and σ 0 on R 0 can be simplified to the effect of the number of tumor-free steady state immune cells E 0 and σ 0 on R 0 . Figure 2 demonstrates the effect of E 0 and σ 0 on R 0 . From Figure 2, it is known that a larger E 0 and a larger σ 0 are required to satisfy the condition R 0 < 1 for the tumor-free steady state. This means that the higher the body's  Table 1.  Table 1.
immune cell content and the greater the ability to recognize and clear tumor cells, the better the tumor immune surveillance ability. We study the global stability of the tumor-free steady state I 0 when the condition R 0 < 1 is satisfied. Let age-dependent tumor cells initial value function n(0, a) = 2 × 10 5 × (a + − a)/a + and effector cell population initial value E 0 = 3.2 × 10 5 . Let µ 0 = 0.003 and other parameters be chosen the same as which in Table 1. It is easy to calculate that R 0 = 0.9339 < 1, which satisfies the condition of Theorem 3.6. Thus, the tumor-free steady state I 0 of system (2.3)−(2.4) is globally asymptotically stable. Figure 3a shows the age distribution of tumor cells when R 0 < 1. Figure 3b  the system requires approximately no less than 7000 hours to reach the tumor-free steady state. This means that even if R 0 < 1 is satisfied, the immune system needs a time long enough to clear the tumor cells when the tumor occurs. Next, we illustrate the existence and the local stability of the tumor steady state I * .
Let µ = 0.001 and the other parameters be chosen as those in Table 1. It is easy to calculate that R 0 = 1.0747 > 1 and thus the positive steady state exists from Theorem 3.9. Choose the tumor population initial value function n(0, a) = 2 × 10 5 × (a + − a)/a + and the effector cell initial value E 0 = 3.2 × 10 5 . Figure 4a shows the age distribution of tumor cells. Figure 4b plots the time series of effector cells E. It can be seen that both the effector cells and the age-dependent tumor cells will tend to a steady state after a long time. Moreover, the fact that E * > E 0 > 0 can also be known in Figures 3b and 4b.
We take different initial values to illustrate the local stability of the tumor steady state. Figures 5a and 6a show the time series of tumor cells with age 0 under three different tumor cell initial value functions and three different effector cell initial values respectively. Figures 5b and 6b plot the curves of the total number of tumor cells N under three different tumor cell initial value functions and three different effector cell initial values respectively. Figures 5c and 6c demonstrate tendencies of number of effector cells E under different initial values as time goes. As can be seen in Figures 5 and 6 for different initial values, the number of newborn tumor cells n(t, 0) tend to stabilize to a positive constant and the number of tumor cells N (t) as well as effector cells E(t) converge to the same value over a long period of time, which implies that the positive steady state I * of the system is unique and stable. The steady state of the system is independent of the initial values. The initial value only affects the time the system reaches the steady state.

Parameter sensitivity analysis
It is of importance to determine that which parameters have significant effects on growth of the number of tumor cells. We investigate effects of parameter changes on the total number of tumor cells N and effector cells E at a sufficiently large time point (t = 16000 hs). We numerically analyze one parameter's sensitivity by assuming the other parameters are unchanged.
From Figures 5 and 6, we have known that system (2.3)−(2.4) tends to steady state after t = 16000hs. Figure 7 shows successively the effects of changes of parameters µ 0 , β 0 , σ 0 , τ , s and d E on the total number of   Figure 7a plots changes of the number of tumor cells N and the number of effector cells E with mortality rate of tumor cells µ 0 goes. One can see that as µ 0 increases, N and E will decrease. N decreases significantly and forces the system shifts its steady state from tumor steady state to tumor-free steady state as µ 0 increases. Figure 7b shows variation tendencies of N and E as the proliferation-dependent coefficient β 0 changes. Obviously, an increase of parameter β 0 contributes to an increase of N and E. But the increase intensity of the effector cells E is not significant compared to the total tumor cells N as β 0 increases. Figure 7c demonstrates that an increase of tumor cells clearance function σ 0 will lead to number decrease of both tumor cells and effector cells. Figure 7d presents that an increase of promotion coefficient of tumor cells to effector cells τ will cause number decrease of tumor cells significantly but number changes of effector cells unobviously. Figure 7e shows changes of N and E as rate of mature effector cells entering into the tumor area s changes. It can be seen that the increase of s will cause N 's decrease but E's increase. Figure 7f illustrates that the larger the mortality rate of effector cells d E the more the tumor cells but the larger the d E the less the effector cells.
The numerical sensitivity analysis on parameters µ 0 , β 0 , σ 0 , τ , s and d E yields that all parameters above except τ can determine the types of system's steady state, which is consistent with the expression of the tumor-free steady state threshold R 0 , that is, One can see from figures that although parameters µ 0 , β 0 , σ 0 , s and d E can all affect the steady state population numbers and steady state types, the sensitivities of parameters are different. Obviously, the system is more sensitive to parameter σ 0 . Moreover, we find that number changes of N and E are very sensitive to changes of parameter τ , though who does not affect the kind of the system steady state.

Conclusion and discussion
We proposed and analyzed an age-structured tumor immune model. First, reasonable assumptions and simplifications were made on the model, and then the existence, uniqueness and ultimate boundedness of the solutions of the simplified system were proved. The threshold condition R 0 < 1 for the existence of the tumor-free steady state was obtained, after which the local and global stabilities of the tumor-free steady state were investigated. The existence and uniqueness of the tumor steady state was also proved, and it was known that the tumor steady state is uniformly persistent when the condition R 0 > 1 is satisfied. The local stability of the tumor steady state was proved under further conditions. We verified the theoretical results and analyzed the sensitivities of some system parameters by numerical simulations, knowing that the system is more sensitive to parameters σ 0 and τ . The parameter settings in the present work are from literatures. If real human experimental data can be acquired, we may be able to estimate the parameter values better and improve the applicability of this paper. This may be an important direction for the future work.
This paper only considered interactions between tumor cells and effector cells. It is more relevant to consider the role of different therapeutic approaches when diagnosing a tumor. Moreover, spatial factors, stochastic perturbations and different objects involved in the tumor micro-environment in the modelling process are not ignored and will be the challenging problems in the future work.
From the second equation of system (3.1), we know that Therefore, there exists T 0 > 0 and 0 < 0 < min{s/d E , 1/(2σ)} such that Thus, there exist T 1 ≥ T 0 and 0 <N < d E /τ such that N (t) ≤N holds for all t ≥ T 1 . Then, for all t ≥ T 1 , we have which implies that holds. According to (A.5) and (A.6), we know that the system (3.1) is ultimately bounded. Fig. 8 . Clearly, 0 < N m2 < N m1 . We will show the ultimate boundedness of solutions of system (3.1) in two steps.
Step I: We firstly show that all solutions starting from regions Γ 1 , Γ 2 and Γ 3 must cross N = N m1 or E = E m and enter into region Γ 4 at least once in a finite time.
It is easy to calculate that in Γ 1 , Thus, all solutions starting from Γ 1 must cross N = N m1 and enter into Γ 4 at least once or cross E = E m and enter into Γ 2 in a finite time. Similarly, all solutions starting from Γ 2 must cross E = E m and enter into Γ 4 at least once or cross N = N m2 and enter into Γ 3 in a finite time. And, all the solutions starting from Γ 3 must cross E = E m and enter into Γ 4 in a finite time.
Step II: We now prove that all solutions that start from region Γ 4 must have an upper bound. It is obvious that solutions starting from Γ 4 are bounded if they always stay in this region.
For the solutions starting from Γ 4 that cross N = N m1 and enter into Γ 1 , we are sure that there is a positive constantŇ such that N (t) ≤Ň for all t ≥ 0. In fact, it follows from dE/dt ≥ s in Γ 1 that the longest time the solution starting from Γ 4 stays in Γ 1 for one time is no more than E m /s. Therefore, to these solutions, the maximum of N is no larger than N m1 · eβ Em/s because of the inequality dN/dt ≤βN . DenoteŇ = N m1 · eβ Em/s , since dN/dt ≤ −N m2 < 0 for all solutions that in Γ 2 , we can conclude that all solutions starting from Γ 4 must have an upper bound in direction of N . That is, N (t) ≤Ň for all t ≥ 0 and (N (0), E(0)) taken from Γ 4 .
Similarly, we can calculate that the time that solution starting from Γ 4 stays in Γ 2 for one time is no longer thanŇ /N m2 and then the maximum of E in Γ 2 is no larger than E m · e (d E +τŇ )Ň /Nm2 :=Ě. And then due to the fact dE/dt ≤ −δ 0 in Γ 3 , the maximum of E in Γ 3 is also no larger thanĚ for all t ≥ 0 and (N (0), E(0)) taken from Γ 4 .
Finally, we have the result that there are two positive constantsŇ andĚ such that lim sup t→∞ N (t) ≤Ň and lim sup t→∞ E(t) ≤Ě for any initial value N (0) = N 0 ≥ 0 and E(0) = E 0 ≥ 0. This completes the proof.
From the definition of the linear operator A and applying the Hille-Yosida Theorem, we know that A is an infinitesimal generator of a C 0 − semigroup {T (t)} t>0 = {e At } t>0 . From Assumption 2.1(i) we know that the nonlinear operator F is Lipschitz continuous. By applying the results given in [5,13,14], we obtain the existence and uniqueness of solutions of the system (2.3)−(2.4).
Appendix C. Proof of local stability of the tumor-free steady state Proof. Let n 1 (t, a) = n(t, a), E 1 (t) = E(t) − E 0 . We can calculate that the linearized equations of (2.3)−(2.4) at I 0 is 1 (t, a), with boundary and initial conditions  H(λ). It is clear that H (λ) < 0, H(0) = R 0 < 1 and lim (λ)→∞ H(λ) = 0 ( (λ) represents the real part of λ). Thus we can know that H(λ) is a monotonically decreasing function and equation (C.8) has no positive real roots if R 0 < 1. If any complex root has nonnegative real part, the modulus on the right hand side of equation (C.8) is always less than 1 when R 0 < 1, which contradicts to equation (C.8) . Therefore, all roots of (C.8) have negative real parts when R 0 < 1. Solving the second equation of (C.1) yields