Global dynamics of a tuberculosis model with fast and slow progression and age-dependent latency and infection

ABSTRACT In this paper, a mathematical model describing tuberculosis transmission with fast and slow progression and age-dependent latency and infection is investigated. It is assumed in the model that infected individuals can develop tuberculosis by either of two pathogenic mechanisms: direct progression or endogenous reactivation. It is shown that the transmission dynamics of the disease is fully determined by the basic reproduction number. By analyzing corresponding characteristic equations, the local stability of a disease-free steady state and an endemic steady state of the model is established. By using the persistence theory for infinite dimensional system, it is proved that the system is uniformly persistent when the basic reproduction number is greater than unity. By constructing suitable Lyapunov functionals and using LaSalle's invariance principle, it is verified that the global dynamics of the system is completely determined by the basic reproduction number.


Introduction
Tuberculosis (TB) is a bacterial disease caused by Mycobacterium tuberculosis, and is usually acquired through airborne infection from active TB cases [4,9]. According to the World Health Organization, one third of the world's population is infected with tuberculosis either latently or actively. Despite effective antimicrobial chemotherapy, tuberculosis infection remains a leading cause of death from an infectious disease [34].
Mathematical modelling has been proven to be important in better understanding the transmission dynamics of TB and evaluating the effectiveness of various control and prevention strategies (see, for example, [1][2][3][5][6][7][10][11][12][18][19][20][21]25,29,30,32]). Tuberculosis is most commonly transmitted from a person suffering from infectious (active) tuberculosis to other persons by infected droplets created when the person with active TB coughs or sneezes. Among generally healthy persons, infection with TB is highly likely to be asymptomatic [10]. In [6], Blower et al. formulated and analysed mathematical models describing the transmission dynamics of untreated tuberculosis epidemics. It was assumed in [6] that infected individuals remain noninfectious until they develop disease by one of two pathogenic mechanisms: direct progression or endogenous reactivation. Two types of tuberculosis contribute to the incidence rate of tuberculosis disease: one type of tuberculosis develops through direct progression soon after infection (fast progression) and the other type of tuberculosis develops through endogenous reactivation in the latently infected individuals (slow progression). In order to assess the intrinsic transmission dynamics of tuberculosis, Blower et al. considered the following mathematical model [6]: In (1), S(t) represents the number of individuals who are not previously exposed to tuberculosis at time t, L(t) represents the number of latently infected individuals who have been infected with Mycobacterium tuberculosis at time t, but have no clinical illness and remain noninfectious, I(t) represents the number of active infectious tuberculosis cases at time t. The assumptions of model (1) were made as follows [6].
(A1) Recruitment to the susceptible population occurs at a constant rate A. The incidence rate of infection is formulated as the product of the number of susceptible population at time t and the per capita force of infection at time t, where the per capita infection force is defined as the per-susceptible risk of becoming infected with Mycobacterium tuberculosis, and is calculated as the product of the number of infectious cases at time t and the transmission coefficient (β) of the pathogen. Where β represents the likelihood that an infectious case will successfully transmit the infection to a susceptible individual. (A2) The two pathogenous mechanisms are assumed by allowing a proportion (p) of the newly infected to develop tuberculosis directly and a proportion (1−p) of the newly infected to enter the latent class. Latently infected individuals either develop tuberculosis slowly at an average rate ν or die at a natural death rate μ before developing tuberculosis. (A3) Active infectious tuberculosis die, either because of tuberculosis at average rate μ T , or because of natural death at an average rate μ.
The qualitative dynamics of model (1) are governed by the basic reproduction number: Where the quantity R 0 represents the average number of secondary infectious cases produced by a typical infectious individual in a susceptible population at a demographic steady state. The first term in (2) gives the new cases resulting from fast progression while the second term represents the cases resulting from slow progression. We note that system (1) was formulated as ordinary differential equations with distinct variables describing the population size of compartments such as susceptible, latently infected and active infectious. It is assumed that all individuals within a compartment behave identically, regardless of how much time they have spent in the compartment. However, an infected individuals may have latent tuberculosis for months, years or even decades before the disease becomes infectious. The risk per unit time of activation appears to be higher in the early stages of latency than in later stages [6]. On the other hand, laboratory studies suggest that the infectivity of infectious individuals be different at the differential age of infection [37]. For tuberculosis infection, the TB bacteria need to develop in the lung to be transmissible through coughing, and their transmissibility depends on their progression in the lung as well as the strength of a host's immune system. Active TB has the highest possibility of developing within the first 2-5 years of infection, while most TB infections remain latent for a long period of time until immune compromise occurs due to aging or co-infection with other illnesses such as HIV (see, for example, [14,25]). In [28], by including the duration that an individual has spent in the exposed compartment and in the infectious compartment as variables, McCluskey presented and studied an SEI epidemiological model with continuous age-structure in the exposed and infectious classes. A Lyapunov functional was used in [28] to show that the endemic steady state is globally stable. Recently, much attention has been paid to the modelling and analysis on infectious disease dynamics with age dependence (see, for example, [8,13,[15][16][17]27,36]).
Motivated by the works of Blower et al. [6] and McCluskey [28], in the present paper, we are concerned with the effect of age structure for both latently infected and active infectious individuals on the transmission dynamics of untreated tuberculosis epidemics with fast and slow progression. To this end, we consider the following differential equation system: with boundary conditions and initial condition In system (3), S(t) represents the number of individuals who are susceptible to tuberculosis at time t; e(θ , t) represents the density of individuals who have been in the latently infected class for duration θ at time t, but have no clinical illness and remain noninfectious; i(a, t) represents the density of individuals who have been in the active infectious class for duration a at time t. The definitions of all parameters and symbols in system (3) are listed in Table 1.
In the sequel, we further make the following assumptions. The rate at which individuals who have been in the latently infected class for duration θ are removed ε(θ) The rate at which individuals who have been in the latently infected class for duration θ progress to the active infectious class β(a) The transmission coefficient of active infectious individuals at age of infection a ν(a) The natural and disease-induced death rate of individuals who have been in the active infectious class for duration a p The proportion of the newly infected to develop tuberculosis directly 1−p The proportion of the newly infected to enter the latent class (H1) β and ε are Lipschitz continuous on R + with Lipschitz coefficients L β and L ε , respectively. (H2) β, ε, μ, ν ∈ L ∞ + (0, ∞) with essential upper boundsβ,ε,μ andν, respectively. (H3) There exists μ 0 ∈ (0, μ S ] such that μ(a) ≥ μ 0 for a ≥ 0 and ν(θ) ≥ μ 0 for θ ≥ 0. (H4) There exist positive constants a β and a ε such that β(a) is positive in a neighbourhood of a β and ε(a) is positive in a neighbourhood of a ε .
Using the theory of age-structured dynamical systems developed in [24,35], we can show that system (3) has a unique solution (S(t), e(·, t), i(·, t)) satisfying the boundary conditions (4) and the initial condition (5). Moreover, it is easy to show that all solutions of system (3) with the boundary conditions (4) and the initial condition (5) are defined on [0, +∞) and remain positive for all t ≥ 0. Furthermore, X is positively invariant and system (3) exhibits a continuous semi-flow : R + × X → X , namely, Given a point (x, ϕ, ψ) ∈ X , we have the norm (x, ϕ, ψ) The primary goal of the present work is to carry out a complete mathematical analysis of system (3) with the boundary conditions (4) and the initial condition (5), and establish its global dynamics. The organization of this paper is as follows. In the next section, we show the asymptotic smoothness of the semi-flow generated by system (3). In Section 3, we calculate the basic reproduction number and discuss the existence of feasible steady states of system (3). In Section 4, by analyzing corresponding characteristic equations, we establish the local asymptotic stability of a disease-free steady state and an endemic steady state of system (3). In Section 5, by means of the persistence theory for infinite dimensional system, we show that system (3) is uniformly persistent when the basic reproduction number is greater than unity. In Section 6, by constructing suitable Lyapunov functionals and using LaSalle's invariance principle, we study the global stability of each of feasible steady states of system (3) with the boundary conditions (4). In Section 7, numerical simulations are carried out to illustrate the feasibility of theoretical results. A brief discussion is given in Section 8 to conclude this work.

Asymptotic smoothness
In order to address the global dynamics of system (3), in this section, we are concerned with the asymptotic smoothness of the semi-flow { (t)} t≥0 generated by system (3).

Boundedness of solutions
In this subsection, we study the boundedness of solutions of system (3) with the boundary conditions (4) and the initial condition (5). Proposition 2.1: Let t be defined as in (6). Then the following statements hold true.

is point dissipative: there is a bounded set that attracts all points in X .
Proof: Let t (X 0 ) = (t, X 0 ) := (S(t), e(·, t), i(·, t)) be any nonnegative solution of system (3) with the boundary conditions (4) and the initial condition (5).
The following results are direct consequences of Proposition 2.1.

Proposition 2.2:
If X 0 ∈ X and X 0 X ≤ K for some K ≥ A/μ 0 , then for all t ≥ 0.

Proposition 2.3: Let C ∈ X be bounded. Then
(1) t (C) is bounded for all t ∈ R + ; (2) t is eventually bounded on C.

Asymptotic smoothness
In order to investigate the global dynamics of system (3), in this subsection, we show the asymptotic smoothness of the semi-flow { (t)} t≥0 . Let (S(t), e(·, t), i(·, t)) be a solution of system (3) with the boundary conditions (4) and the initial condition (5). Integrating the second and the third equations of system (3) along the characteristic line t − a =const., respectively, one has and where and in which

Proposition 2.4:
The functions A 1 (t) and A 2 (t) are Lipschitz continuous on R + . Then On substituting (10) into (13), it follows that By Proposition 2.2, we have that L 2 (t) ≤ (pβK +ε)K. Noting that φ 2 (a) ≤ 1, it follows from (14) that From (10) we have for all a ≥ 0, t ≥ 0, h ≥ 0. Hence, (15) can be rewritten as Noting that 1 − e −x ≤ x for x ≥ 0, it follows from (17) that here the fact that β(a) is Lipschitz continuous on R + was used.
In a similar way, we have This completes the proof.

Proposition 2.5:
The functions L 1 (t) and L 2 (t) are Lipschitz continuous on R + . Then In a similar way, one has This completes the proof.
In order to prove the asymptotic smoothness of the semi-flow generated by system (3), we introduce the following theorems (Theorems 2.46 and B.2 in [31]).

Theorem 2.1:
The semi-flow : R + × X + → X + is asymptotically smooth if there are maps , : x) and the following hold for any bounded closed set C ⊂ X + that is forward invariant under : (2) there exists t C ≥ 0 such that (t, C) has compact closure for each t ≥ t C .

Theorem 2.2:
Let C be a subset of L 1 (R + ). Then C has compact closure if and only if the following assumptions hold: We are now in a position to state and prove a result on the asymptotic smoothness of the semi-flow { (t)} t≥0 generated by system (3).

Theorem 2.3:
The semi-flow generated by system (3) is asymptotically smooth.
Proof: To verify the conditions (1) and (2) in Theorem 2.1, we first decompose the semi-flow into two parts: Clearly, we have = + for t ≥ 0. Let C be a bounded subset of X and K > A/μ 0 the bound for C.
From Proposition 2.2 we see that S(t) remains in the compact set [0, K]. Next, we show thatẽ(θ , t) andĩ(a, t) remain in a pre-compact subset of L 1 + independent of X 0 . It is easy to show thatẽ(θ , t) ≤L 1 e −μ 0 θ ,ĩ(a, t) ≤L 2 e −μ 0 a , wherē Therefore, the assumptions (i),(ii) and (iv) of Theorem 2.2 follow directly. We need only to verify that (iii) of Theorem 2.2 holds. Since we are concerned with the limit as h → 0, we assume that h ∈ (0, t). In this case, we have It follows from (19) and (20) that In a similar way, we have Hence, the condition (iii) of Theorem 2.2 holds. By Theorem 2.1, the asymptotic smoothness of the semi-flow generated by system (3) follows. This completes the proof.
The following result is immediate from Theorem 2.33 in [31] and Theorem 2.3.

Theorem 2.4:
There exists a global attractor A of bounded sets in X .

Steady states and basic reproduction number
In this section, we are concerned with the existence of feasible steady states of system (3) with the boundary conditions (4).
Clearly, system (3) always has a disease-free steady state E 1 (A/μ S , 0, 0). If system (3) with the boundary conditions (4) has an endemic steady state (S * , e * (θ ), i * (a)), then it must satisfy the following equations: We obtain from the second and the third equations of (21) and (11) that and It follows from the fourth and the fifth equations of (21), (22) and (23) that On substituting (24) into the first equation of (21), we have where R 0 is called the basic reproduction number representing the average number of new infections generated by a single newly infectious individual during the full infectious period [33]. It follows from the fourth equation of (21), (23) and (25) that In conclusion, if R 0 > 1, in addition to the disease-free steady state E 1 (A/μ S , 0, 0), system (3) with the boundary conditions (4) has a unique endemic steady state

Local stability
In this section, we study the local stability of the disease-free and endemic steady states of system (3) wih the boundary conditions (4). We first consider the local stability of the disease-free steady state 1 (a, t). Linearizing system (3) at the steady state E 1 , we obtain from (3) and (4) thaṫ 1 (a, t), Looking for solutions of system (26) of the form x 1 (t) = x 11 e λt , y 1 (θ , t) = y 11 (θ ) e λt , z 1 (a, t) = z 11 (a) e λt , where x 11 , y 11 (θ ) and z 11 (a) will be determined later, one obtains the following linear eigenvalue problem: Solving the second and the third equations of system (27) yields and z 11 (a) = z 11 (0) e − a 0 (λ+ν(s)) ds .
On substituting (28) into the fifth equation of system (27), we have that It follows from (30) and the fourth equation of system (27) that On substituting (29) into (31), we obtain the characteristic equation of system (3) at the steady state E 1 of the form: where is a decreasing function. Therefore, if R 0 > 1, Equation (32) has a unique positive root. Accordingly, if R 0 > 1, the steady state E 1 is unstable.
In the following, we claim that if R 0 < 1, the steady state E 1 is locally asymptotically stable. Otherwise, Equation (32) has at least one root λ 1 = a 1 + ib 1 satisfying a 1 ≥ 0. In this case, we have that a contradiction. Therefore, if R 0 < 1, all roots of Equation (32) have negative real parts. Accordingly, the steady state E 1 is locally asymptotically stable if R 0 < 1.
We are now in a position to study the local stability of the endemic steady state E * (S * , e * (θ ), i * (a)) of system (3) with the boundary conditions (4) in the case that R 0 > 1. Letting , and linearizing system (3) at the steady state E * , we have thaṫ (a, t), Looking for solutions of system (33) of the form x(t) = x 2 e λt , y(θ , t) = y 2 (θ ) e λt , z(a, t) = z 2 (a) e λt , where x 2 , y 2 (θ ) and z 2 (a) will be determined later, we obtain the following linear eigenvalue problem: Solving the second and the third equations of system (34) yields and z 2 (a) = z 2 (0) e − a 0 (λ+ν(s)) ds .
It follows from the first equation of system (34) that On substituting (37) into the fourth and the fifth equations of system (34), we have that and It follows from (38) and (39) that On substituting (36) into (40), one obtains the characteristic equation of system (3) at the steady state E * of the form where In the following, we verify that if R 0 > 1, all roots of Equation (41) have real negative parts. Otherwise, Equation (41) has at least one root λ 2 = a 2 + b 2 i satisfying a 2 ≥ 0. In this case, we have a contradiction. Hence, if R 0 > 1, the steady state E * of system (3) is locally asymptotically stable. From what has been discussed above, we have the following result. (5), if R 0 < 1, the disease-free steady state E 1 (A/μ S , 0, 0) is locally asymptotically stable; if R 0 > 1, E 1 is unstable and the endemic steady state E * (S * , e * (θ ), i * (a)) exists and is locally asymptotically stable.

Uniform persistence
In this section, we establish the uniform persistence of the semi-flow { (t)} t≥0 generated by system (3) when the basic reproduction number R 0 > 1. Defineā With the Assumption (H4), we haveā > 0,θ > 0. Denote Following [26], we have the following result. The following result is useful in proving the uniform persistence of the semi-flow { (t)} t≥0 generated by system (3).
Furthermore, we have the following result.

Proposition 5.2:
The semi-flow { (t)} t≥0 has a compact global attractor A 0 ⊂ Y, which attracts any bounded sets of Y + .

Global stability
In this section, we are concerned with the global asymptotic stability of each of feasible steady states of system (3) with the boundary conditions (4). The strategy of proofs is to use LaSalle's invariance principle.
We first give a result on the global stability of the disease-free steady state E 1 (A/μ S , 0, 0) of system (3). e(θ , t), i(a, t)) be any positive solution of system (3) with the boundary conditions (4) and the initial condition (5).

Theorem 6.1:
where the positive constants k 1 and k 2 and the nonnegative kernel functions F 1 (θ ) and F 2 (a) will be determined later. Calculating the derivative of V 1 (t) along positive solutions of system (3) with the boundary conditions (4), we obtain that (a, t) da Using integration by parts, we have from (58) that (a, t) da We choose By direct calculations, one obtains that It follows from (59)-(61) that On substituting (4) into Equation (62), we have that (a, t) da Choosing it follows from (63) that Clearly, if R 0 < 1, then V 1 (t) ≤ 0 holds and V 1 (t) = 0 implies that S(t) = S 0 , e(θ , t) = 0, i(a, t) = 0. Hence, the largest invariant subset of {V 1 (t) = 0} is the singleton E 0 1 (S 0 , 0, 0). By Theorem 4.1, we see that if R 0 < 1, the steady state E 1 is locally asymptotically stable. Therefore, the global asymptotic stability of E 1 follows from LaSalle's invariance principle. This completes the proof.
We are now in a position to study the global asymptotic stability of the endemic steady state E * (S * , e * (θ ), i * (a)) of system (3) with the boundary conditions (4). Denote In order to guarantee the Lyapunov functional in proving the global stability of E * to be well-defined, we make the following assumption: ∞ 0 e − θ 0 (μ(s)+ε(s)) ds ln e 0 (θ ) dθ < +∞. e(θ , t), i(a, t)) be any positive solution of system (3) with the boundary conditions (4) and the initial condition (5). It follows from (10) that

Lemma 6.1: If (H5) holds, then we have
where L 2 (t) is defined in (12). From Proposition 2.2 we see that L 2 (t) is bounded.
Hence, it follows from (66) that Proof: Let (S(t), e(θ , t), i(a, t)) be any positive solution of system (3) with the boundary conditions (4) and the initial condition (5). Define where the function G(x) = x − 1 − ln x for x > 0, and the positive constants c 1 and c 2 will be determined later. If (H5) holds, From Lemma 6.1 we see that all integrals involved in V 2 (t) are finite. Calculating the derivative of V 2 (t) along positive solutions of system (3) with the boundary conditions (4), we have On substituting A = μ S S * + S * ∞ 0 β(a)i * (a) da, ∂e(θ , t)/∂t = −(μ(θ ) + ε(θ))e(θ , t) − ∂e(θ , t)/∂θ and ∂i(a, t)/∂t = −ν(a)i(a, t) − ∂i(a, t)/∂a into Equation (69), it follows that Noting that and we have and We obtain from (70), (73) and (74) that Using integration by parts, it follows from (75) that We derive from (71), (72) and (76) that It is easy to show that We obtain from (77) and (78) that It follows from (4) and (79) that Choosing we have from (24) and (80) that where Noting that one has (84) Noting that It therefore follows from (81), (84) and (85) that Since the function G(x) = x − 1 − ln x ≥ 0 for all x > 0 and G(x) = 0 holds iff x = 1. Hence, V 2 (t) ≤ 0 holds if R 0 > 1. It is readily seen from (86) that V 2 (t) = 0 if and only if for all θ ≥ 0, a ≥ 0. It can be verified that the largest invariant subset of {V 2 (t) = 0} is the singleton E * . By Theorem 4.1, if R 0 > 1, E * is locally asymptotically stable. Therefore, using LaSalle's invariance principle, we see that if R 0 > 1, the global asymptotic stability of E * follows. This completes the proof.

Numerical simulation
In this section, we give some numerical simulations to illustrate the theoretical results in Sections 2 and 3. The backward Euler and linearized finite difference method will be used to discretize the ODEs and PDE in system (3), and the integral will be numerically calculated using Simpson's rule. The transmission coefficient of active infectious individuals at age of infection a is chosen as: The other parameters in system (3) The initial condition is chosen as S(0) = 1.4 × 10 7 , e 0 (θ ) = e −θ , i 0 (a) = e −a . By calculation, we obtain the basic reproduction number R 0 = 2.4827. It is easy to see that in addition to the infection-free steady state E 1 (1.4 × 10 7 , 0, 0), system (3) has an endemic steady state E * (5.6389 × 10 6 , 2.0014 × 10 4 φ 2 (a), 1.1347 × 10 5 φ 1 (θ )) which is locally asymptotically stable. Numerical simulation illustrates this fact (see Figure 1) (we plotted the i(a, t) component only, denote I(t) = 75 0 i(a, t) da).

Discussion
In this work, a mathematical model describing tuberculosis transmission with fast and slow progression and age-dependent latency and infection was investigated. By calculations, the basic reproduction number was obtained. It has been shown that the global dynamics of system (3) with boundary conditions (4) and initial condition (5) is completely determined by the basic reproduction number. By constructing suitable Lyapunov functionals and using LaSalle's invariance principle, it has been verified that if the basic reproduction number is less than unity, the disease-free steady state is globally asymptotically stable and the disease dies out; and if the basic reproduction number is greater than unity, the endemic steady state is globally asymptotically stable and the disease persists. The global stability of the endemic steady state rules out any possibility for the existence of Hopf bifurcations and sustained oscillations in system (3) with boundary conditions (4).
In the following, we present some special cases of system (3) with the boundary conditions (4).
The global dynamics of system (3) with the boundary conditions (90) has been established in [28] by constructing suitable Lyapunov functionals and using LaSalle's invariance principle. It has been shown from Theorems 7.1 and 9.5 in [28] that system (3) with the boundary conditions (90) admits the same global dynamics as system (3) with the boundary conditions (4).
It therefore follows from (3) Hence, system (3) with the boundary conditions (4) and the initial condition (5) reduces to system (1). The global dynamics of system (1) was investigated in [22] by means of Lyapunov functions and LaSalle's invariance principle. It has been verified in [22] that system (1) has the same global dynamics as system (3) with the boundary conditions (4) and the initial condition (5). From the expression of the basic reproduction number R 0 of system (3), we see that the introduction of the age structure has significant influence on the value of the basic reproduction number of tuberculosis transmission. In fact, a direct calculation shows that with the assumption (91), R 0 reduces to R 0 , the basic reproduction number of epidemiological model (1). We note that it is assumed in system (1) that all parameters are positive constants and all individuals within a compartment behave identically, regardless of how much time they have spent in the compartment. For instance, infectious individuals are assumed to be equally infectious during their periodic infectivity and the waiting times in each compartment are assumed to be exponentially distributed. In system (3), we included the duration that an individual has spent in the latently infected compartment and in the active infectious compartment as continuous variables. Hence, R 0 can be used to describe the number of new infections generated by a single newly infectious individual during the full infectious period more accurately than R 0 . This indicates that although ODE models without age structure are usually easier to use than PDE models with age structure, it could lead to the over-valuation or underestimation of the basic reproduction number of tuberculosis transmission.
By choosing appropriate kernel functions, system (3) also contains epidemiological models with time delay, and the global stability results in this work provide the global dynamics for these delayed epidemic models.