Modeling and dynamic analysis of tuberculosis in mainland China from 1998 to 2017: the effect of DOTS strategy and further control

Background Tuberculosis (TB) is one of the most important health topics in the world. Directly observed treatment and short course chemotherapy (DOTS) strategy combines medicine care and modern health system firmly, and it has been carried out by World Health Organization (WHO) since 1997. In the struggle with TB, China has promoted the process of controlling the disease actively, and the full coverage of DOTS strategy has been reached around 2004. Mathematical modeling is a very useful tool to study the transmission of diseases. Understanding the impact of DOTS strategy on the control of TB is important for designing further prevention strategy. Methods We investigate the impact of control strategy on the transmission of TB in China by dynamic model. Then we discuss further control for TB aiming at developing new vaccine and improving treatment. The optimal control problem, minimizing the total number of infectious individuals with the lowest cost, is proposed and analyzed by Pontryagin’s maximum principle. Numerical simulations are provided to illustrate the theoretical results. Results Theoretical analysis for the epidemic model is given. Based on the data reported by National Bureau of Statistics of China (NBSC), the basic reproduction number of each stage is estimated and compared, and they are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {R}_{0}^{1}=1.7885$\end{document}R01=1.7885 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {R}_{0}^{2}=1.0741$\end{document}R02=1.0741, respectively. Optimal control strategy for further control is designed and proved well. An intuitionistic comparison between the optimal control strategy and the current control strategy is given. Conclusions The diagnosis and treatment of TB in China have been promoted a lot and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {R}_{0}$\end{document}R0 is reduced by the full coverage of DOTS strategy. However, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {R}_{0}$\end{document}R0 in China is still greater than 1 now. The relationship between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {R}_{0}$\end{document}R0 and vaccination strategy is shown. Optimal strategy aiming at exposed and infected population is suggested for further control.


Background
Tuberculosis caused by Mycobacterium tuberculosis (Mtb) is a widespread and deadly chronic infectious disease. Generally speaking, the Mtb usually affects the lungs. The disease is transmitted through droplets from the throat and lungs of people with acute infection, and about one-third of the world's population has been infected with Mtb nowadays [1]. The Mtb carriers may become infectious individuals all their lives. The probability of changing from the latent state to active state is only about 1/10 to 1/5 . However, in the statistics of 2017 [2], about 10.0 million of the people get TB, of which 59% are new cases. TB is still one of the top 10 causes of death. The prevention, diagnosis and treatment of TB need adequate funding. Because TB is a chronic disease, the material support must be sustained for a long time. In the recent decade, the finance support for the prevention of TB has increased, but funding gaps still exist. The eradication of TB is a worldwide scale challenge especially in South-East Asia and Africa. Chemotherapy and antibiotics have become powerful measures for the elimination of disease since 20th century. Most deaths of TB can be prevented with early diagnosis and effective treatment which had saved about 54 million of people between 2000 and 2017 [2]. DOTS strategy is a control measure formulated and carried out by WHO to end the transmission of TB. A lower proportion of recent transmissions is observed in the county with long-term DOTS implementation [3] and DOTS strategy is now considered as the most suitable control measure for the cost-effective principle. China is one of the highest TB burden countries in the world, more than 80% infectious are in the mid-western regions of China and the drug-resistance problem is very serious. For the prevention and control of TB, Chinese government has taken a series of measures including the universal coverage of DOTS strategy. The implement of this strategy has made prominent improvement on diagnosis and treatment for TB.
In the past several decades, mathematical model is playing an important role in disease control. The dynamic system can describe the transmission process of TB. Differential equation models are widely used in the study of TB, particularly, complex network models are applied to describe the complexity because of topological characterization. Five different complex network models are established to study the transmission of TB and to help to eradicate the disease [4]. Data fitting gives a proper combination of theoretical analysis and practical situation. Demographic and epidemiological data have been used to explore the effects of population growth, randomness, clustering of contacts and age structure on TB dynamics in [5]. Now there are several works focused on TB, which have been studied through many factors, such as fast and slow progression [6], drug-resistant strains [7], reinfection [8,9], co-infection [10], migration and seasonality [11,12]. For the prevention of TB, the model formulated with isolation [13], treatment [14], vaccination [15] or a combination of different control strategies [16][17][18] are discussed. On the other hand, if the people take initiative on the awareness of prevention and control, the effect will be better. Das et al. [19] investigate the influence of social awareness spread by media in TB transmission dynamics and the corresponding optimal control strategy with minimum cost is also given. To make a better control, optimal control theory is applied in many areas, such as the design of therapy [20], optimal control of the disease among animals [21] and best prevention strategy for TB [22,23]. Nowadays, the impact of DOTS strategy on the control of tuberculosis in China is scarce, moreover, the study of the whole situation of tuberculosis in China is still absent since the reported data is available.
The new reported cases of tuberculosis in China increases sharply around the year of 2004, so the study is divided into two stages. In this paper, we study the TB trends mainly focused on the full coverage of DOTS strategy in China in recent two decades. Theoretical analysis for the dynamic system and the estimation of the basic reproduction number for before and after the 100% coverage of DOTS strategy which is based on the reported TB cases in China are given. We study the efforts of DOTS strategy through the basic reproduction number of the two stages. The full coverage of DOTS strategy has reduced the level of TB obviously, but the disease still exists (R 0 > 1). Then various further controls for TB are analyzed and simulated. Our theoretical results and numerical simulations do great help to design optimal control strategy in practice.

Mathematical model
In this section, we formulate a dynamic model to investigate the current situation of TB transmission in China. In this model, the total size of population at time t is divided into five subclasses: susceptible, vaccinated, exposed, infectious and recovered denoted by S, V , E, I and R, respectively. Because the vaccination strategy for TB in China focuses on newborns, the subclass V is regarded as the newborns vaccinated successfully. During the protection period of BCG, people cannot get infected, even if they contact with infectious individuals. But the time-efficiency of BCG is very limited, the vaccinated individuals may lose immunity and become susceptible. The susceptible class is increased at recruitment rate and at the rate of k which is from vaccinated class by losing temporary immunity. The parameter p ∈[ 0, 1] indicates the fraction of the newborns vaccinated successfully. G(β) is the transmission coefficient. In the first stage, i.e., the time from 1998 to 2003, we take G(β) = β; in the second stage, i.e., the year from 2004 to 2017, G(β) is taken as ωβ in which ω measures the increasing of the reported rate due to the improvement of diagnosis. Bilinear incidence G(β)SI is used to represent the rate at which the susceptible individuals are infected. d is the natural death rate coefficient, σ is the disease-induced death rate, ε is the rate of progression to active TB, and γ 1 and γ 2 denote the treatment rate for latent class and infectious class, respectively.
The above model description can be constructed as follows: Let N denote the total population size Then, it follows from (1) that As t → ∞ in (3), the total population size N → /d. It implies that the positively invariant set for system (1) is The dynamics behaviors of system (1) on the set is sufficient to study in the following parts.

Dynamics analysis
There is a unique infection-free equilibrium of system (1): Then by the principle of the next generation matrix in [24], we can obtain the basic reproduction number where ρ represents the spectral radius of matrix FV −1 and the matrices F and V are given by ε/(ε+d+γ 1 ) is the fraction of latent population becoming infected among those exiting the E class and 1/(d +δ +γ 2 ) represents the mean duration of exiting the I class. We have the following results, and the relevant proof can be found in Appendix.

Estimation and analysis of R 0
By the above analysis, we know that the basic reproduction number is an important threshold for the control of disease. In order to estimate the value of basic reproduction number, we use the data of annual reported TB cases which have been released by the National Data [25] of China to do data fitting. In this paper, we study the reported data form 1998 to 2017. Notice that, during the period of 2004 to 2005, China has completed the 100% coverage of DOTS strategy for TB, and the diagnosis and treatment have been improved a lot. We divides the time into two stages (before and after the full coverage of DOTS strategy) to study the transmission and control of TB in China. The parameters are estimated with the reported data by least square method (LSM).
For the first stage, the time period is from 1998 to 2003 in which the coverage rate of DOTS in China is about 60%. By fitting system (1) to the annual reported TB cases, we obtain the estimations for the transmission rate β and the recovery rate of infectious γ 2 , which are listed in Table 1 and the goodness of fit is shown in Fig. 1, stars represent the newly reported TB cases by year, the red solid line shows the fit based on the true data and the dashed line shows the estimation of new reported cases in 2004 and 2005. In our estimate, the new reported cases in 2004 should be 7.0 × 10 5 , but actually the collected data in 2004 is 9.7 × 10 5 . It implies that the diagnosis of tuberculosis has been improved a lot. Using the estimated parameter values, we calculate the basic reproduction number R 1 0 of the first stage is 1.7885.  [12,28] ε Rate of progression to infec-6 year −1 [29] tious stage from the exposed For the second stage, we fit system (1) to the annual reported TB cases from 2004 to 2017. The DOTS coverage rate increased rapidly to 95% in 2004 and achieved 100% coverage since 2005. Compared with the first stage, not only for the control of TB, the integral medical level of China has promoted a lot, one of the most obvious reactions is the average lifetime has risen to about 72 years. To make a better data fit, we take the natural death rate as 0.0139year −1 in the second stage. With LSM, the transmission rate and the recovery rate of infectious are  Table 1. The goodness of fit is shown in Fig. 2 with red solid line and the black dashed line shows the trend of new reported cases under current control strategy. According to the parameter values we estimated, the basic reproduction number of the second stage is R 2 0 = 1.0741. Judging from the data we collected, the newly reported TB cases increase at a relatively stable speed in the first stage, but proliferate fast in the prophase of the second stage and then decrease year by year. It is closely related to the full coverage of DOTS strategy. In the prophase of the 100% coverage of DOTS strategy, the diagnosis has developed rapidly, the data we have studied comes from the annual summary of the number of TB patients in different levels of hospitals. Due to the improvement of diagnosis and the strengthening of screening, the newly reported TB cases increase first after the full coverage of DOTS strategy. It is reflected in the estimation of parameter G(β) in our model. The parameter G(β) in the second stage is larger than in the first stage. And the corresponding treatment measures have been improved a lot, it makes the newly reported TB cases decreased since 2009. It can be seen from the significantly increase recovery rate of the infectious (γ 2 ). A more intuitionistic expression is that the period in which the infected individuals need to recover (1/γ 2 ) is shorten sharply.
Comparing the basic reproduction number of the two stages, we can clearly see that under the full coverage of DOTS strategy, the control of TB has improved a lot and it is strongly related to the promotion of diagnosis and treatment. The decrease of the basic reproduction number indicates that the level of transmission for TB is declining in population. Although the newly reported TB cases in China have fallen in recent years and it seems that the disease is gradually vanished, however, the basic reproduction number is still greater than 1, the elimination of TB cannot achieve in this situation. The estimation of R 0 is very much in line with previous studies [32][33][34]. In order to eradicate the disease, we need to take a further control.

Further control for TB
According to the above results, although the newly reported TB cases in China have decreased since 2009 and after the full coverage of DOTS strategy the basic reproduction number has reduced obviously, the basic reproduction number is still greater than 1. Further control for TB must focus on the development of new vaccine and new method for diagnosis and treatment. First, we consider that if the new vaccine comes into use mainly from two sides: increase the vaccination rate and protection period.
To explore the effect of successful vaccination (p) on the basic reproduction number R 0 , we calculate the derivative of R 0 with respect to the parameter p as follows Because all the parameters are non-negative, ∂R 0 /∂p < 0. Figure 3a illustrates how the basic reproduction number  Table 1. Therefore, increasing the fraction of successful vaccinated infants helps to reduce TB. Similarly, to investigate the effect of extending the protection period of vaccination (1/k) on the basic reproduction number R 0 , direct calculation shows that Figure 3b shows how the basic reproduction number R 0 depends on the various weights of k with different p: as k decreases, the basic reproduction number decreases.
Since ∂R 0 /∂k > 0, increasing the protection period (1/k) can do active effect on TB control. All the other parameters are shown in Table 1. Now, we consider how to give an optimal control for TB through the development of treatment. We mainly discuss on the estimation of parameter G(β) and γ 2 in the two stages and the corresponding R 0 . It is related to the diagnosis and treatment measures nowadays. Although the full coverage of DOTS strategy strengthens the isolation and treatment of the infectious individuals, the treatment rate for latent individuals still stays in a relatively low level. So in our previous study, we fix the recovery rate of exposed (γ 1 ) as the same value through the two stages, and for a further control of TB, there is still a big space for improvement. Then we consider both γ 1 and γ 2 as control functions in the following strategy. An optimal control model for TB is formulated by extending the system above to include control functions. For continuous system, Pontryagin's maximum principle is an important tool to solve optimal control problems. The objective function usually reflects the bioscience goal, and it has been used to describe optimal therapy and treatment strategy [35,36]. In our model, the main aim is to minimize the total number of infected individuals over a finite time interval with the lowest cost possible. The dynamical control system is It is necessary to do early intervention on the exposed population. In system (8), two control functions u 1 (t) and u 2 (t), represent the treatment rates for latent class and infectious class, respectively. We consider state system (8) of ODEs in R 5 with the set of admissible control functions given by and a i , b i (i = 1, 2) are fixed positive constants.
The objective functional is defined as Our target is to minimize the infected individuals (including latent and infectious individuals) as well as the costs required to control TB by treating latent and infectious individuals. The constants A i (i = 1, 2) are the weight constants for class E and I respectively, B i (i = 1, 2) are the measures of the relevant cost of the interventions associated with u 1 (t) and u 2 (t).
An optimal control (u * The Hamiltonian H associated with (8), (9) and (10) is given by where λ i (i = 1, · · · , 5) are adjoint variables, which can be determined by solving the following system of differential equations: satisfying the transversality condition By calculating, we have And the characterized control set of (u * 1 , u * 2 ) is which is determined by

Numerical simulations
In the foundation of the analysis and research above, we study the optimal control problem for TB. In this section, we perform numerical simulations of the optimal system. To obtain the optimal solutions, we use an iterative method with the help of Runge-Kutta fourth order procedure. The detail algorithm is as follows: Step 1: We solve system (8) with an initial guess for the control (u 1 (t), u 2 (t)) over the time interval [ 0, t f ] using a forward fourth order Runge-Kutta scheme.
Step 2: With the application of the current iteration solution of the state variables (8) and the transversality conditions λ i (t f ) = 0, i = 1, · · · , 5, the backward Runge-Kutta fourth order procedure is used to solve adjoint system (14) in the same time interval. The controls are updated by combining the state and adjoint values.
Step 3: Then repeat the iterative process. This situation continues until the values of the unknowns at the previous iteration are very close to the ones at the present iteration.
The time interval for which the optimal control is taken as 20 years here. Based on the goodness of fit for the second stage in "Results" section, we forecast the value of S, V , E, I and R in 2020 and take them as the initial values of our optimal control which are given by S(0) = 965740160, V (0) = 36381737, E(0) = 115575, I(0) = 1575311 and R(0) = 289842548 respectively. Here, we suppose that the vaccine has developed, the protection period of the vaccine is taken as 20 years and the fraction of vaccination is raised to 90%. Other parameters follow the second stage in Table 1. For the latent population, because of its asymptomatic, the diagnose of the latent individuals is harder than the infectious individuals and the cost of treatment in latent population is much lower than the infectious population. To measure the difference between latent and infectious population, the weights in the objective function are A 1 = 1, A 2 = 10 3 , B 1 = 10 3 and B 2 = 10 6 [16]. The bounds on the control are based on the effectiveness of the intervention strategy. It is reasonable to assume that the upper bound for u 1 (t) is 1 and lower bound is 0; the upper bound for u 2 (t) is 0.75 and lower bound is 0.
We compare our optimal control strategy with the results of system (1) under current control strategy which is the situation of the second stage for the following 20 years. In Fig. 4, the population with optimal control actions are shown in blue solid lines and the red dashed lines represent the population under current control strategy. From Fig. 4 we can clearly see that the optimal control strategy gives a better result, that is in the case of optimal control, the decrease is more obvious in the exposed and infectious individuals than in the case of current control strategy.

Conclusions and discussion
In this paper, we propose mathematical models and use numerical simulations to describe the transmission of TB in China. The aim of system (1) is to reconstruct the past two decades epidemic situation of TB and to forecast the further spread of TB. A compartmental nonlinear deterministic epidemic model is formulated. Theoretic analysis for its global dynamical behaviors is given. R 0 < 1 ensures the global stability of the infection-free equilibrium. We examine the model with the reported data from 1998 to 2017 in China. The new reported cases of TB change dramatically around the year of 2004 due to the universal coverage of DOTS strategy, so we study the control of TB in China from two different stages, i.e., before and after the full coverage of DOTS strategy. The diagnosis and treatment a b c d Fig. 4 Comparison between the current TB control strategy and optimal control strategy level are the most important considered factors in our two stage-study. The basic reproduction number for each stage is 1.7885 and 1.0741 respectively. Although the DOTS strategy plays an important role to eliminate TB in China, which is mainly reflected in the basic reproduction number of the two stages, the disease still exists.
Because vaccines and treatment measures are developing fast nowadays, we study further control for TB from the two sides. R 0 is reduced by improving the fraction of vaccination and increasing the protection period of vaccine (see Fig. 3). On the other hand, to improve the treatment for early exposed and infected population is also explored. An optimal control problem to minimize the total number of infected individuals with the lowest cost is proposed and analyzed by Pontryagin's maximum principle. Numerical simulations of the optimal control system provide valuable results to illustrate the theoretical analysis. The optimal strategy is compared with the current strategy in Fig. 4. The discussion of further control for TB indicates that vaccination and treatment strategies are effective in reducing the transmission of TB. Our study provides theoretical support for making further control of TB and help the design of optimal control strategy in practice.