Modeling and analysis of a multilayer solid tumour with cell physiological age and resource limitations

Abstract

We study an avascular spherical solid tumour model with cell physiological age and resource constraints in vivo. We divide the tumour cells into three components: proliferating cells, quiescent cells and dead cells in necrotic core. We assume that the division rate of proliferating cells is nonlinear due to the nutritional and spatial constraints. The proportion of newborn tumour cells entering directly into quiescent state is considered, since this proportion can respond to the therapeutic effect of drug. We establish a nonlinear age-structured tumour cell population model. We investigate the existence and uniqueness of the model solution and explore the local and global stabilities of the tumour-free steady state. The existence and local stability of the tumour steady state are studied. Finally, some numerical simulations are performed to verify the theoretical results and to investigate the effects of different parameters on the model.

KEYWORDS:

• Nonlinear age-structured model
• proliferating cells
• resource limitation
• steady state

MATHEMATICS SUBJECT CLASSIFICATION:

• 92B05

1. Introduction

Tumor growth model originated from the empirical equation was first proposed by Gompertz in 1825 to simulate tumour growth [1], it remains highly relevant in a wide range of biomedical fields due to its good fit and simplicity with experimental data [2,3]. Later, a large number of papers aimed at understanding and developing the mechanisms of Gompertz's model equation [4–7]. Over the next numerous years, a large number of mathematical models were proposed and studied for tumours in vivo or in vitro [8–13]. In which Gyllenberg and Webb's model divided cancer cell populations into proliferating and quiescent cells, which greatly contributed to the development of mathematical models of tumours [11]. Single-population tumour models provide a concise description of tumour growth, while multiple-population tumour models provide a more detailed description of the different physiological behaviours of different types of subpopulations within a tumour population.

Population age has been widely considered in various fields of bio-mathematics. Population models that consider age structure better reflects the real world, and corresponding models such as demographic models [14,15], epidemic models [16–18], microscopic virus models [19,20] and cell population models [8,21] are presented. In the cell population model, cells are divided by age to represent the period of cell cycle. In 2002, Dyson et al. [22] took the age structure into account in a linear two-component cell model, and he divided cell populations into proliferative and quiescent states to describe the different proliferative behaviours of cells. His model well characterized the early exponential growth of the cells and demonstrated that the final age distribution tended to be consistent.

Solid tumours tend to exhibit different growth rates in vascular and avascular growth environments. Initially, solid tumours in the avascular environment usually exhibit exponential growth based on the relative abundance of nutrients and space, followed by the development of a linear growth stage due to areas of quiescent cell development and core necrosis [23–25], In the middle to later stage avascular tumour growth exhibits nonlinear proliferation. Ward et al. [26] constructed a model of avascular tumour growth and investigated the exponential (initial) and linear (intermediate) stages of growth. Based on the improvement of Gyllenberg and Webb's model [11] on the assumption of immediate clearance of dead cells, Alzahrani et al. [27] further extended the model to a three-compartment form, after which the above model was further improved by considering a resource-constrained form [28], and by studying the steady state of the tumour, the relationship between tumour size and the clearance rate of dead cells at positive steady state was found. Their work takes into account the prevalence of necrotic nuclear properties in advanced avascular tumours. Necrotic cores are mostly found in advanced avascular tumours, and the multilayer solid tumour model is more suitable for middle to later stage spherical tumour growth with avascular and resource-limited environment in vivo.

If solid tumour cells growing in an avascular environment do not produce enough vascular growth signalling protein (e.g. angiopoietin), then the tumour will only grow to the maximum size supported by the avascular environment and will not enter a vascular state. In avascular multicellular tumour spheroids, proliferating cells usually enter a quiescent state and stop dividing due to resource and space constraints, and eventually quiescent cells die and form necrotic cores inside the tumour. Avascular tumours do not immediately remove the dead cells from the necrotic cores, forming multilayered solid tumour structures. The quiescent layer and necrotic core tend to be common in the later stages of avascular multicellular tumour spheroids [29]. Quiescent cells usually increase with tumour size, usually in a nonlinear fashion. The proliferation rate of tumour proliferating cells tends to exhibit nonlinear characteristics in the middle and late stages of avascular spheroid tumours. Liu et al. [30–32] considered nonlinear proliferation functions in age-structured tumour models and analysed the stability of tumour models with nonlinear proliferation functions. Due to the nonlinear effects of resource limitation on cell populations in an avascular tumour environment, nonlinear models can better describe tumour growth in this situation. We consider nonlinear proliferation functions to the cell population model as in [22,33,34] and consider necrotic cores in tumour cells to represent the multilayered structure of solid spherical tumours in vivo. The physiological age change of tumour cells may differ from the real time change, and we emphasize the evolution speed κ of tumour cells in the cell cycle ($\kappa =$ tumour cell age/time) [22,30,33–35]. The ratio of physiological age to time is different for different tumour cells. For a certain class of tumours, benign tumours may have a slower physiological age growth and thus have a smaller value of κ. Cycle-specific drugs are the most commonly used anti-cancer drugs, and we describe the therapeutic effect of cycle-specific drugs in terms of the proportion ℓ of newborn tumour cells entering the quiescent state [30,34].

The rest of this paper is structured as follows. In Section 2, we present a three-component tumour model with nonlinear age-structured proliferation functions. In Section 3, we rewrite the model and establish conditions for the existence, uniqueness and boundedness of the model. The existence and stability of tumour-free stable state is discussed in Section 4. The existence and local stability of the tumour steady state is studied in Section 5. In Section 6, we verify the results with numerical simulations and analyse the effects of parameter changes on the subpopulation size of the tumour. We give a short discussion at the end of the paper.

2. Multilayer age-structured tumour model

We construct a three-layer solid tumour model with nonlinear age-structured proliferation functions. We divide the tumour cells into three components being proliferating cells, quiescent cells and dead cells in necrotic cores not removed in time. Due to the three-layer structure of the avascular multicellular tumour spheroids in the late phase, we assume that death proliferating cells in the outer layer is cleared immediately, while dead quiescent cells enter the necrotic core and are cleared at a clearance rate d. The key developments in this model are based on previous research including [8,22,30,34]. Proliferating cells exhibit different physiological behaviours at different ages, whereas quiescent cells are not sensitive to their age stage, so in our model we only consider the age structure of proliferating cells.

The ratio of the change in physiological age of proliferating cells to the real time change is assumed to be a constant value κ, i.e. $\kappa =\mathrm{d}a/\mathrm{d}t$. If $\kappa =2$, it means that the physiological age a change is twice the real time t change, which may indicate more active tumour cells growth. We define the nonlinear division rate of proliferating tumour cells as $b(a,N(t))$, which is related to the age of proliferating tumour cells a and the total number of tumour cells $N(t)$. The proliferating cell division rate $b(a,N(t))$ decreases as the population size $N(t)$ increases. Tumor cells need time to prepare substances required for cell proliferation, meaning that cell proliferation behaviour is correlated with cell age a. For the relatively limited amount of material in the avascular tumour environment, the larger the total number of tumour cells N, the less material is available for one tumour cell proliferation, so an increase for total number of tumour cells will lead to a decrease for cell proliferation function. We also assume that the age of the newly proliferating cells is 0. The newborn cells enter the proliferative phase again in proportion $1-\ell$ and enter the quiescent phase directly in proportion ℓ with the age of 0 [8,33,35]. This parameter ℓ reflects the effect of the therapeutic dose of drugs such as erlotinib [36], which may affects the proportion of newborn cancer cells entering into quiescence [34].

Since proliferating cells exhibit different physiological states at different ages, we consider the age-based change in the death rate of proliferating cells defined as ${\mu }_p(a)$ at age a and the death rate of quiescent cells defined as a constant ${\mu }_q>0$. The constant d>0 is defined as the clearance rate of dead cells in necrotic care. As Liu [30], we consider $[0,a^{+}]$ as the cell age interval of cancer. We assume that due to spatial and nutritional constraints there is no quiescent cells entry into proliferating cells, but has age-based transition from proliferating cells to quiescent cells. This transfer rate is defined as $\vartheta (a)$ at age a with $a\in [0,a^{+}]$.

The model adopts the following three-compartment system (1) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{d}Q}{\mathrm{d}t}=2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a-{\mu }_{q}Q(t),\\ \frac{\mathrm{d}D}{\mathrm{d}t}={\mu }_{q}Q(t)-\mathrm{d}D(t)\end{array}$(1) with boundary and initial conditions (2) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ P(0,a)=P_0(a)\ge 0,\\ Q(0)=Q_0\ge 0,\\ D(0)=D_0\ge 0,\end{array}$(2) where $Q(t),D(t)$ are numbers of quiescent and dead cells at time t, $P(t,a)$ is the age-specific density of proliferating cell population at time t with age a and age range of $[0,a^{+}]$. ${\int }_{a_1}^{a_2}P(t,a)\mathrm{d}a$ represents the total number of proliferating cells physiological age from $a_1$ to $a_2$ at time t. $N(t)={\int }_0^{a^{+}}P(t,a)\mathrm{d}a+Q(t)+D(t)$ is the total number of cells at time t including proliferating ones, quiescent ones and dead ones.

For the rest of this this paper, we let the parameters ${\mu }_q,d,\kappa ,\ell$ are positive and assumption (H1) hold.

(H1) $b(\cdot ,N)$, ${\mu }_p(\cdot )$, $\vartheta (\cdot )$ $\in L_{+}^{\mathrm{\infty }}[0,a^{+}]$, where $N\ge 0$.

The nonnegative assumption on the age-specific parameters ${\mu }_p,b$ and ϑ satisfies the biological significance. (H1) guarantees the existence and uniqueness of the model solutions.

Following this, we first do some preparatory work and prove some fundamental properties of the model.

3. Model transformation and fundamental properties

In this section, we rewrite the model to a convenient form to prove the existence, uniqueness and boundedness of the solution of the model.

3.1. Model transformation

Let (3) $Q(t)={\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a,D(t)={\int }_0^{+\mathrm{\infty }}D(t,a)\mathrm{d}a.$(3) The total number of proliferation cells at time t is $P(t)={\int }_0^{+\mathrm{\infty }}P(t,a)\mathrm{d}a$. In this way we make the proliferation cells, the quiescent cells and the dead cells have the same structures. According to the biological significance of $P(t,a)$ we know that $P(t,a)=0$ when $a\in (a^{+},+\mathrm{\infty })$. Similarly, we assume that Q and D satisfy $Q(t,a)=0$ and $D(t,a)=0$ when $a\in (a^{+},+\mathrm{\infty })$.

Based on (3(3) $Q(t)={\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a,D(t)={\int }_0^{+\mathrm{\infty }}D(t,a)\mathrm{d}a.$(3) ), we rewrite system (1(1) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{d}Q}{\mathrm{d}t}=2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a-{\mu }_{q}Q(t),\\ \frac{\mathrm{d}D}{\mathrm{d}t}={\mu }_{q}Q(t)-\mathrm{d}D(t)\end{array}$(1) )–( 2(2) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ P(0,a)=P_0(a)\ge 0,\\ Q(0)=Q_0\ge 0,\\ D(0)=D_0\ge 0,\end{array}$(2) ) to be (4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) with initial conditions (5) $\{\begin{array}{l}P(0,a)=P_0(a)\in L_{+}^1((0,+\mathrm{\infty }),\mathbb{R}),\\ Q(0,a)=Q_0(a)\in L_{+}^1((0,+\mathrm{\infty }),\mathbb{R}),\\ D(0,a)=D_0(a)\in L_{+}^1((0,+\mathrm{\infty }),\mathbb{R})\end{array}$(5) and boundary conditions (6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) Obviously, the existence, uniqueness and boundedness of solutions between system (1(1) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{d}Q}{\mathrm{d}t}=2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a-{\mu }_{q}Q(t),\\ \frac{\mathrm{d}D}{\mathrm{d}t}={\mu }_{q}Q(t)-\mathrm{d}D(t)\end{array}$(1) )–( 2(2) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ P(0,a)=P_0(a)\ge 0,\\ Q(0)=Q_0\ge 0,\\ D(0)=D_0\ge 0,\end{array}$(2) ) and new system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) are consistent. We then investigate the existence, uniqueness and boundedness of solutions of system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ).

3.2. Existence and uniqueness of system solutions

We study the fundamental properties of system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) including the existence and uniqueness of the system solutions. Define a Banach space $X=L^1(0,a^{+})\times L^1(0,a^{+})\times L^1(0,a^{+})$ with norm $\Vert \varphi \Vert =\Vert {\varphi }_1\Vert +\Vert {\varphi }_2\Vert +\Vert {\varphi }_3\Vert$ for $\varphi (a)=({\varphi }_1(a),{\varphi }_2(a),{\varphi }_3(a){)}^T\in X$, where $\Vert \cdot \Vert$ is the norm in $L^1$ and ${\nu }^T$ is the transpose of the vector ν. Consider the initial-boundary value problem of system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) as an abstract Cauchy problem on the Banach space X.

Let A be a linear operator defined by (7) $(A\varphi )(a):={\left(\begin{array}{ccc}-\kappa \frac{\mathrm{d}}{\mathrm{d}a}{\varphi }_1(a),& -\kappa \frac{\mathrm{d}}{\mathrm{d}a}{\varphi }_2(a)-{\mu }_p{\varphi }_2(a),& \kappa \frac{\mathrm{d}}{\mathrm{d}a}{\varphi }_3(a)-\mathrm{d}{\varphi }_3(a)\end{array}\right)}^T.$(7) The domain $D(A)\in X\to X$ is given as follows (8) $\begin{array}{rl}D(A)& =\{\varphi \in X_{+}:=L_{+}^1(0,a^{+})\times L_{+}^1(0,a^{+})\times L_{+}^1(0,a^{+}):\\ & {\varphi }_1,{\varphi }_2,{\varphi }_3\in AC[0,a^{+}],\varphi (0)=({\varphi }_1(0),{\varphi }_2(0),{\varphi }_3(0){)}^T\},\end{array}$(8) where $AC[0,a^{+}]$ denotes the set of absolutely continuous functions on interval $[0,a^{+}]$, $L_{+}^1(0,a^{+})$ is positive cone of $L^1(0,a^{+})$, (9) $\begin{array}{rl}{\varphi }_1(0)& =2(1-\ell ){\int }_0^{a^{+}}b(a,N(t)){\varphi }_1(a)\mathrm{d}a,\\ {\varphi }_2(0)& =\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t)){\varphi }_1(a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a){\varphi }_1(a)\mathrm{d}a],\\ {\varphi }_3(0)& =\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}{\varphi }_2(a)\mathrm{d}a].\end{array}$(9) We define a nonlinear operator $F:X_{+}\to X$ as (10) $(F\varphi )(a)=\left(\begin{array}{c}-[b(a,N){\varphi }_1(a)+\vartheta (a){\varphi }_1(a)+{\mu }_p(a){\varphi }_1(a)]\\ 0\\ 0\end{array}\right).$(10) Based on assumption (H1), we can know that the operator F is Lipschitz continuous.

Let $v(t)=(P(t,\cdot ),Q(t,\cdot ),D(t,\cdot ){)}^T\in X$. The system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) can be considered as a semilinear Cauchy problem in Banach space X as follows (11) $\frac{\mathrm{d}v(t)}{\mathrm{d}t}=Av(t)+F(v(t)),v(0)=v_0\in X,$(11) where $v_0(a)=(P_0(a),Q_0(a),D_0(a){)}^T$. Using the same method as [37], the following theorem can be obtained.

Theorem 3.1

For any initial value $v_0\in X_{+}$, there exists a unique mild solution $S(t)v_0$ in $X_{+}$ for the abstract Cauchy problem (11(11) $\frac{\mathrm{d}v(t)}{\mathrm{d}t}=Av(t)+F(v(t)),v(0)=v_0\in X,$(11) ), and the positive invariant set of the semiflow $\{S(t){\}}_{t\ge 0}$ is $X_{+}$.

3.3. Boundedness of the model

Boundedness of the system is important for studying tumour models. In the following, we discuss the boundedness of system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ). We make the following assumptions about the nonlinear age-structured proliferation function.

(H2) Assume that the nonlinear proliferation rate $b(a,N(t))$ has the separated variable form $b(a)\mathrm{\Psi }(N(t))$. $\mathrm{\Psi }(x)$ satisfies $\mathrm{\Psi }(0)=1,\underset{x\to +\mathrm{\infty }}{lim}\mathrm{\Psi }(x)=0$, $\mathrm{\Psi }(x)$ is continuously differentiable in $[0,+\mathrm{\infty })$. In addition, $\frac{\mathrm{d}\mathrm{\Psi }(x)}{\mathrm{d}t}\le 0$ and $0\le \mathrm{\Psi }(x)x\le \hat{M}$ when $x\in [0,+\mathrm{\infty })$, where $\hat{M}$ is a positive constant. $b(a)$ satisfies $0\le b(a)\le \hat{b}$ with a positive constant $\hat{b}$.

Based on assumption (H2) we can get the following theorem.

Theorem 3.2

System (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) is ultimately bounded if assumptions (H1) and (H2) are satisfied.

Proof.

Let $\tilde{N}(t)={\int }_0^{+\mathrm{\infty }}P(t,a)\mathrm{d}a+Q(t)$ be the total number of all living cancer cells. From system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) we can know that (12) $\begin{array}{rl}\frac{\mathrm{d}\tilde{N}(t)}{\mathrm{d}t}& =\kappa p(t,0)-{\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a-{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a\\ & -{\int }_0^{a^{+}}{\mu }_p(a)P(t,a)\mathrm{d}a\\ & +2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a-{\mu }_{q}Q(t)\\ & =[2(\kappa -1)(1-\ell )+1]{\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a\\ & -{\int }_0^{a^{+}}{\mu }_p(a)P(t,a)\mathrm{d}a-{\mu }_{q}Q(t)\\ & \le [2(\kappa -1)(1-\ell )+1]{\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a-\hat{\mu }\tilde{N}(t),\end{array}$(12) where $\hat{\mu }=min\{{\mu }_q,\underset{0\le a\le a^{+}}{min}{\mu }_p(a)\}$. From assumption (H2) we get (13) ${\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a\le \mathrm{\Psi }(N(t))\hat{b}{\int }_0^{a^{+}}P(t,a)\mathrm{d}a\le \hat{b}\hat{M}.$(13) Let $\overline{M}=[2(\kappa -1)(1-\ell )+1]\hat{b}\hat{M}$, therefore, (14) $\underset{t\to +\mathrm{\infty }}{lim}\tilde{N}(t)\le \frac{\overline{M}}{\hat{\mu }}.$(14) Similarly, for t large enough, we have (15) $\frac{\mathrm{d}D}{\mathrm{d}t}\le {\mu }_q\frac{\overline{M}}{\hat{\mu }}-\mathrm{d}D(t),\underset{t\to +\mathrm{\infty }}{lim}D(t)\le \frac{{\mu }_q\overline{M}}{\hat{\mu }d}.$(15) In addition, (16) $\underset{t\to +\mathrm{\infty }}{lim}N(t)\le \frac{\overline{M}}{\hat{\mu }}(\frac{{\mu }_q}{d}+1):=M.$(16) Therefore,$\underset{t\to +\mathrm{\infty }}{lim}N(t)\le M$, M is a positive constant. It is easy to see system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) is bounded if assumptions (H1) and (H2) are satisfied.

4. Existence and global stability of tumour-free steady state

In this section, we study the local stability of the tumour-free steady state by using a linear system and then study the global stability of the tumour-free steady state under two different conditions.

The steady state $\overline{E}(a):=(\overline{P}(a),\overline{Q}(a),\overline{D}(a))$ of system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) needs to satisfy the time-independent system (17) $\{\begin{array}{l}\frac{\mathrm{d}\overline{P}}{\mathrm{d}a}=-\overline{b}(a,N)\overline{P}(a)-\overline{\vartheta }(a)\overline{P}(a)-{\overline{\mu }}_p(a)\overline{P}(a),\\ \frac{\mathrm{d}\overline{Q}}{\mathrm{d}a}=-{\overline{\mu }}_q\overline{Q}(a),\\ \frac{\mathrm{d}\overline{D}}{\mathrm{d}a}=-\overline{d}\overline{D}(a)\end{array}$(17) with initial conditions (18) $\begin{array}{rl}\overline{P}(0)& =2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a,\\ \overline{Q}(0)& =2\ell {\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{P}(a)\mathrm{d}a,\\ \overline{D}(0)& ={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{Q}(a)\mathrm{d}a,\end{array}$(18) where (19) $\overline{b}(a,N)=\frac{b(a,N)}{\kappa },\overline{\vartheta }(a)=\frac{\vartheta (a)}{\kappa },{\overline{\mu }}_p(a)=\frac{{\mu }_p(a)}{\kappa },\overline{d}=\frac{\mathrm{d}}{\kappa },{\overline{\mu }}_q=\frac{{\mu }_q}{\kappa }$(19) with $N={\int }_0^{a^{+}}(\overline{P}(a)+\overline{Q}(a)+\overline{D}(a))\mathrm{d}a$. The existence of tumour cell subpopulations is consistent with the initial value (18(18) $\begin{array}{rl}\overline{P}(0)& =2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a,\\ \overline{Q}(0)& =2\ell {\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{P}(a)\mathrm{d}a,\\ \overline{D}(0)& ={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{Q}(a)\mathrm{d}a,\end{array}$(18) ) when equation (17(17) $\{\begin{array}{l}\frac{\mathrm{d}\overline{P}}{\mathrm{d}a}=-\overline{b}(a,N)\overline{P}(a)-\overline{\vartheta }(a)\overline{P}(a)-{\overline{\mu }}_p(a)\overline{P}(a),\\ \frac{\mathrm{d}\overline{Q}}{\mathrm{d}a}=-{\overline{\mu }}_q\overline{Q}(a),\\ \frac{\mathrm{d}\overline{D}}{\mathrm{d}a}=-\overline{d}\overline{D}(a)\end{array}$(17) )–( 18(18) $\begin{array}{rl}\overline{P}(0)& =2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a,\\ \overline{Q}(0)& =2\ell {\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{P}(a)\mathrm{d}a,\\ \overline{D}(0)& ={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{Q}(a)\mathrm{d}a,\end{array}$(18) ) is satisfied. This means that either all three cell subpopulations are present or not present at the steady state. The system (17(17) $\{\begin{array}{l}\frac{\mathrm{d}\overline{P}}{\mathrm{d}a}=-\overline{b}(a,N)\overline{P}(a)-\overline{\vartheta }(a)\overline{P}(a)-{\overline{\mu }}_p(a)\overline{P}(a),\\ \frac{\mathrm{d}\overline{Q}}{\mathrm{d}a}=-{\overline{\mu }}_q\overline{Q}(a),\\ \frac{\mathrm{d}\overline{D}}{\mathrm{d}a}=-\overline{d}\overline{D}(a)\end{array}$(17) )–( 18(18) $\begin{array}{rl}\overline{P}(0)& =2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a,\\ \overline{Q}(0)& =2\ell {\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{P}(a)\mathrm{d}a,\\ \overline{D}(0)& ={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{Q}(a)\mathrm{d}a,\end{array}$(18) ) has a tumour-free stable state ${\overline{E}}_0:=(0,0,0)$ and no boundary stable states.

4.1. Local stability of the tumour-free steady state

To discuss the local stability of the tumour-free steady state, we hypothesize that

(H3) $b(\cdot ,N)$ is differentiable for N.

Based on assumption (H3), we can derive the tumour-free steady state linear system as (20) $\{\begin{array}{l}\frac{\mathrm{\partial }x}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }x}{\mathrm{\partial }a}=-(b(a)+\vartheta (a)+{\mu }_p(a))x(t,a),\\ \frac{\mathrm{\partial }y}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }y}{\mathrm{\partial }a}=-{\mu }_{q}y(t,a),\\ \frac{\mathrm{\partial }z}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }z}{\mathrm{\partial }a}=-\mathrm{d}z(t,a)\end{array}$(20) with the boundary conditions (21) $\{\begin{array}{l}x(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a)x(t,a)\mathrm{d}a,\\ y(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a)x(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)x(t,a)\mathrm{d}a],\\ z(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{a^{+}}y(t,a)\mathrm{d}a].\end{array}$(21) Assume that the separating variable forms of the solutions are (22) $x(t,a)=e^{\lambda t}\overline{x}(a),y(t,a)=e^{\lambda t}\overline{y}(a),z(t,a)=e^{\lambda t}\overline{z}(a).$(22) Then, we get (23) $\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=-\overline{\lambda }\overline{x}(a)-(\overline{b}(a)+\overline{\vartheta }(a)+{\overline{\mu }}_p(a))\overline{x}(a),\\ \frac{\mathrm{d}y}{\mathrm{d}t}=-\overline{\lambda }\overline{y}(a)-{\overline{\mu }}_q\overline{y}(a),\\ \frac{\mathrm{d}z}{\mathrm{d}t}=-\overline{\lambda }\overline{z}(a)-\overline{d}\overline{z}(a),\\ \overline{x}(0)=2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a)\overline{x}(a)\mathrm{d}a,\\ \overline{y}(0)=2\ell {\int }_0^{a^{+}}\overline{b}(a)\overline{x}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{x}(a)\mathrm{d}a,\\ \overline{z}(0)={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{y}(a)\mathrm{d}a.\end{array}$(23) Solving the first equation of (23(23) $\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=-\overline{\lambda }\overline{x}(a)-(\overline{b}(a)+\overline{\vartheta }(a)+{\overline{\mu }}_p(a))\overline{x}(a),\\ \frac{\mathrm{d}y}{\mathrm{d}t}=-\overline{\lambda }\overline{y}(a)-{\overline{\mu }}_q\overline{y}(a),\\ \frac{\mathrm{d}z}{\mathrm{d}t}=-\overline{\lambda }\overline{z}(a)-\overline{d}\overline{z}(a),\\ \overline{x}(0)=2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a)\overline{x}(a)\mathrm{d}a,\\ \overline{y}(0)=2\ell {\int }_0^{a^{+}}\overline{b}(a)\overline{x}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{x}(a)\mathrm{d}a,\\ \overline{z}(0)={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{y}(a)\mathrm{d}a.\end{array}$(23) ) and substitute it into the forth equation, we get (24) $\overline{x}(0)=2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a)\overline{x}(0)e^{-\overline{\lambda }a-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a,$(24) where $b(a):=b(a,N){|}_{N=0}$ and $\overline{\epsilon }(\xi ):=\overline{b}(\xi )+\overline{\vartheta }(\xi )+{\overline{\mu }}_p(\xi )$. Because of $\overline{x}(0)\ne 0$ ($\overline{x}(0)=0$ leads to the tumour-free steady state), from (24(24) $\overline{x}(0)=2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a)\overline{x}(0)e^{-\overline{\lambda }a-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a,$(24) ) we get the following characteristic equation: (25) $1=2\kappa (1-\vartheta ){\int }_0^{a^{+}}\overline{b}(a)e^{-\overline{\lambda }a-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a.$(25) Define the right-hand side of Equation (25(25) $1=2\kappa (1-\vartheta ){\int }_0^{a^{+}}\overline{b}(a)e^{-\overline{\lambda }a-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a.$(25) ) by $ϝ(\overline{\lambda })$. It is not difficult to see that $ϝ(\overline{\lambda })$ is a continuous monotonically decreasing function with respect to $\overline{\lambda }$ and satisfies $\underset{Re(\overline{\lambda })\to +\mathrm{\infty }}{lim}ϝ(\overline{\lambda })=0$. Thus Equation (25(25) $1=2\kappa (1-\vartheta ){\int }_0^{a^{+}}\overline{b}(a)e^{-\overline{\lambda }a-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a.$(25) ) has and only has a unique real root ${\overline{\lambda }}^{\ast }$. Moreover, for the special case $\overline{\lambda }=0$ we have $ϝ(0)=2\kappa (1-\vartheta ){\int }_0^{a^{+}}\overline{b}(a)e^{-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a.$If condition (26) $2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a)e^{-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a<1$(26) is satisfied, we can obtain from the monotonicity of $ϝ(\overline{\lambda })$ that all the root real parts of $ϝ(\overline{\lambda })$ are less than 0 when Equation (25(25) $1=2\kappa (1-\vartheta ){\int }_0^{a^{+}}\overline{b}(a)e^{-\overline{\lambda }a-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a.$(25) ) holds.

Thus for the tumour-free steady state, the following theorem can be obtained.

Theorem 4.1

Let (H1) and (H3) hold. Moreover, if Equation (26(26) $2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a)e^{-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a<1$(26) ) is satisfied, then the tumour-free steady state ${\overline{E}}_0:=(0,0,0)$ is locally asymptotically stable.

4.2. Global stability of tumour-free steady state

Here, we discuss the global stability of the tumour-free stable state ${\overline{E}}_0$ under two different assumptions. We first make the following assumption.

(H4) Assume $b(a,N(t))=b(a)\mathrm{\Psi }(N(t))$ with $\mathrm{\Psi }(0)=1,\underset{x\to +\mathrm{\infty }}{lim}\mathrm{\Psi }(x)=0$. Moreover, $\mathrm{\Psi }(x)$ is continuously differentiable and $\frac{\mathrm{d}\mathrm{\Psi }(x)}{\mathrm{d}x}<0$ in $x\in [0,+\mathrm{\infty })$.

Remark 1

The nonlinear proliferation function $b(a,N(t))$ is related to the age a and the total number of tumour cells N. Assume that $b(a,N(t))$ has a separable form $b(a)\mathrm{\Psi }(N(t))$, where $b(a)$ is the proliferation function of tumour cells with respect to age a when there is no restriction from the number of tumour cells N. $\mathrm{\Psi }(N(t))$ expresses the restriction on proliferation by the total number of tumours N, and the larger the N, the smaller the proliferation ability. Hence, we assume that $\mathrm{\Psi }(0)=1,\underset{x\to +\mathrm{\infty }}{lim}\mathrm{\Psi }(x)=0$.

If the assumption (H4) is satisfied, the following theorem can be obtained.

Theorem 4.2

If assumptions (H1) and (H4) hold and (27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) is satisfied, the tumour-free stable state ${\overline{E}}_0$ is globally asymptotically stable.

Proof.

From system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ), we have (28) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)p(t,a)-{\mu }_p(a)P(t,a),\\ P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ P(0,a)=P_0(a).\end{array}$(28) Rewrite System (28(28) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)p(t,a)-{\mu }_p(a)P(t,a),\\ P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ P(0,a)=P_0(a).\end{array}$(28) ) as a Cauchy problem for the Banach space $X:=\mathbb{R}\times L^1(0,a^{+})$ as follows. (29) $\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}t}(\genfrac{}{}{0ex}{}{0}{u(t,\cdot )})=B(\genfrac{}{}{0ex}{}{0}{u(t,\cdot )})+(\genfrac{}{}{0ex}{}{2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a}{-b(a,N(t))P(t,\cdot )}),\\ (\genfrac{}{}{0ex}{}{0}{u(0,\cdot )})=(\genfrac{}{}{0ex}{}{0}{P(0,\cdot )})\in X,\end{array}$(29) where the linear operator $B:D(B)\subset X\to X$ is defined as $B(\genfrac{}{}{0ex}{}{0}{\varphi })=(\genfrac{}{}{0ex}{}{-\varphi (0)}{-\kappa {\varphi }^{\prime }-{\mu }_p(\cdot )\varphi -\vartheta (\cdot )\varphi })$and $D(B):=\{0\}\times AC[0,a^{+}]$.

Construct linear comparison equations (30) $\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}t}(\genfrac{}{}{0ex}{}{0}{z(t,\cdot )})=B(\genfrac{}{}{0ex}{}{0}{z(t,\cdot )})+(\genfrac{}{}{0ex}{}{2(1-\ell ){\int }_0^{a^{+}}b(a)z(t,a)\mathrm{d}a}{0}),\\ (\genfrac{}{}{0ex}{}{0}{z(0,\cdot )})=(\genfrac{}{}{0ex}{}{0}{P(0,\cdot )})\in X.\end{array}$(30) Let $\hat{u}(t,\cdot )$ and $\hat{z}(t,\cdot )$ be a integral solution of systems (29(29) $\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}t}(\genfrac{}{}{0ex}{}{0}{u(t,\cdot )})=B(\genfrac{}{}{0ex}{}{0}{u(t,\cdot )})+(\genfrac{}{}{0ex}{}{2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a}{-b(a,N(t))P(t,\cdot )}),\\ (\genfrac{}{}{0ex}{}{0}{u(0,\cdot )})=(\genfrac{}{}{0ex}{}{0}{P(0,\cdot )})\in X,\end{array}$(29) ) and (30(30) $\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}t}(\genfrac{}{}{0ex}{}{0}{z(t,\cdot )})=B(\genfrac{}{}{0ex}{}{0}{z(t,\cdot )})+(\genfrac{}{}{0ex}{}{2(1-\ell ){\int }_0^{a^{+}}b(a)z(t,a)\mathrm{d}a}{0}),\\ (\genfrac{}{}{0ex}{}{0}{z(0,\cdot )})=(\genfrac{}{}{0ex}{}{0}{P(0,\cdot )})\in X.\end{array}$(30) ), respectively. When assumption (H4) holds, applying the comparability principle we get (31) $0\le \hat{u}(t,\cdot )\le \hat{z}(t,\cdot )$(31) for all $t\ge 0$. For linear problem (30(30) $\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}t}(\genfrac{}{}{0ex}{}{0}{z(t,\cdot )})=B(\genfrac{}{}{0ex}{}{0}{z(t,\cdot )})+(\genfrac{}{}{0ex}{}{2(1-\ell ){\int }_0^{a^{+}}b(a)z(t,a)\mathrm{d}a}{0}),\\ (\genfrac{}{}{0ex}{}{0}{z(0,\cdot )})=(\genfrac{}{}{0ex}{}{0}{P(0,\cdot )})\in X.\end{array}$(30) ), assume that (32) $z(t,a)=e^{{}^{\lambda t}}\overline{z}(a).$(32) The linear abstract problem (30(30) $\{\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}t}(\genfrac{}{}{0ex}{}{0}{z(t,\cdot )})=B(\genfrac{}{}{0ex}{}{0}{z(t,\cdot )})+(\genfrac{}{}{0ex}{}{2(1-\ell ){\int }_0^{a^{+}}b(a)z(t,a)\mathrm{d}a}{0}),\\ (\genfrac{}{}{0ex}{}{0}{z(0,\cdot )})=(\genfrac{}{}{0ex}{}{0}{P(0,\cdot )})\in X.\end{array}$(30) ) becomes (33) $\{\begin{array}{l}\frac{\mathrm{d}\overline{z}(a)}{\mathrm{d}a}=-\overline{\lambda }\overline{z}(a)-{\overline{\mu }}_p(a)\overline{z}(a)-\overline{\vartheta }(a)\overline{z}(a),\\ \overline{z}(0)=2(1-\ell ){\int }_0^{a^{+}}b(a)\overline{z}(a)\mathrm{d}a,\end{array}$(33) where $\overline{\lambda }=\frac{\lambda }{\kappa }$, ${\overline{\mu }}_p(a)$ and $\overline{\vartheta }(a)$ are defined in (19(19) $\overline{b}(a,N)=\frac{b(a,N)}{\kappa },\overline{\vartheta }(a)=\frac{\vartheta (a)}{\kappa },{\overline{\mu }}_p(a)=\frac{{\mu }_p(a)}{\kappa },\overline{d}=\frac{\mathrm{d}}{\kappa },{\overline{\mu }}_q=\frac{{\mu }_q}{\kappa }$(19) ). Solving the first equation of Equation (33(33) $\{\begin{array}{l}\frac{\mathrm{d}\overline{z}(a)}{\mathrm{d}a}=-\overline{\lambda }\overline{z}(a)-{\overline{\mu }}_p(a)\overline{z}(a)-\overline{\vartheta }(a)\overline{z}(a),\\ \overline{z}(0)=2(1-\ell ){\int }_0^{a^{+}}b(a)\overline{z}(a)\mathrm{d}a,\end{array}$(33) ) and substituting it to the second equation yields (34) $\overline{z}(0)=2(1-\ell ){\int }_0^{a^{+}}b(a)\overline{z}(0)e^{-{\int }_0^a\overline{\lambda }+{\overline{\mu }}_p(\xi )+\overline{\vartheta }(\xi )\mathrm{d}\xi }\mathrm{d}a,$(34) where $\overline{z}(0)\ne 0$ ($\overline{z}(0)=0$ leads to the tumour-free steady state). From Equation (34(34) $\overline{z}(0)=2(1-\ell ){\int }_0^{a^{+}}b(a)\overline{z}(0)e^{-{\int }_0^a\overline{\lambda }+{\overline{\mu }}_p(\xi )+\overline{\vartheta }(\xi )\mathrm{d}\xi }\mathrm{d}a,$(34) ) we get the following characteristic equation: (35) $1=2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\overline{\lambda }+{\overline{\mu }}_p(\xi )+\overline{\vartheta }(\xi )\mathrm{d}\xi }\mathrm{d}a.$(35) If condition (27(27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) ) is satisfied, then all roots of Equation (35(35) $1=2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\overline{\lambda }+{\overline{\mu }}_p(\xi )+\overline{\vartheta }(\xi )\mathrm{d}\xi }\mathrm{d}a.$(35) ) only have negative real parts.

Therefore, we get (36) $0\le \underset{t\to \mathrm{\infty }}{lim}supP(t,a)=\underset{t\to \mathrm{\infty }}{lim}sup\hat{u}(t,a)\le \underset{t\to \mathrm{\infty }}{lim}sup\hat{z}(t,a)=\underset{t\to \mathrm{\infty }}{lim}{e}^{\overline{\lambda }t}\overline{z}(a)=0.$(36) Then, for t large enough, we get $\frac{\mathrm{d}Q}{\mathrm{d}t}=-{\mu }_{q}Q(t),Q(0)=Q_0\ge 0⟹Q(t)=Q_0{e}^{-{\mu }_{q}t}.$And hence $0\le \underset{t\to +\mathrm{\infty }}{lim}supQ(t)\le \underset{t\to +\mathrm{\infty }}{lim}sup{Q}_0{e}^{-{\mu }_{q}t}=0.$In the same way, we get $0\le \underset{t\to +\mathrm{\infty }}{lim}supD(t)\le \underset{t\to +\mathrm{\infty }}{lim}sup{D}_0{e}^{-\mathrm{d}t}=0.$Let $N(t,a)=P(t,a)+Q(t,a)+D(t,a)$, then we obtain $0\le \underset{t\to \mathrm{\infty }}{lim}supP(t,a)+\underset{t\to \mathrm{\infty }}{lim}supQ(t,a)+\underset{t\to \mathrm{\infty }}{lim}supD(t,a)\le \underset{t\to \mathrm{\infty }}{lim}supN(t,a)=0.$End of this proof.

When we choose assumption (H2), the global stability of the tumour-free steady state is also obtained under specific conditions.

Theorem 4.3

When assumptions (H1) and (H2) are true, and condition (37) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{1}{k}{\int }_0^a{\gamma }_{0}b(\xi )+\sigma (\xi )+{\mu }_p(\xi )\mathrm{d}\xi }\mathrm{d}a<1,$(37) is satisfied, the tumour-free stable state ${\overline{E}}_0$ is globally asymptotically stable, where ${\gamma }_0:=\psi (M)$ and M is the ultimate upper bound of the system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ).

Proof.

Transform the system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) into the Cauchy problem in Banach space $X:=\mathbb{R}\times \mathbb{R}\times \mathbb{R}\times L^1(0,a^{+})\times L^1(0,a^{+})\times L^1(0,a^{+})$ as (38) $\{\begin{array}{ll}& \frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}0\\ 0\\ 0\\ u_1(t,\cdot )\\ u_2(t,\cdot )\\ u_3(t,\cdot )\end{array}\right)=H\left(\begin{array}{c}0\\ 0\\ 0\\ u_1(t,\cdot )\\ u_2(t,\cdot )\\ u_3(t,\cdot )\end{array}\right)\\ & +\left(\begin{array}{c}2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))u_1(t,a)\mathrm{d}a\\ \frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))u_1(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)u_1(t,a)\mathrm{d}a]\\ \frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}{u}_2(t,a)\mathrm{d}a\\ -[b(a,N(t))u_1(t,\cdot )+\vartheta (a)u_1(t,\cdot )+{\mu }_p(a)u_1(t,\cdot )]\\ -{\mu }_p{u}_2(t,\cdot )\\ -\mathrm{d}{u}_3(t,\cdot )\end{array}\right),\\ & (0,0,0,u_1(0,\cdot ),u_2(0,\cdot ),u_3(0,\cdot ){)}^T=(0,0,0,P(0,\cdot ),Q(0,\cdot ),D(0,\cdot ){)}^T\in X.\end{array}$(38) Define the linear operator H from X to X as (39) $H(0,0,0,{\varphi }_1,{\varphi }_2,{\varphi }_3{)}^T:=(-{\varphi }_1(0),-{\varphi }_2(0),-{\varphi }_3(0),-k{{\varphi }_1}^{\prime },-k{{\varphi }_2}^{\prime },-k{{\varphi }_3}^{\prime }{)}^T$(39) and $D(H)$ is given as $D(H):=\{0\}\times \{0\}\times \{0\}\times AC[0,a^{+}]\times AC[0,a^{+}]\times AC[0,a^{+}]$. Similar to the approach in Theorem 4.2, we construct the linear abstract equations (40) $\{\begin{array}{ll}& \frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}0\\ 0\\ 0\\ {\mathit{z}}_1(t,\cdot )\\ {\mathit{z}}_2(t,\cdot )\\ {\mathit{z}}_3(t,\cdot )\end{array}\right)=H\left(\begin{array}{c}0\\ 0\\ 0\\ {\mathit{z}}_1(t,\cdot )\\ {\mathit{z}}_2(t,\cdot )\\ {\mathit{z}}_3(t,\cdot )\end{array}\right)\\ & +\left(\begin{array}{c}2(1-\ell ){\int }_0^{a^{+}}b(a){\mathit{z}}_1(t,a)\mathrm{d}a\\ \frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a){\mathit{z}}_1(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a){\mathit{z}}_1(t,a)\mathrm{d}a]\\ \frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}{\mathit{z}}_2(t,a)\mathrm{d}a\\ -[b(a){\gamma }_0{\mathit{z}}_1(t,\cdot )+\vartheta (a){\mathit{z}}_1(t,\cdot )+{\mu }_p(a){\mathit{z}}_1(t,\cdot )]\\ -{\mu }_p{\mathit{z}}_2(t,\cdot )\\ -\mathrm{d}{\mathit{z}}_3(t,\cdot )\end{array}\right),\\ & (0,0,0,{\mathit{z}}_1(0,\cdot ),{\mathit{z}}_2(0,\cdot ),{\mathit{z}}_3(0,\cdot ){)}^T=(0,0,0,P(0,\cdot ),Q(0,\cdot ),D(0,\cdot ){)}^T\in X.\end{array}$(40) Let the integral solutions of system (38(38) $\{\begin{array}{ll}& \frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}0\\ 0\\ 0\\ u_1(t,\cdot )\\ u_2(t,\cdot )\\ u_3(t,\cdot )\end{array}\right)=H\left(\begin{array}{c}0\\ 0\\ 0\\ u_1(t,\cdot )\\ u_2(t,\cdot )\\ u_3(t,\cdot )\end{array}\right)\\ & +\left(\begin{array}{c}2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))u_1(t,a)\mathrm{d}a\\ \frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))u_1(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)u_1(t,a)\mathrm{d}a]\\ \frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}{u}_2(t,a)\mathrm{d}a\\ -[b(a,N(t))u_1(t,\cdot )+\vartheta (a)u_1(t,\cdot )+{\mu }_p(a)u_1(t,\cdot )]\\ -{\mu }_p{u}_2(t,\cdot )\\ -\mathrm{d}{u}_3(t,\cdot )\end{array}\right),\\ & (0,0,0,u_1(0,\cdot ),u_2(0,\cdot ),u_3(0,\cdot ){)}^T=(0,0,0,P(0,\cdot ),Q(0,\cdot ),D(0,\cdot ){)}^T\in X.\end{array}$(38) ) and system (41(41) ${\hat{u}}_1(t,\cdot )\le {\hat{\mathit{z}}}_1(t,\cdot ),{\hat{u}}_2(t,\cdot )\le {\hat{\mathit{z}}}_2(t,\cdot ),{\hat{u}}_3(t,\cdot )\le {\hat{\mathit{z}}}_3(t,\cdot ),t\ge 0.$(41) ) be $({\hat{u}}_1(t,\cdot ),{\hat{u}}_2(t,\cdot ),{\hat{u}}_3(t,\cdot ))$ and $({\hat{\mathit{z}}}_1(t,\cdot ),{\hat{\mathit{z}}}_2(t,\cdot ),{\hat{\mathit{z}}}_3(t,\cdot ))$, respectively, and using the principle of comparability, we can know that (41) ${\hat{u}}_1(t,\cdot )\le {\hat{\mathit{z}}}_1(t,\cdot ),{\hat{u}}_2(t,\cdot )\le {\hat{\mathit{z}}}_2(t,\cdot ),{\hat{u}}_3(t,\cdot )\le {\hat{\mathit{z}}}_3(t,\cdot ),t\ge 0.$(41) Define the operator C from $D(H)$ to X as $C\left(\begin{array}{c}0\\ 0\\ 0\\ {\varphi }_1\\ {\varphi }_2\\ {\varphi }_3\end{array}\right)=\left(\begin{array}{c}2(1-\ell ){\int }_0^{a^{+}}b(a){\varphi }_1(a)\mathrm{d}a\\ \frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a){\varphi }_1(a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a){\varphi }_1(a)\mathrm{d}a]\\ \frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}{\varphi }_2(a)\mathrm{d}a\\ -[b(\cdot ){\gamma }_0{\varphi }_1+\vartheta (\cdot ){\varphi }_1+{\mu }_p(\cdot ){\varphi }_1]\\ -{\mu }_p{\varphi }_2\\ -\mathrm{d}{\varphi }_3\end{array}\right).$Let the eigenvector of the operator H + C be defined as ω. The characteristic equation satisfied by system (41(41) ${\hat{u}}_1(t,\cdot )\le {\hat{\mathit{z}}}_1(t,\cdot ),{\hat{u}}_2(t,\cdot )\le {\hat{\mathit{z}}}_2(t,\cdot ),{\hat{u}}_3(t,\cdot )\le {\hat{\mathit{z}}}_3(t,\cdot ),t\ge 0.$(41) ) can be obtained as follows. (42) $1=2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{\omega a}{\kappa }-\frac{1}{\kappa }{\int }_0^a{\gamma }_{0}b(\xi )+\vartheta (\xi )+{\mu }_p(\xi )\mathrm{d}\xi }\mathrm{d}a.$(42) If condition (37(37) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{1}{k}{\int }_0^a{\gamma }_{0}b(\xi )+\sigma (\xi )+{\mu }_p(\xi )\mathrm{d}\xi }\mathrm{d}a<1,$(37) ) is satisfied, we can obtain that all eigenvectors ω have only negative real parts. Hence, we get $\begin{array}{rl}& \underset{t\to +\mathrm{\infty }}{lim}sup\Vert {\hat{u}}_1(t,\cdot )\Vert \le \underset{t\to +\mathrm{\infty }}{lim}sup\Vert {\hat{\mathit{z}}}_1(t,\cdot )\Vert =\underset{t\to +\mathrm{\infty }}{lim}\mid P(0)e^{\frac{\omega t}{k}}\mid =0,\\ & \underset{t\to +\mathrm{\infty }}{lim}sup\Vert {\hat{u}}_2(t,\cdot )\Vert \le \underset{t\to +\mathrm{\infty }}{lim}sup\Vert {\hat{\mathit{z}}}_2(t,\cdot )\Vert =\underset{t\to +\mathrm{\infty }}{lim}\mid Q(0)e^{\frac{\omega t}{k}}\mid =0,\\ & \underset{t\to +\mathrm{\infty }}{lim}sup\Vert {\hat{u}}_3(t,\cdot )\Vert \le \underset{t\to +\mathrm{\infty }}{lim}sup\Vert {\hat{\mathit{z}}}_3(t,\cdot )\Vert =\underset{t\to +\mathrm{\infty }}{lim}\mid D(0)e^{\frac{\omega t}{k}}\mid =0,\end{array}$where $P(0):={\int }_0^{a^{+}}P(0,a)\mathrm{d}a,Q(0):={\int }_0^{a^{+}}Q(0,a)\mathrm{d}a,D(0):={\int }_0^{a^{+}}D(0,a)\mathrm{d}a$. That is $\begin{array}{rl}& \underset{t\to +\mathrm{\infty }}{lim}sup{\int }_0^{a^{+}}P(t,a)\mathrm{d}a=0,\\ & \underset{t\to +\mathrm{\infty }}{lim}sup{\int }_0^{a^{+}}Q(t,a)\mathrm{d}a=0,\underset{t\to +\mathrm{\infty }}{lim}sup{\int }_0^{a^{+}}D(t,a)\mathrm{d}a=0.\end{array}$Since $P(t,a)$, $Q(t,a)$ and $D(t,a)$ are all nonnegative, we have $\underset{t\to +\mathrm{\infty }}{lim}supP(t,a)=0,\underset{t\to +\mathrm{\infty }}{lim}supQ(t,a)=0,\underset{t\to +\mathrm{\infty }}{lim}supD(t,a)=0.$This proof is completed.

We proved the existence and local stability of the tumour-free steady state. Conditions for the global stability of the tumour-free steady state were obtained under different assumptions. From the conditions that need to be satisfied in Theorems 4.2 and 4.3, we know that condition (27(27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) ) is stronger, and when Equation (27(27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) ) holds, Equations (20(20) $\{\begin{array}{l}\frac{\mathrm{\partial }x}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }x}{\mathrm{\partial }a}=-(b(a)+\vartheta (a)+{\mu }_p(a))x(t,a),\\ \frac{\mathrm{\partial }y}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }y}{\mathrm{\partial }a}=-{\mu }_{q}y(t,a),\\ \frac{\mathrm{\partial }z}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }z}{\mathrm{\partial }a}=-\mathrm{d}z(t,a)\end{array}$(20) ) and  (37(37) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{1}{k}{\int }_0^a{\gamma }_{0}b(\xi )+\sigma (\xi )+{\mu }_p(\xi )\mathrm{d}\xi }\mathrm{d}a<1,$(37) ) also hold. Theorems 4.2 and 4.3 prove the global asymptotic stability of the tumour-free stable state under different assumptions. Assumption (H2) has more requirements relative to assumption (H4) that is, assumption (H2) is more strongly conditioned, and assumption (H4) is necessarily satisfied when assumption (H2) is satisfied. The interval between the two theorems lies in the difference in the assumptions and the strength of the threshold conditions.

5. Existence and local stability of the tumour steady state

In this section, we first discuss the existence of the positive steady state of system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) and then study the local stability of it.

5.1. Existence of the tumour steady state

Define the tumour stable state of system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) as ${\overline{E}}^{\ast }(a)=({\overline{P}}^{\ast }(a),{\overline{Q}}^{\ast }(a),{\overline{D}}^{\ast }(a))$. ${\overline{E}}^{\ast }(a)$ should satisfy the time-independent system (17(17) $\{\begin{array}{l}\frac{\mathrm{d}\overline{P}}{\mathrm{d}a}=-\overline{b}(a,N)\overline{P}(a)-\overline{\vartheta }(a)\overline{P}(a)-{\overline{\mu }}_p(a)\overline{P}(a),\\ \frac{\mathrm{d}\overline{Q}}{\mathrm{d}a}=-{\overline{\mu }}_q\overline{Q}(a),\\ \frac{\mathrm{d}\overline{D}}{\mathrm{d}a}=-\overline{d}\overline{D}(a)\end{array}$(17) )–( 18(18) $\begin{array}{rl}\overline{P}(0)& =2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a,\\ \overline{Q}(0)& =2\ell {\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{P}(a)\mathrm{d}a,\\ \overline{D}(0)& ={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{Q}(a)\mathrm{d}a,\end{array}$(18) ). A simple calculation of the time-independent system (17(17) $\{\begin{array}{l}\frac{\mathrm{d}\overline{P}}{\mathrm{d}a}=-\overline{b}(a,N)\overline{P}(a)-\overline{\vartheta }(a)\overline{P}(a)-{\overline{\mu }}_p(a)\overline{P}(a),\\ \frac{\mathrm{d}\overline{Q}}{\mathrm{d}a}=-{\overline{\mu }}_q\overline{Q}(a),\\ \frac{\mathrm{d}\overline{D}}{\mathrm{d}a}=-\overline{d}\overline{D}(a)\end{array}$(17) )–( 18(18) $\begin{array}{rl}\overline{P}(0)& =2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a,\\ \overline{Q}(0)& =2\ell {\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{P}(a)\mathrm{d}a,\\ \overline{D}(0)& ={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{Q}(a)\mathrm{d}a,\end{array}$(18) ) yields (43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) where $\overline{\epsilon }(a,N)=\overline{b}(a,N)+\overline{\sigma }(a)+{\overline{\mu }}_p(a)$. Substituting the second equation of (18(18) $\begin{array}{rl}\overline{P}(0)& =2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a,\\ \overline{Q}(0)& =2\ell {\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{P}(a)\mathrm{d}a,\\ \overline{D}(0)& ={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{Q}(a)\mathrm{d}a,\end{array}$(18) ) and Equation (43(43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) ) into the integral solution of the second equation of (17(17) $\{\begin{array}{l}\frac{\mathrm{d}\overline{P}}{\mathrm{d}a}=-\overline{b}(a,N)\overline{P}(a)-\overline{\vartheta }(a)\overline{P}(a)-{\overline{\mu }}_p(a)\overline{P}(a),\\ \frac{\mathrm{d}\overline{Q}}{\mathrm{d}a}=-{\overline{\mu }}_q\overline{Q}(a),\\ \frac{\mathrm{d}\overline{D}}{\mathrm{d}a}=-\overline{d}\overline{D}(a)\end{array}$(17) ) yields (44) $\begin{array}{rl}{\overline{Q}}^{\ast }(a)& ={\overline{Q}}^{\ast }(0)e^{-{\overline{\mu }}_{q}a}\\ & =[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N){\overline{P}}^{\ast }(\xi )\mathrm{d}\xi +{\int }_0^{a^{+}}\overline{\vartheta }(\xi ){\overline{P}}^{\ast }(\xi )\mathrm{d}\xi ]e^{-{\overline{\mu }}_{q}a}\\ & ={\overline{P}}^{\ast }(0)[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_{q}a}.\end{array}$(44) In the same way, we get (45) $\begin{array}{rl}{\overline{D}}^{\ast }(a)& ={\overline{D}}^{\ast }(0)e^{-\overline{d}a}\\ & ={\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}{\overline{Q}}^{\ast }(\varpi )\mathrm{d}\varpi \\ & ={\overline{P}}^{\ast }(0){\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi .\end{array}$(45) From Equations (43(43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) )–( 45(45) $\begin{array}{rl}{\overline{D}}^{\ast }(a)& ={\overline{D}}^{\ast }(0)e^{-\overline{d}a}\\ & ={\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}{\overline{Q}}^{\ast }(\varpi )\mathrm{d}\varpi \\ & ={\overline{P}}^{\ast }(0){\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi .\end{array}$(45) ) we can know that the existence of positive solutions of ${\overline{Q}}^{\ast }(a)$, ${\overline{D}}^{\ast }(a)$ and ${\overline{P}}^{\ast }(a)$ depend on ${\overline{P}}^{\ast }(0)$. Thus the existence of tumour steady state for system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) is consistent with the existence of the positive solution of Equation (43(43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) ). By using Equation (43(43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) ) and the boundary condition ${\overline{P}}^{\ast }(0)$ can be represented as (46) ${\overline{P}}^{\ast }(0)={\overline{P}}^{\ast }(0)U(N),$(46) where (47) $U(N)=2(1-\vartheta ){\int }_0^{a^{+}}b(a,N)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi }\mathrm{d}a.$(47) Equation (46(46) ${\overline{P}}^{\ast }(0)={\overline{P}}^{\ast }(0)U(N),$(46) ) has a unique non-zero solution when $U(N)=1$. The biological significance of $U(N)=1$ is that at tumour steady state a tumour cell produces an average of one tumour cell during its entire life cycle.

Theorem 5.1

Let assumption (H1) and N>0 hold, then (48) $U(N)=1,$(48) is a sufficient and necessary condition for the existence of the tumour steady state. When in this case, the tumour steady state ${\overline{E}}^{\ast }(a)$ can be represented by N in the form of (43(43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) )–( 45(45) $\begin{array}{rl}{\overline{D}}^{\ast }(a)& ={\overline{D}}^{\ast }(0)e^{-\overline{d}a}\\ & ={\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}{\overline{Q}}^{\ast }(\varpi )\mathrm{d}\varpi \\ & ={\overline{P}}^{\ast }(0){\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi .\end{array}$(45) ) with (49) $\begin{array}{rl}{\overline{P}}^{\ast }(0)& =N({\int }_0^{a^{+}}(e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi }+{\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_{q}a}\\ & +{\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & {+{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi )\mathrm{d}a)}^{-1}.\end{array}$(49)

Proof.

Combining Equations (43(43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) )–( 45(45) $\begin{array}{rl}{\overline{D}}^{\ast }(a)& ={\overline{D}}^{\ast }(0)e^{-\overline{d}a}\\ & ={\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}{\overline{Q}}^{\ast }(\varpi )\mathrm{d}\varpi \\ & ={\overline{P}}^{\ast }(0){\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi .\end{array}$(45) ) we can deduce that (50) $\begin{array}{rl}N& ={\int }_0^{a^{+}}({\overline{P}}^{\ast }(a)+{\overline{Q}}^{\ast }(a)+{\overline{D}}^{\ast }(a))\mathrm{d}a\\ & ={\int }_0^{a^{+}}({\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi }+{\overline{P}}^{\ast }(0){\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_{q}a}\\ & +{\overline{P}}^{\ast }(0){\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi )\mathrm{d}a\\ & ={\overline{P}}^{\ast }(0){\int }_0^{a^{+}}(e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi }+{\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_{q}a}\\ & +{\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi )\mathrm{d}a.\end{array}$(50) Then, it can be seen that Equation (49(49) $\begin{array}{rl}{\overline{P}}^{\ast }(0)& =N({\int }_0^{a^{+}}(e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi }+{\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_{q}a}\\ & +{\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & {+{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi )\mathrm{d}a)}^{-1}.\end{array}$(49) ) holds. In order to complete the proof, it can be sufficiently showed that the combination (43(43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) )–( 45(45) $\begin{array}{rl}{\overline{D}}^{\ast }(a)& ={\overline{D}}^{\ast }(0)e^{-\overline{d}a}\\ & ={\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}{\overline{Q}}^{\ast }(\varpi )\mathrm{d}\varpi \\ & ={\overline{P}}^{\ast }(0){\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi .\end{array}$(45) ) and (48(48) $U(N)=1,$(48) ) is equivalent to Equation (46(46) ${\overline{P}}^{\ast }(0)={\overline{P}}^{\ast }(0)U(N),$(46) ). The analysis in this section just examines it is right. This proof is completed.

5.2. Local stability of the tumour steady state

To study the tumour steady state ${\overline{E}}^{\ast }(a)$ of the system we enhanced the assumption (H3) as follows.

(H5) $b(a,N)$ is differentiable relative to N, and the partial derivative $\frac{\mathrm{\partial }b(a,N)}{\mathrm{\partial }N}$ of $b(a,N)$ with respect to N is bounded in the interval $[0,a^{+}]$.

The tumour steady state ${\overline{E}}^{\ast }(a)$ local perturbation form is defined as $(x(t,a),y(t,a),z(t,a))$, i.e. $x(t,a)=P(t,a)-{\overline{P}}^{\ast }(a),y(t,a)=Q(t,a)-{\overline{Q}}^{\ast }(a),z(t,a)=D(t,a)-{\overline{D}}^{\ast }(a).$Let $n(t)=N(t)-N^{\ast }$, where $N^{\ast }={\int }_0^{a^{+}}({\overline{P}}^{\ast }(a)+{\overline{Q}}^{\ast }(a)+{\overline{D}}^{\ast }(a))\mathrm{d}a$. It can be seen that $N^{\ast }$ is a constant and depend on the steady state. Then, $n(t)={\int }_0^{a^{+}}(x(t,a)+y(t,a)+z(t,a))\mathrm{d}a.$Based on basic system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ), the perturbed form of the tumour steady state needs to satisfy (51) $\{\begin{array}{l}\frac{\mathrm{\partial }x}{\mathrm{\partial }t}+k\frac{\mathrm{\partial }x}{\mathrm{\partial }a}=-b(a,N^{\ast })x(t,a)-[\vartheta (a)+{\mu }_p(a)]x(t,a)\\ -b_N(a,N^{\ast }){\overline{P}}^{\ast }(a)n(t)-{\mathrm{\Gamma }}_0(t,a),\\ \frac{\mathrm{\partial }y}{\mathrm{\partial }t}+k\frac{\mathrm{\partial }y}{\mathrm{\partial }a}=-{\mu }_{q}y(t,a),\\ \frac{\mathrm{\partial }z}{\mathrm{\partial }t}+k\frac{\mathrm{\partial }z}{\mathrm{\partial }a}=-\mathrm{d}z(t,a),\\ x(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N^{\ast })x(t,a)\mathrm{d}a+2(1-\ell )\delta (N^{\ast })n(t)\\ +2(1-\ell ){\int }_0^{a^{+}}{\mathrm{\Gamma }}_0(t,a)\mathrm{d}a,\\ y(t,0)=\frac{2\ell }{\kappa }{\int }_0^{a^{+}}b(a,N^{\ast })x(t,a)\mathrm{d}a+\frac{1}{\kappa }{\int }_0^{a^{+}}\vartheta (a)x(t,a)\mathrm{d}a\\ +\frac{2\ell }{\kappa }\delta (N^{\ast })n(t)+\frac{2\ell }{\kappa }{\int }_0^{a^{+}}{\mathrm{\Gamma }}_0(t,a)\mathrm{d}a,\\ z(t,0)=\frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}y(t,a)\mathrm{d}a,\end{array}$(51) where $\begin{array}{rl}{\mathrm{\Gamma }}_0(t,a)& =b_N(a,N^{\ast })n(t)x(t,a)+\mathrm{\Lambda }(a,n)({\overline{P}}^{\ast }(a)+x(t,a)),\\ \delta (N^{\ast })& ={\int }_0^{a^{+}}{b}_N(a,N^{\ast }){\overline{P}}^{\ast }(a)\mathrm{d}a,\end{array}$and $\mathrm{\Lambda }(a,n)=b(a,N^{\ast }+n(t))-b(a,N^{\ast })-b_N(a,N^{\ast })n(t).$The local properties are determined by the lower order terms, and we neglect the higher order terms in system (48(48) $U(N)=1,$(48) ) to obtain the following linear system. (52) $\{\begin{array}{l}\frac{\mathrm{\partial }x}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }x}{\mathrm{\partial }a}=-b(a,N^{\ast })x(t,a)-[\vartheta (a)+{\mu }_p(a)]x(t,a)-b_N(a,N^{\ast }){\overline{P}}^{\ast }(a)n(t),\\ \frac{\mathrm{\partial }y}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }y}{\mathrm{\partial }a}=-{\mu }_{q}y(t,a),\\ \frac{\mathrm{\partial }z}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }z}{\mathrm{\partial }a}=-\mathrm{d}z(t,a),\\ x(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N^{\ast })x(t,a)\mathrm{d}a+2(1-\ell )\delta (N^{\ast })n(t),\\ y(t,0)=\frac{2\ell }{\kappa }{\int }_0^{a^{+}}b(a,N^{\ast })x(t,a)\mathrm{d}a+\frac{1}{\kappa }{\int }_0^{a^{+}}\vartheta (a)x(t,a)\mathrm{d}a+\frac{2\ell }{\kappa }\delta (N^{\ast })n(t),\\ z(t,0)=\frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}y(t,a)\mathrm{d}a.\end{array}$(52) Based on assumption (23(23) $\{\begin{array}{l}\frac{\mathrm{d}x}{\mathrm{d}t}=-\overline{\lambda }\overline{x}(a)-(\overline{b}(a)+\overline{\vartheta }(a)+{\overline{\mu }}_p(a))\overline{x}(a),\\ \frac{\mathrm{d}y}{\mathrm{d}t}=-\overline{\lambda }\overline{y}(a)-{\overline{\mu }}_q\overline{y}(a),\\ \frac{\mathrm{d}z}{\mathrm{d}t}=-\overline{\lambda }\overline{z}(a)-\overline{d}\overline{z}(a),\\ \overline{x}(0)=2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a)\overline{x}(a)\mathrm{d}a,\\ \overline{y}(0)=2\ell {\int }_0^{a^{+}}\overline{b}(a)\overline{x}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{x}(a)\mathrm{d}a,\\ \overline{z}(0)={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{y}(a)\mathrm{d}a.\end{array}$(23) ), we can obtain (53) $\{\begin{array}{l}\frac{\overline{x}(a)}{\mathrm{d}a}=-\overline{\lambda }\overline{x}(a)-\overline{b}(a,N^{\ast })\overline{x}(a)-[\overline{\vartheta }(a)+{\overline{\mu }}_p(a)]\overline{x}(a)-{\overline{b}}_N(a,N^{\ast }){\overline{P}}^{\ast }(a)n^{\ast },\\ \frac{\overline{y}(a)}{\mathrm{d}a}=-\overline{\lambda }\overline{y}(a)-{\overline{\mu }}_q\overline{y}(a),\\ \frac{\overline{z}(a)}{\mathrm{d}a}=-\overline{\lambda }\overline{z}(a)-\overline{d}\overline{z}(a),\\ \overline{x}(0)=2(1-\ell )\kappa {\int }_0^{a^{+}}\overline{b}(a,N^{\ast })\overline{x}(a)\mathrm{d}a+2(1-\ell )\kappa \overline{\delta }(N^{\ast })n^{\ast },\\ \overline{y}(0)=2\ell {\int }_0^{a^{+}}\overline{b}(a,N^{\ast })\overline{x}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{x}(a)\mathrm{d}a+2\ell \overline{\delta }(N^{\ast })n^{\ast },\\ \overline{z}(0)=\frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}\overline{y}(a)\mathrm{d}a,\end{array}$(53) where $\begin{array}{rl}\overline{\lambda }& =\frac{\lambda }{\kappa },\overline{b}(a,N^{\ast })=\frac{b(a,N^{\ast })}{\kappa },{\overline{b}}_N(a,N^{\ast })=\frac{b_N(a,N^{\ast })}{\kappa },\\ n^{\ast }& ={\int }_0^{a^{+}}\overline{x}(a)+\overline{y}(a)+\overline{z}(a)\mathrm{d}a,\overline{\delta }(N^{\ast })={\int }_0^{a^{+}}{\overline{b}}_N(a,N^{\ast })\overline{P}(a)\mathrm{d}a=\frac{\delta (N^{\ast })}{\kappa }.\end{array}$Let $\overline{\epsilon }(a,N^{\ast })=\overline{b}(a,N^{\ast })+\overline{\sigma }(a)+{\overline{\mu }}_p(a)$. Solving the integral solution of system (53(53) $\{\begin{array}{l}\frac{\overline{x}(a)}{\mathrm{d}a}=-\overline{\lambda }\overline{x}(a)-\overline{b}(a,N^{\ast })\overline{x}(a)-[\overline{\vartheta }(a)+{\overline{\mu }}_p(a)]\overline{x}(a)-{\overline{b}}_N(a,N^{\ast }){\overline{P}}^{\ast }(a)n^{\ast },\\ \frac{\overline{y}(a)}{\mathrm{d}a}=-\overline{\lambda }\overline{y}(a)-{\overline{\mu }}_q\overline{y}(a),\\ \frac{\overline{z}(a)}{\mathrm{d}a}=-\overline{\lambda }\overline{z}(a)-\overline{d}\overline{z}(a),\\ \overline{x}(0)=2(1-\ell )\kappa {\int }_0^{a^{+}}\overline{b}(a,N^{\ast })\overline{x}(a)\mathrm{d}a+2(1-\ell )\kappa \overline{\delta }(N^{\ast })n^{\ast },\\ \overline{y}(0)=2\ell {\int }_0^{a^{+}}\overline{b}(a,N^{\ast })\overline{x}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{x}(a)\mathrm{d}a+2\ell \overline{\delta }(N^{\ast })n^{\ast },\\ \overline{z}(0)=\frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}\overline{y}(a)\mathrm{d}a,\end{array}$(53) ), we get (54) $\begin{array}{rl}& \overline{x}(a)=e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }(\overline{x}(0)e^{-\overline{\lambda }a}-{\overline{P}}^{\ast }(0)n^{\ast }{\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi ),\\ & \overline{y}(a)=\overline{y}(0)e^{-{\int }_0^a(\overline{\lambda }+{\overline{\mu }}_q)\mathrm{d}a},\\ & \overline{z}(a)=\overline{z}(0)e^{-{\int }_0^a(\overline{\lambda }+\overline{d})\mathrm{d}a}.\end{array}$(54) The first equation of (51(51) $\{\begin{array}{l}\frac{\mathrm{\partial }x}{\mathrm{\partial }t}+k\frac{\mathrm{\partial }x}{\mathrm{\partial }a}=-b(a,N^{\ast })x(t,a)-[\vartheta (a)+{\mu }_p(a)]x(t,a)\\ -b_N(a,N^{\ast }){\overline{P}}^{\ast }(a)n(t)-{\mathrm{\Gamma }}_0(t,a),\\ \frac{\mathrm{\partial }y}{\mathrm{\partial }t}+k\frac{\mathrm{\partial }y}{\mathrm{\partial }a}=-{\mu }_{q}y(t,a),\\ \frac{\mathrm{\partial }z}{\mathrm{\partial }t}+k\frac{\mathrm{\partial }z}{\mathrm{\partial }a}=-\mathrm{d}z(t,a),\\ x(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N^{\ast })x(t,a)\mathrm{d}a+2(1-\ell )\delta (N^{\ast })n(t)\\ +2(1-\ell ){\int }_0^{a^{+}}{\mathrm{\Gamma }}_0(t,a)\mathrm{d}a,\\ y(t,0)=\frac{2\ell }{\kappa }{\int }_0^{a^{+}}b(a,N^{\ast })x(t,a)\mathrm{d}a+\frac{1}{\kappa }{\int }_0^{a^{+}}\vartheta (a)x(t,a)\mathrm{d}a\\ +\frac{2\ell }{\kappa }\delta (N^{\ast })n(t)+\frac{2\ell }{\kappa }{\int }_0^{a^{+}}{\mathrm{\Gamma }}_0(t,a)\mathrm{d}a,\\ z(t,0)=\frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}y(t,a)\mathrm{d}a,\end{array}$(51) ) requires the expression of ${\overline{P}}^{\ast }(0)={\overline{P}}^{\ast }(a)e^{{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }$. Substituting the three equations of (54(54) $\begin{array}{rl}& \overline{x}(a)=e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }(\overline{x}(0)e^{-\overline{\lambda }a}-{\overline{P}}^{\ast }(0)n^{\ast }{\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi ),\\ & \overline{y}(a)=\overline{y}(0)e^{-{\int }_0^a(\overline{\lambda }+{\overline{\mu }}_q)\mathrm{d}a},\\ & \overline{z}(a)=\overline{z}(0)e^{-{\int }_0^a(\overline{\lambda }+\overline{d})\mathrm{d}a}.\end{array}$(54) ) into the expression for $n^{\ast }$ yields $\begin{array}{rl}{n}^{\ast }& ={\int }_0^{a^{+}}{e}^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }(\overline{x}(0)e^{-\overline{\lambda }a}-{\overline{P}}^{\ast }(0)n^{\ast }{\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\\ & +\overline{y}(0)e^{-(\overline{\lambda }+{\overline{\mu }}_q)a}+\overline{z}(0)e^{-(\overline{\lambda }+\overline{d})a}\mathrm{d}a.\end{array}$In the same way, substituting $\overline{x}(0)$, $\overline{y}(0)$ and $\overline{z}(0)$ into the above equation, we get $n^{\ast }=\frac{\overline{x}(0)\times ({\int }_0^{a^{+}}{e}^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi -\overline{\lambda }a}\mathrm{d}a+\mathbb{A}(\overline{\lambda })\mathbb{C}(\overline{\lambda }))}{1+{\overline{P}}^{\ast }(0)\mathbb{B}(\overline{\lambda })+\overline{P^{\ast }}(0)\mathbb{A}(\overline{\lambda })\mathbb{D}(\overline{\lambda })},$where $\begin{array}{rl}& \mathbb{A}(\overline{\lambda })={\int }_0^{a^{+}}{e}^{-(\overline{\lambda }+{\overline{\mu }}_q)}\mathrm{d}a\times (1+\frac{{\mu }_q}{\kappa }{\int }_0^{a^{+}}{e}^{-(\overline{\lambda }+\overline{d})a}\mathrm{d}a),\\ & \mathbb{B}(\overline{\lambda })={\int }_0^{a^{+}}(e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }{\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a,\\ & \mathbb{C}(\overline{\lambda })=\frac{\ell }{(1-\ell )\kappa }+{\int }_0^{a^{+}}\overline{\vartheta }(a)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi -\overline{\lambda }a}\mathrm{d}a,\\ & \mathbb{D}(\overline{\lambda })={\int }_0^{a^{\ast }}(\overline{\vartheta }(a)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }{\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a.\end{array}$Let $\mathbb{E}(\overline{\lambda })=\frac{{\int }_0^{a^{+}}{e}^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi -\overline{\lambda }a}\mathrm{d}a+\mathbb{A}(\overline{\lambda })\mathbb{C}(\overline{\lambda })}{1+{\overline{P}}^{\ast }(0)\mathbb{B}(\overline{\lambda })+\overline{P^{\ast }}(0)\mathbb{A}(\overline{\lambda })\mathbb{D}(\overline{\lambda })}.$We get $n^{\ast }=\overline{x}(0)\times \mathbb{E}(\overline{\lambda }).$Then, substituting $n^{\ast }$ into the expression of $\overline{x}(a)$ and $\overline{x}(0)$. We get (55) $\overline{x}(a)=\overline{x}(0)\times e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\times (e^{-\overline{\lambda }a}-{\overline{P}}^{\ast }(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi ),$(55) and (56) $\begin{array}{rl}\overline{x}(0)& =\overline{x}(0)\times [2(1-\ell )\kappa {\int }_0^{a^{+}}(\overline{b}(a,N^{\ast })e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\\ & \times (e^{-\overline{\lambda }a}-{\overline{P}}^{\ast }(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a\\ & +2(1-\ell )\kappa \mathbb{E}(\overline{\lambda })\overline{\delta }(N^{\ast })].\end{array}$(56) Then we get characteristic equation: (57) $\begin{array}{rl}1& =2(1-\ell )\kappa {\int }_0^{a^{+}}(\overline{b}(a,N^{\ast })e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\\ & \times (e^{-\overline{\lambda }a}-\overline{P^{\ast }}(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a+2(1-\ell )\kappa \mathbb{E}(\overline{\lambda })\overline{\delta }(N^{\ast }).\end{array}$(57) Finally, if Equation (57(57) $\begin{array}{rl}1& =2(1-\ell )\kappa {\int }_0^{a^{+}}(\overline{b}(a,N^{\ast })e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\\ & \times (e^{-\overline{\lambda }a}-\overline{P^{\ast }}(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a+2(1-\ell )\kappa \mathbb{E}(\overline{\lambda })\overline{\delta }(N^{\ast }).\end{array}$(57) ) is satisfied, the system solution of the form (22(22) $x(t,a)=e^{\lambda t}\overline{x}(a),y(t,a)=e^{\lambda t}\overline{y}(a),z(t,a)=e^{\lambda t}\overline{z}(a).$(22) ) exists. If all the characteristic roots satisfying Equation (57(57) $\begin{array}{rl}1& =2(1-\ell )\kappa {\int }_0^{a^{+}}(\overline{b}(a,N^{\ast })e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\\ & \times (e^{-\overline{\lambda }a}-\overline{P^{\ast }}(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a+2(1-\ell )\kappa \mathbb{E}(\overline{\lambda })\overline{\delta }(N^{\ast }).\end{array}$(57) ) have only negative real parts, the perturbation form $(x(t,a),y(t,a),z(t,a))$ near the tumour steady state ${\overline{E}}_{\ast }(a)$ tends to $(0,0,0)$ as time tends to positive infinity. The following theorem is obtained.

Theorem 5.2

Let (H1), (H5) and (48(48) $U(N)=1,$(48) ) be satisfied. In addition, if all solutions $\overline{\lambda }$ of Equation (57(57) $\begin{array}{rl}1& =2(1-\ell )\kappa {\int }_0^{a^{+}}(\overline{b}(a,N^{\ast })e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\\ & \times (e^{-\overline{\lambda }a}-\overline{P^{\ast }}(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a+2(1-\ell )\kappa \mathbb{E}(\overline{\lambda })\overline{\delta }(N^{\ast }).\end{array}$(57) ) have negative real parts, the tumour steady state ${\overline{E}}^{\ast }(a)$ of system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) is locally asymptotically stable.

Remark 2

We determine the local stability of the tumour steady state by analysing the sign of the eigenvalue of the characteristic Equation (57(57) $\begin{array}{rl}1& =2(1-\ell )\kappa {\int }_0^{a^{+}}(\overline{b}(a,N^{\ast })e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\\ & \times (e^{-\overline{\lambda }a}-\overline{P^{\ast }}(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a+2(1-\ell )\kappa \mathbb{E}(\overline{\lambda })\overline{\delta }(N^{\ast }).\end{array}$(57) ) that satisfies the tumour steady state perturbation form [38]. Based on our assumptions in this paper, if no solution $\overline{\lambda }$ of (57(57) $\begin{array}{rl}1& =2(1-\ell )\kappa {\int }_0^{a^{+}}(\overline{b}(a,N^{\ast })e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\\ & \times (e^{-\overline{\lambda }a}-\overline{P^{\ast }}(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a+2(1-\ell )\kappa \mathbb{E}(\overline{\lambda })\overline{\delta }(N^{\ast }).\end{array}$(57) ) with $Re\overline{\lambda }\ge 0$, the tumour steady state of the system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) will be locally asymptotically stable. The expression of (57(57) $\begin{array}{rl}1& =2(1-\ell )\kappa {\int }_0^{a^{+}}(\overline{b}(a,N^{\ast })e^{-{\int }_0^a\overline{\epsilon }(\xi ,N^{\ast })\mathrm{d}\xi }\\ & \times (e^{-\overline{\lambda }a}-\overline{P^{\ast }}(0)\mathbb{E}(\overline{\lambda }){\int }_0^a{\overline{b}}_N(\xi ,N^{\ast })e^{-\overline{\lambda }(a-\xi )}\mathrm{d}\xi )\mathrm{d}a+2(1-\ell )\kappa \mathbb{E}(\overline{\lambda })\overline{\delta }(N^{\ast }).\end{array}$(57) ) does not lead to the monotonicity of $\overline{\lambda }$ and it is more difficult to get a mathematical criterion for determining the sign of $\overline{\lambda }$. But, in Section 6, the existence and stability of the tumour steady state is discussed by using numerical simulations.

6. Numerical simulations

We conduct some numerical simulations to illustrate the local and global stabilities of the trivial steady state ${\overline{E}}_0$ as well as the existence and local stability of the positive steady state ${\overline{E}}^{\ast }$. There are four key parameters: the proportion ℓ of newly born tumour cells entering the quiescent state, the evolution speed κ, the dead rate µ, dead cell clearance rate d should be more concerned. Different cells have different maximum survival times, and in the numerical simulations we chose $a^{+}=72$ hours as the maximum tumour survival age and applied it to all examples. The age at which tumour cell proliferation begins is $\overline{a}=50$ hours.

For the nonlinear age-structured tumour proliferation functions we assume a separation variable of the form $b(a,N)=b(a)\mathrm{\Psi }(N)$, where the age-structured proliferation function takes (58) $b(a)=\{\begin{array}{l}0,a\le \overline{a},\\ \frac{1}{\varrho }\frac{(a-\overline{a}{)}^2}{2{\varrho }^2+2\varrho (a-\overline{a})+(a-\overline{a}{)}^2},a>\overline{a}\end{array}$(58) in the interval $a\in [0,a^{+}]$ and the restriction term takes $\mathrm{\Psi }(N)={10}^7/(N+{10}^7)$ for all $N\ge 0$ (see [30,34]), where $\varrho =2$ is a variance constant. Assuming that the mortality rate differs before and after the age at which the tumour begins to proliferate, the mortality rate of the proliferating cells is select as a function as follows. (59) ${\mu }_p(a)=\{\begin{array}{l}\mu ,a\le \overline{a},\\ \mu +\mu (a-\overline{a}{)}^2,a>\overline{a}\end{array}$(59) for all $a\in [0,a^{+}]$. Let the quiescent cells death rate be ${\mu }_q=\mu$. Let (60) $\vartheta (a)=\{\begin{array}{l}\vartheta \frac{a_0-a}{a_0},a<a_0,\\ 0,a\ge a_0,\end{array}$(60) where $\vartheta =0.01$ [30]. $a_0$ is the maximum age of transformation from proliferating cells to quiescent cells. The mitotic interval is generally divided into $G_1$ phase, S phase and $G_2$ phase, but only in $G_1$ phase, the transformation from proliferating cells to quiescent cells occurs. The duration of $G_1$ phase varies from cell to cell in different environments. Billy et al. [39] characterized the cell cycle at the single cell level in proliferating cell populations and differentiated $G_1$ phase cells based on fluorescence image modelling techniques. Basse et al. [40] used flow cytometry to analyse the DNA content of human melanoma cells at multiple time points and gave the proportion of $G_1$ phase cells. For convenience, we take $a_0=40$ hours with the maximum cell survival age as $a^{+}=72$ hours. The forms of b, Ψ, µ and ϑ are based on biological significance and the assumptions of this article.

We study the effect of parameters on the conditions for determining the local stability of the tumour-free steady state, and investigate the global stability of the tumour-free steady state under suitable parameters. Take $\mu =0.002$ [41], the initial functions $P(0,a)=(a^{+}-a)\times {10}^5/a^{+}$, $Q(0)=P(0)$ [35] and $D(0)=0$ [27]. We first study the effect of the proportion ℓ on the condition (26(26) $2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a)e^{-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a<1$(26) ) in Theorem 4.1 for determining the local stability of the tumour-free steady state at different values of κ. Let $J(\mu ,\kappa ,\ell )=2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\overline{\epsilon }(\xi )\mathrm{d}\xi }\mathrm{d}a.$The tumour-free steady state ${\overline{E}}_0$ is locally stable when J<1.

Figure 1(a) shows the effect of the proportion ℓ on J at three different κ values. Note that J<1 always holds for $\kappa =0.5$. The effect of κ and ℓ on J is illustrated in Figure 1(b).

It is known from Figure 1 that for larger values of κ a larger value of ℓ is required to ensure the local stability condition J<1 for the tumour-free steady state. This means that a larger dose of the drug erlotinib may prevent newborn cells from entering the proliferative state. It is interesting to note that when the cell evolution speed κ is less than 0.8717, any proportion ℓ in $[0,1]$ can contribute the inequality J<1, i.e. the extinction of the tumour cells.

Figure 1. (a) Effect of the parameter ℓ on the threshold condition J for stable tumour-free steady state ${\overline{E}}_0$ at three different values of κ. (b) The effect of parameters ℓ and κ on J. ($\mu =0.002$).

Here we study the effect of µ on J.

Figure 2(a) shows the effect of the parameter µ on J for different values of κ. Figure 2(b) shows the effect of parameters µ and ℓ on J for different values of κ. From Figure 2(a ,b) we can conclude that in order to guarantee the condition J<1 with a larger κ, a larger µ and a larger ℓ are required. It can be seen from Figure 2 that the mortality rate µ is smaller in order of magnitude and has a greater effect on J, meaning that the system is more sensitive to the parameter µ.

Now, we study the global stability of the trivial steady state ${\overline{E}}_0$. The initial value functions are chosen to be $P(0,a)=4\times {10}^5\times (a^{+}-a)/a^{+},Q(0)=P(0)$ and $D(0)=0$ in all cases without special instructions. From the previous results (Theorems 4.2 and 4.3) we can know that the global stability of the tumour steady state is independent of clearance rate d. From condition (27(27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) ) in Theorem 4.2 and choosing parameter values $\mu =0.0062$, $\kappa =1$, $\ell =0.3$ respectively, we can obtain $J_1(\mu ,\kappa ,\ell )=2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa })\mathrm{d}\xi }\mathrm{d}a=0.9994<1.$It is easy to see that these parameter values satisfy the global stability condition of ${\overline{E}}_0$. The age distribution of the proliferating cells in Figure 3(a) shows that the number of proliferating cells reaches 0 in less than 500 hours, and Figure 3(b) shows this change from the number of proliferating cells at different ages. Figure 3(c) plots the changes in the number of the three tumour cell subpopulations P, Q and D as well as the changes in the total number of total cells N. Figure 3(d) shows the effect of different initial value functions on the total number of cells N when condition (27(27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) ) is satisfied. It can be seen that the change in the initial value functions is only relevant for the time to reach the tumour-free steady state of the system.

Figure 2. (a) The effect of the parameter µ on J for κ taking three different values. (b) The effect of parameters ℓ and µ on J when κ takes three different values.

Figure 3. (a) Age distribution of proliferating tumour cells $P(t,a)$. (b) Trends in proliferating cells at six different ages. (c) The time series of three types of cells P, Q, D and the total number N. (d) The time series of N under four different initial value functions. The parameters of these four figures take values of $\mu =0.0062$, $\ell =0.3$ and $\kappa =1$, which satisfy the global stability condition (27(27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) ) of Theorem 4.2 for the tumour-free steady state.

Choosing a uniform upper bound $M=2\times {10}^7$ on the total number of cells and by condition (37(37) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{1}{k}{\int }_0^a{\gamma }_{0}b(\xi )+\sigma (\xi )+{\mu }_p(\xi )\mathrm{d}\xi }\mathrm{d}a<1,$(37) ) in Theorem 4.3 we define $J_2(\mu ,\kappa ,\ell )=2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{1}{\kappa }{\int }_0^a({\gamma }_{0}b(\xi )+\vartheta (\xi )+{\mu }_p(\xi ))\mathrm{d}\xi }\mathrm{d}a.$If we choose the parameter values $\kappa =1,\ell =0.3$ and $\mu =0.0043$, by calculation we can get $J_2=0.9940<1$, which satisfies condition (37(37) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{1}{k}{\int }_0^a{\gamma }_{0}b(\xi )+\sigma (\xi )+{\mu }_p(\xi )\mathrm{d}\xi }\mathrm{d}a<1,$(37) ).

Similar to Figures 3 and 4 plots the global stability of the tumour steady state ${\overline{E}}_0$ under a superior condition. Condition (37(37) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{1}{k}{\int }_0^a{\gamma }_{0}b(\xi )+\sigma (\xi )+{\mu }_p(\xi )\mathrm{d}\xi }\mathrm{d}a<1,$(37) ) requires a lower death rate µ and a longer time t to reach the trivial steady state ${\overline{E}}_0$ than condition (27(27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) ). It can be seen that $\mu >0.062$ is needed to satisfy $J_1<1$ if $\kappa =1$ and $\ell =0.3$, yet only $\mu =0.0043$ is needed to satisfy $J_2<1$. Numerical simulations are used to show again that condition (27(27) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-{\int }_0^a\left(\frac{\vartheta (\xi )+{\mu }_p(\xi )}{\kappa }\right)\mathrm{d}\xi }\mathrm{d}a<1,$(27) ) is stronger than condition (37(37) $2(1-\ell ){\int }_0^{a^{+}}b(a)e^{-\frac{1}{k}{\int }_0^a{\gamma }_{0}b(\xi )+\sigma (\xi )+{\mu }_p(\xi )\mathrm{d}\xi }\mathrm{d}a<1,$(37) ).

Figure 4. (a) The age distribution of proliferating cells $P(t,a)$. (b) The time series of $P(t,a)$ at six fixed ages corresponding to Theorem 4.3. (c) The time series of four types of cells P, Q, D and the total number N. (d) The time series of the total number N under four different initial values. Where parameters $\kappa =1$, $\ell =0.3$, and $\mu =0.0043$ are taken. These four figures show the global stability of the trivial steady state ${\overline{E}}_0$ under conditions of Theorem 4.3.

Next we discuss the existence and local stability of the tumour steady state. The values of the parameters are taken as $\mu =0.001,\kappa =1$ and $\ell =0.2$. With these parameter values, the tumour steady state existence condition $U(N)=1$ is satisfied when N reaches $4.5737\times {10}^6$.

Figures 5–7 describe the changes in the number of subpopulations in the tumour steady state from different perspectives. The initial value functions in Figures 5 and 6 is chosen as $P(0,a)=2.5\times {10}^5\times (a^{+}-a)/a^{+}$, $Q(0)=P(0)$ and $D(0)=0$. Figure 5(a) demonstrates the age distribution of proliferating tumour cells $P(t,a)$ under conditions that satisfy the existence of tumour steady state. Figure 5(b) shows the trends of proliferating cells at different ages. Figure 5(c) shows the cross-sections of Figure 5(a) at different times, i.e. the age distribution of proliferating cells at different times. It can be seen that the age distribution of proliferating cells tends to be consistent with increasing time. Figure 6(a) shows the time series of three types of cells P, Q, D and the total number N. From Figure 6(a) we can see that the number of the dead cells D is much lower than the number of other cells, so we show the time series of the number of the death cells D in Figure 6(b) alone.

Figure 5. (a) Age distribution of proliferating cells in tumour steady state. (b) Trends in proliferating cells at different ages. (c) Age distribution of proliferating cells at six different time points.

Figure 6. (a) The time series of three types of cells P, Q, D and the total number N. (b) The time series of the number of the death cells D.

Figure 7. (a) Effect of different initial value functions on proliferating cells of age 0. (b) Time series of N under different initial value functions.

From (43(43) ${\overline{P}}^{\ast }(a)={\overline{P}}^{\ast }(0)e^{-{\int }_0^a\overline{\epsilon }(\xi ,N)\mathrm{d}\xi },$(43) )–( 45(45) $\begin{array}{rl}{\overline{D}}^{\ast }(a)& ={\overline{D}}^{\ast }(0)e^{-\overline{d}a}\\ & ={\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}{\overline{Q}}^{\ast }(\varpi )\mathrm{d}\varpi \\ & ={\overline{P}}^{\ast }(0){\overline{\mu }}_q{e}^{-\overline{d}a}{\int }_0^{a^{+}}[2\ell {\int }_0^{a^{+}}\overline{b}(\xi ,N)e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi \\ & +{\int }_0^{a^{+}}\overline{\vartheta }(\xi )e^{-{\int }_0^{\xi }\overline{\epsilon }(\eta ,N)\mathrm{d}\eta }\mathrm{d}\xi ]e^{-{\overline{\mu }}_q\varpi }\mathrm{d}\varpi .\end{array}$(45) ) we can know that steady state system (17(17) $\{\begin{array}{l}\frac{\mathrm{d}\overline{P}}{\mathrm{d}a}=-\overline{b}(a,N)\overline{P}(a)-\overline{\vartheta }(a)\overline{P}(a)-{\overline{\mu }}_p(a)\overline{P}(a),\\ \frac{\mathrm{d}\overline{Q}}{\mathrm{d}a}=-{\overline{\mu }}_q\overline{Q}(a),\\ \frac{\mathrm{d}\overline{D}}{\mathrm{d}a}=-\overline{d}\overline{D}(a)\end{array}$(17) )–( 18(18) $\begin{array}{rl}\overline{P}(0)& =2\kappa (1-\ell ){\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a,\\ \overline{Q}(0)& =2\ell {\int }_0^{a^{+}}\overline{b}(a,N)\overline{P}(a)\mathrm{d}a+{\int }_0^{a^{+}}\overline{\vartheta }(a)\overline{P}(a)\mathrm{d}a,\\ \overline{D}(0)& ={\overline{\mu }}_q{\int }_0^{a^{+}}\overline{Q}(a)\mathrm{d}a,\end{array}$(18) ) depends continuously on the number of the newborn proliferating cells with age zero. We study the effect of different initial value functions on the total number of cells and the number of proliferating cells at age 0.

It can be seen from Figure 7(a ,b) that even though the initial value functions are different, the total number of cells N and the proliferating cells with age 0 tend to be the same for a long enough time, which indicates the stability of the tumour steady state. It shows that the initial value functions have no effect on the steady state of the system at a long enough time.

Figures 8 and 9 show the effect of parameters κ, ℓ and d on the number of P, Q, D and N at sufficiently long time points t = 4500 hours under the same mortality rate $\mu =0.001$, respectively. Figure 7, we can know that the system is stable when the time is taken as 4500 hours. From Figure 8(a) we can observe that all tumour subpopulations P, Q, D and the total number of tumour cells N in the steady state increase with increasing κ. This may indicate that larger solid spherical tumours have a greater rate of cell evolution κ. Figure 8(b) shows the number of cells P, Q, D and N strictly decreases with the increase in the parameter ℓ.

Figure 8. (a) The effect of increasing the parameter κ at a sufficiently large time point $t_0=4500h$ on P, Q, D and N. (b) The effect of increasing the parameter ℓ on P, Q, D and N at the same time. ($\mu =0.001$, $\ell =0.2$ and d = 0.01).

Figure 9. (a)–(b) The effect of the increase of d on P, Q, D and N when the system reaches the steady state ($t_0=4500h$). The parameter values are chosen as $\mu =0.001$, $\ell =0.2$ and $\kappa =1$.

The tumour steady state is transformed into the tumour-free steady state if ℓ is large enough. We can also learn that $J_2\le 1$ is satisfied when $\kappa =1$, $\mu =0.001$, d = 0.001, $M=2\times {10}^7$, $\ell >0.6343$, i.e. system (4(4) $\{\begin{array}{l}\frac{\mathrm{\partial }P}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }P}{\mathrm{\partial }a}=-b(a,N(t))P(t,a)-\vartheta (a)P(t,a)-{\mu }_p(a)P(t,a),\\ \frac{\mathrm{\partial }Q}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }Q}{\mathrm{\partial }a}=-{\mu }_{q}Q(t,a),\\ \frac{\mathrm{\partial }D}{\mathrm{\partial }t}+\kappa \frac{\mathrm{\partial }D}{\mathrm{\partial }a}=-\mathrm{d}D(t,a)\end{array}$(4) )–( 6(6) $\{\begin{array}{l}P(t,0)=2(1-\ell ){\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a,\\ Q(t,0)=\frac{1}{\kappa }[2\ell {\int }_0^{a^{+}}b(a,N(t))P(t,a)\mathrm{d}a+{\int }_0^{a^{+}}\vartheta (a)P(t,a)\mathrm{d}a],\\ D(t,0)=\frac{1}{\kappa }[{\mu }_q{\int }_0^{+\mathrm{\infty }}Q(t,a)\mathrm{d}a].\end{array}$(6) ) has a globally stable tumour-free steady state. From Figure 9(a ,b) we can obtain that with the increase of the clearance rate d, Q and P will initially increase, while D and N will initially decrease. The increase of the clearance rate d initially had a significant effect on the number of dead cells D (see Figure 9(b)). When the clearance rate d>0.7, the change of the number of three kind cells P,Q,D and the total cells' number N will be very small.

We know from numerical simulations that controlling solid tumour growth requires a smaller κ, as well as a larger ℓ. This means that controlling tumour growth requires a smaller rate of cell evolution and a larger proportion of newly born tumour cells that allowed to enter the quiescent state. The initial value functions affect the arrival time to the steady state, but have no effect on the system when the time is long enough. The dead cell clearance rate d does not have a significant impact on the system than κ and ℓ, while D is still sensitive to the initial change in parameter d.

7. Discussion

We analysed a multilayer age-structured tumour model with nonlinear proliferation functions and necrotic cores. First, we discussed some fundamental results such as the existence, uniqueness (Theorem 3.2) and boundedness (Theorem 3.2) of solutions of the model. Then we obtained local stability conditions (Theorem 4.1) for the tumour-free steady state under some assumptions. The global stability of the tumour-free steady state is proved under two different assumptions (Theorems 4.2 and 4.3) using the comparability principle. The existence (Theorem 5.1) and local stability (Theorem 5.2) of the tumour steady state have also been demonstrated. Numerical simulations visualize the effect of parameter changes on the threshold conditions of the tumour-free steady state and analyse the effect of initial values and parameter changes on the system.

Our three-component model and mathematical theoretical results can provide some suggestions on these tumours with avascular growth. A greater proportion ℓ and a greater death rate µ can reduce tumour size and even eliminate tumours. A smaller evolution speed κ contributes to tumour control. Moreover, The death cell clearance rate d does not affect the extinction of the tumour population. Compared with the results in article [28], the dead cell clearance rate d has an effect on the tumour size in the positive steady state, but the effect of this is very small.

In this paper, the selection of the parameters is based on literature but not distortion. If supported by experimental data, we may estimate the parameter values, and hence can greatly improve the application of this paper. This will be an important work we will pursue in the future.

This model does not consider the spatial factors and the bidirectional transformation of the proliferating and quiescent states. Considering spatial factors and bidirectional transfer into the model is very challenging and meaningful problem. In addition, it would be interesting and worthwhile to consider the four stages of the proliferating cell cycle and tumour immunity into future work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was supported by the National Natural Science Foundation of China [grant number 12271068], the Innovation Development Joint Funds of the Chongqing Natural Science Foundation [grant number 2022NSCQ-LZX0360], Chongqing Municipal Education Commission [grant number KJQN202000730], Joint Training Base Construction Project for Graduate Students in Chongqing [grant number JDLHPYJD2021016] and Group Building Scientific Innovation Project for universities in Chongqing [grant number CXQT21021].