key: cord-0698078-4xj45teg
authors: Guo, Zhong-Kai; Xiang, Hong; Huo, Hai-Feng
title: Analysis of an age-structured tuberculosis model with treatment and relapse
date: 2021-04-02
journal: J Math Biol
DOI: 10.1007/s00285-021-01595-1
sha: b7b4438c82b32624584037147460fe3341afe542
doc_id: 698078
cord_uid: 4xj45teg

A new tuberculosis model consisting of ordinary differential equations and partial differential equations is established in this paper. The model includes latent age (i.e., the time elapsed since the individual became infected but not infectious) and relapse age (i.e., the time between cure and reappearance of symptoms of tuberculosis). We identify the basic reproduction number [Formula: see text] for this model, and show that the [Formula: see text] determines the global dynamics of the model. If [Formula: see text] , the disease-free equilibrium is globally asymptotically stable, which means that tuberculosis will disappear, and if [Formula: see text] , there exists a unique endemic equilibrium that attracts all solutions that can cause the spread of tuberculosis. Based on the tuberculosis data in China from 2007 to 2018, we use Grey Wolf Optimizer algorithm to find the optimal parameter values and initial values of the model. Furthermore, we perform uncertainty and sensitivity analysis to identify the parameters that have significant impact on the basic reproduction number. Finally, we give an effective measure to reach the goal of WHO of reducing the incidence of tuberculosis by 80% by 2030 compared to 2015.

Tuberculosis (TB), one of the top 10 causes of death worldwide, is an infectious disease caused by the bacillus mycobacterium tuberculosis. The bacillus can attack any part of the body such as the lungs (pulmonary TB), kidney, spine and brain (extrapulmonary TB). The disease is spread through the air from one person to another, when people who are sick with pulmonary TB expel the bacillus into the air, for example, by coughing, speaking and so on, people nearby may breathe in these bacillus and become infected. The extrapulmonary TB is usually not infectious. About 5-10% of people who infected with the bacillus will develop TB during their lifetime. However, people with weak immune systems have a much higher risk of TB, such as people living with HIV, malnutrition or diabetes, or people who smoke. TB occurs in every part of the world. In 2017, there were an estimated 10 million new cases of TB worldwide, the majority of them were concentrated in the African, South-East Asia and Western Pacific regions. The following eight countries accounted for two thirds of the new TB cases: India, China, Indonesia, the Philippines, Pakistan, Nigeria, Bangladesh and South Africa (WHO 2019; CDC 2019).

Prevention and treatment have always been considered as important means of controlling the spread of TB, they mainly include: treatment of latent TB infection, vaccination of children with the bacille Calmette-Guérin (BCG) vaccine, people with TB take the drugs exactly as prescribed for 6-9 months, and so on. These measures are considered to have made a major contribution to the reduction in TB burden since 2000. TB incidence has fallen 1.5% per year since 2000, for a total reduction of 18% worldwide. However, as the WHO notes, TB kills millions of people worldwide each year, which is unacceptable in an era when you can diagnose and cure nearly every person with TB. At this rate, the WHO 2030 goal (the new TB cases will reduce by 80% by 2030 compared to 2015) will not be achieved. If the world is to end TB, it needs to scale up services and invest in research. Unfortunately, the gap between the funding available and the requirement has widened over the past 2 years, this is not a good sign for TB control. Therefore, increasing investment and developing the most effective control measures are crucial for the current control of TB (WHO 2019; CDC 2019) .

Mathematical modeling has always been a very useful and important tool in studying the transmission mechanism of epidemic disease, drinking or other behavior, through mathematical analysis and numerical simulation of the model, effective control measures can be designed and evaluated (Ullah et al. 2019; Wang et al. 2019; Agusto and Khan 2018; Rebaza 2019; Ngeleja et al. 2018; Zhao et al. 2019; Huo et al. 2019b; Ma 2019; Shi and Wang 2019; Cai et al. 2018; Du and Feng 2018; Zhang and Xing 2019; Meng and Wu 2018; Zhang et al. 2018; Huang and Tian 2019; Li and Huang 2019; Huo et al. 2019a; Zhang et al. 2019a, b; Huo et al. 2020; Huang 2020) . Due to the complexity of the pathogenesis of tuberculosis, there have been many articles to establish different models to study TB and incorporate certain factors, such as treatment, drug resistance, HIV/TB co-infection, relapse, vaccination and so on (Whang et al. 2011; Zhao et al. 2017; Ozcaglar et al. 2012; Liu and Zhang 2011; Ren 2017; Choi et al. 2015) . Liu and Zhang (2011) argued that vaccination can only reduce infection, not eliminate it entirely, and treatment can also only reduce the infectivity of treated TB cases (compared to untreated cases). According to this, they investigated the effects of vaccination and treatment on the spread of TB. Agusto and Adekunle (2014) formulated a mathematical model of the transmission of TB-HIV/AIDS co-infection. By using treatment of infected individuals with TB as the system control variable, they investigated optimal strategies for controlling the spread of the disease. Most of the models are qualitative analysis of tuberculosis transmission, and few are quantitative analysis based on actual data. Recently, more and more researchers are beginning to combine data to study the spread of diseases (Jing et al. 2020; Ullah et al. 2019; Zhao et al. 2017; Ghosh et al. 2018; Liu et al. 2020) . Ullah et al. (2019) studied a TB model and used the actual data to estimate the model parameters, then they analyzed the effect of various parameters to define the importance of these parameters to disease transmission. Based on tuberculosis data from China, Zhao et al. (2017) evaluated the parameters by the Least Square method. Then, they gave two effective measures to achieve WHO's goal through numerical simulation.

During the spread of tuberculosis, it is worth noting that after becoming infected, some people develop TB disease soon (within weeks), before their immune system can fight the TB bacteria. Other people may get sick years later, when their immune system becomes weak for some reason. Many people with TB infection never develop TB disease. Iannelli and Milner (2017) argued that age structure is important when modeling a long-lasting disease, because the chance of infected individuals becoming infectious individuals and recovering patients relapsing is variable.

Motivated by these factors, in this paper, we formulate a new age-structured TB model to study the effects of relapse and treatment on transmission dynamics of TB. We not only introduce the treatment individuals as a class but also assume that the treatment class is also infectious. The article is organized as follows. In Sect. 2, we introduce an age-structured TB model and present some basic properties. In Sect. 3, we define the basic reproductive number and prove the local stability of the disease-free equilibrium and the unique endemic equilibrium. In Sect. 4, we present the uniform persistence result and prove the global stability of the disease-free equilibrium and the unique endemic equilibrium. In Sect. 5, we perform data fitting and sensitivity analysis of the basic reproductive number, and assess the feasibility of the WHO End TB Strategy by 2030. In Sect. 6, we give a brief discussion.

The model divides the total population at time t into five mutually-exclusive subgroups: susceptible class, latent class (those who are infected but not yet infectious), infectious class, treatment and recovered class denoted by S(t), e(t, a), I (t), T (t) and r (t, θ), respectively. Here the parameter a denotes the latent age of the infected individuals, θ denotes the relapse age of the recovered individuals. We assume that people during treatment may still be re-infected due to contact with TB patients. If some people have been treated successfully, but in the process of treatment, they re-contact with Fig. 1 Flowchart of the TB transmission tuberculosis patients, they will enter the latent class, and those who are re-infected in other situations will enter the infectious class. Burman et al. (2009) argued that recurrent TB can develop in some patients after they complete therapy. This shows that people who have recovered may have a relapse. The flow among those subgroups is shown in the following flowchart ( Fig. 1) .

Base on the above notations and the flowchart (Fig. 1) , we formulate the TB model as follows:

with the following boundary conditions

and initial conditions

where e 0 (a), r 0 (θ ) ∈ L 1 + (0, +∞), and s 0 , i 0 , t 0 ∈ R + . The meanings of the parameters in (1) are explained in Table 1 .

Throughout this paper, we make the following assumptions and notations. (A1) : , μ, β, ρ 1 , η, γ, ρ, α, μ I , μ t > 0; The recruitment rate of the susceptible class μ

The natural death rate of the population β

The transmission coefficient of TB ρ

Parameters that reduce the rate of transmission due to treatment δ(a)

The rate at which latent individuals progress into recovered class due to treatment σ (a)

The rate at which latent individuals progress into infectious class γ

The treatment rate of infectious class μ I The death rate due to TB μ t The death rate due to treatment α

The rate at which treatment individuals progress into recovered class due to complete treatment ρ 1

The rate of developing a new infections due to exposure to tuberculosis during treatment η

The proportion of individuals who become the latent class due to re-infection ω(θ)

The age-dependent relapse rate (A2) : ω(θ), δ(a), σ (a) ∈ L ∞ + (0, +∞) with essential upper boundsω,δ,σ > 0 , respectively;

For a, θ ≥ 0, we denote k 1 (a) = e − a 0 (μ+δ(s)+σ 1 (s))ds , k 2 (θ ) = e − θ 0 (μ+ω(s))ds , K 1 = +∞ 0 σ 1 (a)k 1 (a)da,

Base on above assumptions, we can verify the local existence of unique and nonnegative solution of the model (1) with the nonnegative initial conditions (see Webb 1985 and Iannelli (1994) ), thus we obtain the following proposition.

Proposition 1 Let x 0 ∈ Y , then there exists > 0 and neighborhood B 0 ⊂ Y with x 0 ∈ B 0 such that there exists a unique continuous function,

According to the comparison principle, we have

which yields

Boundedness is a direct consequence of nonnegativity of solution. Then we have the following proposition.

Proposition 2 Let x 0 ∈ Y , then there exists a unique continuous semiflow,

is the solution to (1) with (0, x 0 ) = x 0 , and (4),(5) are satisfied for t ∈ R + . The following set is positively invariant for system (1)

From Proposition 2 and inequality (4), we obtain the following proposition.

Proposition 3 (i) The solution of (1), (t, ·), is point dissipative, and attracts all points in Y ;

Integrating the equations for e(t, a), r (t, θ) in (1) along the characteristic lines, t − a=const., t − θ =cost., respectively, we deduce that

The continuous semiflow { (t, ·)} t≥0 is said to be asymptotically smooth, if each positively invariant bounded closed set is attracted by a nonempty compact set.

We will use the following two lemmas proposed by Smith and Thieme (2011) to prove the asymptotic smoothness of the semiflow.

x), and the following hold for any bounded closed set A ⊂ Y that is forward invariant under :

Since Y is an infinite dimensional space, space L 1 + (0, +∞) is a component of Y . For infinite dimensional space, we cannot deduce precompactness only from boundedness. Thus, to derive the precompactness of Y , we apply the following lemma:

Lemma 2 Let B be a bounded subset of L 1 + (0, +∞). Then B has compact closure if and only if the following conditions hold:

We are now ready to prove a result on the semiflow generated by system (1) is asymptotically smooth.

Proof We first decompose into the following two parts: U 1 , U 2 defined respectively by

Thus, lim t→+∞ diam U 2 (t, A ) = 0. In the following we will show that U 1 (t, A ) has compact closure for each t ≥ 0. From Proposition 3, we know that S(t), I (t), T (t) remain in the compact set [0, r ] for all t ≥ 0. Next, we will show that e(t, a) and r (t, θ) remain in a precompact subset of L 1 + (0, +∞) which is independent of x. According to

and (2), it is easy to show that 0 ≤ e(t, a) ≤ (β(1 + ρ)r + ηρ 1 )re −μa .

Therefore, the conditions (i),(ii) and (iv) of Lemma 2 are satisfied. Now, we only need to check the condition (iii) of Lemma 2.

Then

Thus, the condition (iii) of Lemma 2 holds, then we can get that e(t, a) satisfies the conditions of Lemma 2. In a similar way, r (t, θ) also satisfies the conditions of Lemma 2. Therefore, we obtain U 1 (t, A ) has compact closure for each t ≥ 0, Using Lemma 1, we know semiflow is asymptotically smooth. This completes the proof.

Combining Proposition 3 and is asymptotically smooth, as well as the existence theory of global attractors, the following result is immediate from Theorem 2.6 in Magal and Zhao (2005) and Theorem 2.4 in D' Agata et al. (2006) .

The semiflow has a global attractor A in Y , which attracts any bounded subset of Y .

System (1) always has a disease free equilibrium E 0 = ( μ , 0, 0, 0 L 1 (0,+∞) , 0 L 1 (0,+∞) ).

One can obtain the reproduction number from system (1)

R 0 represents the average number of secondary infected individuals produced by a primary infected individual. Next, we investigate the endemic equilibria of system (1). A endemic equilibrium (S * , I * , T * , e * (a), r * (θ )) of system (1) satisfies the following equations:

From (8) we can easily find that if R 0 > 1, system (1) has a endemic equilibrium E * = (S * , I * , T * , e * (a), r * (θ )), where

Summarizing the discussions above, we have the following theorem.

Theorem 3 System (1) always has the disease free equilibrium E 0 . Moreover, apart from E 0 , if R 0 > 1, system (1) exists a unique endemic equilibrium E * .

Now we consider the linearized system of (1) at an equilibrium

, then removing the bar, we obtain the following linearized system:

with the following boundary conditions

Let

with initial conditions

We obtain from the system (10) that

Combined with the third equation of system (10), and by simple calculation, we have

We obtain the characteristic equation of system (1) at an equilibrium E as follows:

Next we will analyze the local stability of the equilibria. Rigorous proof of local stability require more thorough spectral analysis, which be referred to in Liu et al. (2017) . Liu et al. (2017) formulated an age-structured model as an abstract nondensely defined Cauchy problem, and Lemma 3.4 in Liu et al. (2017) shows that point spectrum and spectrum are equal. Thus, the growth rate of solutions is given by the point spectrum, so we only need to study the eigenvalues of the characteristic equation of system (1).

Proof We first consider the local stability of the disease-free equilibrium E 0 . The characteristic equation corresponding to E 0 is

Hence, if R 0 < 1, all roots of f (λ) = 1 have negative real parts, then the disease-free equilibrium E 0 is locally stable if R 0 < 1.

Next, we consider the local stability of the endemic equilibrium E * . The characteristic equation corresponding to E * is

If R 0 > 1, the endemic equilibrium E * is locally stable. Otherwise, f (λ * ) = 1 has at least one root λ * = a 2 + ib 2 satisfying a 2 ≥ 0. We notice that

Hence, if R 0 > 1, all roots of f (λ) = 1 have negative real parts, then the endemic equilibrium E * is locally stable if R 0 > 1. This completes the proof.

In this subsection, our purpose is to show that system (1) is uniformly persistent when R 0 > 1. Define

Theorem 5 The sets M 0 and ∂ M 0 are forward invariant under the semiflow (t, ·). Also, the disease-free equilibrium E 0 of system (1) is globally asymptotically stable for the semiflow (t, ·) restricted to ∂ M 0 .

This implies the fact that (t, M 0 ) ⊂ M 0 , i.e. M 0 is forward invariant under the semiflow (t, ·).

Next, we will prove ∂ M 0 is forward invariant under the semiflow (t, ·). Let (0, x 0 ) ∈ ∂ M 0 , we consider the following system

Since

where

By the formulation (6), we havê

Substituting (15) into the first two equations of (14) yield

where

dθ.

Since

The system (16) can be written

Thus, (16) has a unique solutionÎ (t) ≡ 0,T (t) ≡ 0 for t ≥ 0.

From (14), (15), we know thatê(t, s) = 0,r (t, s) = 0 for 0 ≤ s ≤ t, thus,

Similarly, we have r (t, θ) L 1 + = 0. By using (13) we can obtain that I (t) = 0, T (t) = 0, e(t, a) L 1

Finally, we prove the disease-free equilibrium E 0 of system (1) is globally asymptotically stable for the semiflow (t, ·) restricted to ∂ M 0 . In ∂ M 0 , system (1) can be written as follows

Obviously, lim t→+∞ S(t) = μ . Hence, the equilibrium E 0 is globally asymptotically stable restricted to ∂ M 0 . The proof of Theorem 5 is complete.

Proof Since the disease-free equilibrium E 0 is globally asymptotically stable restricted to ∂ M 0 , applying Theorem 4.2 in Hale and Waltman (1989) , we only need to prove

By way of contradiction, we assume that there exists a x 0 ∈ M 0 such that x 0 ∈ W s (E 0 ). Then we can find a list of

Denote (t, x n ) = (S n (t), I n (t), T n (t), e n (t, ·), r n (t, ·)). Then for all t ≥ 0, we have

and (t, x n ) ⊂ M 0 . From the system (1), we know that there exists t 0 ≥ 0 such that I (t) > 0 for t ≥ t 0 or T (t) > 0 for t ≥ t 0 , so without loss of generality, we can take t 0 = 0 and I (0) > 0. Since R 0 > 1, we can choose sufficiently large n such that μ > 1 n and

Now we construct the following system

Using similar analysis as in Sect. 2, we can get existence, uniqueness and nonnegative of solution to system (18). By the comparison principle, we know

By use of Volterra formulation (6), we havê

Substituting (20) into the first two equations of (18) yield

We can claim that at least one ofÎ (t),T (t) is unbounded. Otherwise, we can use Laplace transform in the first two inequations of (21) yield

where

After a simple calculation, we obtain

then there exists ε > 0 such that

for all λ ∈ [0, ε). According to (23), we have L[Î ](λ) < 0 for all λ ∈ (0, ε). But this contradicts the nonnegative ofÎ (t)(t ≥ 0). Hence, at least one ofÎ (

is unbounded. This contradicts the Proposition 3. Therefore, W s (E 0 ) ∩ M 0 = ∅. By Theorem 4.2 in Hale and Waltman (1989) , we get that semiflow { (t, ·)} t≥0 generated by system (1) is uniformly persistent. By Theorem 3.7 in Magal and Zhao (2005) , we get that there exists a compact subset A 0 ⊂ M 0 which is a global attractor for { (t, ·)} t≥0 in M 0 .

Theorem 7 The disease-free equilibrium E 0 of system (1) is globally asymptotically stable if R 0 < 1.

is non-negative and continuous in (0, +∞) with a unique root at x = 1. Using similar arguments to the proof of Lemma 4.2 in Browne and Pilyugin (2013), we can get that any solution in A satisfies that S(t) > 0 for t ∈ R (In Theorem 2, we have proved that A exists). Next, we construct the following Lyapunov function L = L 0 + L 1 + L 2 + L 3 + L 4 on the global attractor A, by the compactness of A, we can easily deduce L is bounded on A, where

Calculating the derivative of L 0 ,L 1 , L 2 , L 3 , L 4 along solutions of system (1), respectively. We can deducė

Notice that if R 0 < 1, then d L dt ≤ 0 holds. According to the proof of the Theorem 4.3 in Guo et al. (2019) , we can easily obtain that A = {E 0 }. This proves that E 0 is globally asymptotically stable. The proof is complete.

If R 0 > 1, we have know that the system (1) is uniformly persistent and has a global attractor A 0 in M 0 . Let x ∈ A 0 , we can find a complete orbit { (t, ·)} t∈R through x in A 0 . By similar analytical method used in [ Browne and Pilyugin (2013) , subsection 3.2], system (1) can be written as

In order to study the global stability of E * , we first prove the following lemma.

Lemma 3 There exist , M > 0 such that all solutions in A 0 for t ∈ R, the following inequalities are satisfied

Next we prove that I (t) > 0, T (t) > 0 for t ∈ R. First we assume that there exists t 0 ∈ R such that I (t 0 ) = 0 and T (t 0 ) = 0, then according to I (t), T (t) equations in (24), we can deduce that I (t) = 0, T (t) = 0 for t ≤ t 0 . Next, form (24) we know On the other hand, we assume that there exists t 0 ∈ R such that one of I (t 0 ), T (t 0 ) is positive, so without loss of generality we assume that I (t 0 ) = 0, T (t 0 ) > 0. From

This contradicts the A 0 ⊂ M 0 . Hence, I (t) > 0, T (t) > 0 for all t ∈ R. In addition, from (24), we deduce that e(t, a) > 0, r (t, θ) > 0 for t ∈ R, a, θ ∈ R + . According to the compactness of A 0 , we know there exist , M > 0 such that

The proof is complete.

Theorem 8 If R 0 > 1, then endemic equilibrium E * of system (1) is globally asymptotically stable in M 0 .

Proof We define the following Lyapunov function V = V 1 + V 2 + V 3 + V 4 + V 5 on the global attractor A 0 (In Theorem 6, we have proved that A 0 exists), using Lemma 3, we can easily deduce V is bounded on A 0 , where

and A(a) = +∞ a (σ (u) + K 2 δ(u))e − u a (μ+σ (s)+δ(s))ds da,

Calculating the derivative of V along a solution in A 0 . When calculate the time derivative of V 1 , we notice that = β S * (I * + ρT * ) + μS * . 

we have

Thus, 0) ) − 1 + S I e * (0) S * I * e(t, 0) ) 0) ) − 1 + ST e * (0) S * T * e(t, 0) ) 0) ))) + (K 1 + K 2 K 3 )(ηρ 1 T * − ηρ 1 T e * (0) e(t, 0) ) (t, 0) ))) (t, 0) ))) (t, 0) ))) It follows from the non-negativity of g, we know that dV dt ≤ 0, and according to the proof of the Theorem 4.3 in Guo et al. (2019) , we can easily obtain that A 0 = {E * }. This proves that E * is globally asymptotically stable. The proof is complete.

In this section, we firstly use the system (1) to simulate the annual new TB cases data of China from 2007 to 2018, then we predict whether current tuberculosis control measures in China will be able to achieve the WHO End TB Strategy in 2030. To do so, we need to estimate the parameters of system (1). Some of the parameters in our model are distributed, which is different from ODE models, and we believe that the parameter values may be different due to different research groups and different time units. According to the WHO (2019) report in 2018, the TB incidence rate is falling, and the treatment success rate is also falling. Through the analysis of these data, it can be found that the parameter values in system (1) (Table 2 ) come from the Chinacdc (2019). Next, we estimate the unknown parameters and initial valueŝ θ = (β, δ, ρ, σ 1 , σ 2 , δ 1 , δ 2 , ω 1 , ω 2 , η, ρ 1 , μ I , μ T , α, e(0), I (0), T (0), r (0)) of system (1), where we assume δ(a) = δ 1 e −δ 2 a , σ (a) = σ 1 e −σ 2 a , ω(a) = ω 1 e −ω 2 a , e 0 (a) = e(0)δ 2 e −δ 2 a , r 0 (θ ) = r (0)ω 2 e −ω 2 θ . Let P(t,θ) be the number of new TB cases from model (1) at the t th year, then P(t,θ) can be written as where X (t) represents the cumulative number of people infected with TB by the t th year, and satisfies the following equation

Thus, we use P(t,θ) to simulate the new TB cases per year. We are attracted by the advantages of algorithms mentioned in Seyedali Mirjalili (2014) when reading it, and we find that the convergence was good and the error was small after combining with the model (1) (the mean absolute percentage error of TB cases is 0.75%). By Grey Wolf Optimizer (GWO) algorithm, we estimate the optimal parameters and initial values for model (1) ( Table 3) .

The comparison between the new TB cases in China from 2007 to 2018 and the simulation of P(t,θ) from model (1) is given in Fig. 2 . Moreover, we also find the basic reproduction numbers R 0 = 1.0608. Moualeu et al. (2015) , it can be found that the transmission rate of the Moualeu et al. (2015) is β 3 = 6.33563 × 10 −6 , while that of this paper is β = 1.1479 × 10 −10 . Analysis shows that the Moualeu et al. (2015) studies the transmission of TB in Cameroon from 1994 to 2010, and the new cases of TB in Cameroon increased rapidly during this period. In this paper, we used new TB cases in China from 2007 to 2018, during which the number of new TB cases in China decreased gradually. Therefore, the transmission rate of this paper is much lower than that of Moualeu et al. (2015) . This is also why some parameters in this paper need to be estimated. Parameter estimation is very important for studying the transmission of TB. Arregui et al. (2018) and Dowdy et al. (2013) also provide some methods to analyze the parameters in tuberculosis models.

The GWO algorithm mimics the leadership hierarchy and hunting mechanism of grey wolves in nature. Four types of grey wolves such as alpha, beta, delta, and omega are employed for simulating the leadership hierarchy. In addition, the three main steps of hunting, searching for prey, encircling prey, and attacking prey, are (1) implemented. Twenty nine test functions were employed in order to benchmark the performance of the proposed algorithm in terms of exploration, exploitation, local optima avoidance, and convergence. The goal of this paper is to minimize the target function

where r (t) represents the actual number of new TB cases at the t th year.

The outputs of system (1) is governed by the system input parameters and the initial values, but some of these parameters and initial values are obtained by data fitting, which may exhibit some uncertainty in their selection. The purpose of uncertainty analysis (UA) (Samsuzzoha et al. 2013; Ghosh et al. 2018; Marino et al. 2008) is to determine the reliability of parameter estimates. To ensure the reliability of the estimates of β, ρ 1 , γ, μ I , μ t , α, ρ, η obtained through the GWO, we employ Markov Chain Monte Carlo (MCMC) method with Delayed Rejection and Adaptive Metropolis (DRAM) algorithm (DRAM is an algorithm combining DR and AM in MCMC method, and the efficiency of the combination is demonstrated with various examples in Haario et al. (2006) ). The convergence of the chain is confirmed by using the Geweke's Z-scores (According to the MCMC program (MCMC 2019), the closer the Geweke's Z-score is to 1, the better the convergence of Markov Chain) . The mean, Year 5 The new TB cases 10 5 Fig. 3 The fitting results of the number of new TB cases reported from 2007 to 2018. The solid black line represents the fitted data, and the red dots represent the actual data. The areas from the darkest to the lightest correspond to the 50%, 90%, 95% and 99% posterior limits of the model uncertainty standard deviation and 95% confidence interval of the estimated parameters are shown in Table 4 , and the fitting result can be seen in Fig. 3 . Sensitivity analysis (SA) is to identify critical parameters that have significant impact on the basic reproduction numbers R 0 and to quantify how parameters uncertainty impact R 0 . Now, a global sensitivity analysis is usually implemented by using sampling-based methods. We will use partial rank correlation coefficient (PRCC) method to study SA. For the parameters in Table 3 , we fix σ 1 , σ 2 , δ 1 , δ 2 , ω 1 , ω 2 , e(0), I (0), T (0), r (0), other parameters are considered to obey normal distribution, their mean and standard deviation are shown in Table 4 . The sampling method is Latin hypercube sampling(LHS). We draw 1000 samples and obtain distribution histogram of the basic reproduction number R 0 (see Fig. 4 ). From Fig. 4 , we know that the distribution of the basic reproduction number in the range [0.6838, 1.5672], and the average is 1.0660. Combined with the analysis in Sect. 4, we can conclude that under current control measures, the possibility of extinction of TB is small. We calculate the PRCC between the parameters and the basic reproduction number R 0 (see Fig. 5 ).

From the PRCC values, we can know that β, μ t , ρ have the most important impact on R 0 , β represents the transmission rate coefficient of TB, and β is in direct proportion to bc, b represents the number of times an infected person has contact with other person per unit time, c represents the probability of infection per contact. ρ represents reduce the rate of transmission due to treatment. μ t represents the death rate due to treatment. 

We explain the specific steps of the LHS-PRCC scheme by studying the relationship between parameters parameters x 1 , x 2 , x 3 and output y = f (x 1 , x 2 , x 3 ).

Step 1 Take uniform distribution as an example, the parameter distributions are divided into N equal probability intervals, which are then sampled. N represents the sample size (see Fig. 6 ).

Step 2 We set up the LHS matrix by assembling the samples from each probability density functions. Each row of the LHS matrix represents a unique combination of parameter values sampled.

Step 3 Next, X and Y are rank transformed to X R and Y R in order from small to large.

Step 4 We calculate the residuals X 1 = X R1 −X R1 and Y 1 = Y R −Ŷ R , whereX R1 andŶ R are the following linear regression models:

Then

is the PRCC between x 1 and y. In the same way, we can calculate the PRCC between x 2 and y, and the PRCC between x 1 and y.

Over the past few decades, China has made great efforts to control the spread of TB, such as vaccination of newborns with BCG and active treatment of tuberculosis patients and so on. However, based on the above theoretical analysis and numerical simulation results, by using the current TB control measures, China may not reach the WHO's goal by 2030. To do so, China should give more feasible control measures.

In the above analysis, we know that β, μ t , ρ are the most important factors for TB control. Next, we adjust β, ρ to predict whether an optimal control measure can be given to achieve the WHO's goal by 2030. We use the parameter values in Table 2 as the baseline to compare the following control effects. First, we only consider changing the value of parameter β (see Fig. 7 ), we can find that the value of β needs to be reduced by 90% to reach the WHO 2030 target. Next, we only consider changing the value of parameter ρ (see Fig. 8 ), we can find the value of ρ needs to be reduced by 95% to reach the WHO 2030 target. Finally, we consider changing both parameter β and parameter ρ (see Fig. 9 ), we can find that if we can reduce parameter β by 70%, and reduce parameter ρ by 70%, then we will reach the WHO 2030 target. Parameters β and ρ represent the TB transmission coefficient and the reduction of infectiousness of treated individuals infected with TB, respectively. Some researchers mentioned that media coverage and public education campaigns can reduce the transmission coefficient β of disease (Barbara et al. 2020) . Barbara et al. (2020) suggests that the following measure can accelerate progress toward TB elimination: maintaining awareness of both the incidence of TB disease for the public through news releases and posting of information online. Take the COVID-19 pandemic as an example, media coverage and public education campaigns play very important role in the process of controlling the pandemic. People know more about the COVID-19 and enhance their self-protecting awareness by the media report about the COVID-19, and will change their behaviours and take correct precautions such as frequent hand-washing, wearing masks, reducing the party, keeping social distances, and even quarantining themselves at home to avoid contacting with others (CDC 2021). Therefore, if the media strengthen the publicity of tuberculosis and COVID-19, it will certainly have a good effect (Saunders and Evans 2020; The Lancet Infectious Diseases 2021). Reducing the contact between treated patients and susceptible individuals can reduce the parameter ρ.

Currently, the treatment of TB patients in China includes hospitalization and home treatment, and they are not forced to be hospitalized, so these people may still infect susceptible individuals. Based on the above analysis, China should strengthen the publicity of tuberculosis knowledge so that susceptible people can know how to avoid being infected. Also, patients who are treated with TB should be given longer hospital stays to reduce their chances of contact with susceptible people. By using these control measures, China may reach the WHO End TB Strategy in 2030. 

In this study, we used an age-structured mathematical model with treatment and relapse to study the transmission dynamics of TB. Sufficient conditions were derived for the global asymptotic stability of TB-free equilibrium and the endemic equilibrium. We estimated model parameters by fitting the annual new TB cases data of China, and concluded that, by using current tuberculosis control measures, China may not reach the WHO's goal by 2030. In order to achieve the WHO's goal, China needs to develop more practical control measures. PRCC values of the basic reproduction number, R 0 , with respect to some important model parameters demonstrate that the control measures for the spread of TB include reduction of the TB transmission coefficient β and reduction of the TB transmission coefficient βρ of treated individuals infected with TB. These can be achieved through media coverage, public education campaigns and increased hospital stays for TB patients.

TB is prevalent in many countries including China (Creswell et al. 2014) . The transmission mechanism of TB in different countries is the same, so the TB model in this paper can be used in the research of TB in other countries, and some feasible measures to control the transmission of TB can also be put forward by using local data of TB and some population data such as birth rate and death rate.

Optimal control of a two-strain tuberculosis-HIV/AIDS co-infection model

Optimal control strategies for dengue transmission in Pakistan

Data-driven model for the assessment of mycobacterium tuberculosis transmission in evolving demographic structures

Propagation phenomena for partially degenerate nonlocal dispersal models in time and space periodic habitats

Uniqueness and stability of time-periodic pyramidal fronts for a periodic competition-diffusion system

Essential components of a public health tuberculosis prevention, control, and elimination program: recommendations of the advisory council for the elimination of tuberculosis and the national tuberculosis controllers association

Global analysis of age-structured within-host virus model

Relapse associated with active disease caused by Beijing strain of mycobacterium tuberculosis

Environmental variability in a stochastic epidemic model

Optimal intervention strategy for prevention tuberculosis using a smokingtuberculosis model

Tuberculosis in BRICS: challenges and opportunities for leadership within the post-2015 agenda

Asymptotic behavior in nosocomial epidemic models with antibiotic resistance

Data needs for evidence-based decisions: a tuberculosis modeler's "wish list

Existence and asymptotic behaviors of traveling waves of a modified vector-disease model

A simple SI-type model for HIV/AIDS with media and self-imposed psychological fear

Global dynamics of an age-structured malaria model with prevention

Dram: efficient adaptive MCMC

Persistence in infinite-dimensional systems

Quasilinear elliptic equations with exponential nonlinearity and measure data

Marcinkiewicz estimates for solution to fractional elliptic Laplacian equation

Modelling and analysis of an alcoholism model with treatment and effect of twitter

Dynamics for an SIRS epidemic model with infection age and relapse on a scale-free network

Dynamics of an edge-based SEIR model for sexually transmitted diseases

Mathematical theory of age-structured population dynamics

The Basic Approach to Age-Structured Population dynamics

Modeling the effects of meteorological factors and unreported cases on seasonal influenza out breaks in Gansu province, China

Stability and bifurcation for a single-species model with delay weak kernel and constant rate harvesting

Global stability for a tuberculosis model

Oscillations in age-structured models of consumer-resource mutualisms

Predicting the cumulative number of cases for the COVID-19 epidemic in China from early data

Spatiotemporal dynamics of a diffusive Leslie-Gower prey-predator model with strong Alee effect

Global attractors and steady states for uniformly persistent dynamical systems

A methodology for performing global uncertainty and sensitivity analysis in systems biology

Bifurcation and control in a singular phytoplankton-zooplankton-fish model with nonlinear fish harvesting and taxation

Dynamics analysis of a predator-prey system with harvesting prey and disease in prey species

Parameter identification for a tuberculosis model in Cameroon

Epidemiological models of mycobacterium tuberculosis complex infections

Global stability of a multipatch disease epidemics model

Global stability in a tuberculosis model of imperfect treatment with age-dependent latency and relapse

Uncertainty and sensitivity analysis of the basic reproduction number of a vaccinated epidemic model of influenza

Covid-19, tuberculosis and poverty: preventing a perfect storm

Grey wolf optimizer

Klein-Gordon-Zakharov system in energy space: blow-up profile and subsonic limit

Providence The Lancet Infectious Diseases (2021) Tuberculosis and malaria in the age of covid-19

Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa

Optimal control and cost-effectiveness analysis of a Zika virus infection model with comprehensive interventions

Theory of nonlinear age-dependent population dynamics

A dynamic model for tuberculosis transmission and optimal treatment strategies in South Korea

Modeling the effects of health education and early therapy on tuberculosis transmission dynamics

Extremal solutions for nonlinear first-order impulsive integro-differential dynamic equations

Dynamic behavior of a stochastic SIQS epidemic model with Levy jumps

Dynamic behavior of a stochastic SIR epidemic model with vertical transmission

Dynamics of tuberculosis with fast and slow progression and media coverage

Analysis of transmission and control of tuberculosis in Mainland China, 2005-2016, based on the age-structure mathematical model

Traveling wave solutions of a diffusive SEIR epidemic model with nonlinear incidence rate

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

When calculate the time derivative of V 3 , we notice that 

Conflict of interest The authors declare that they have no conflict of interest.