Dynamics and optimal control of tuberculosis model with the combined effects of vaccination, treatment and contaminated environments

1.   Introduction

Tuberculosis (TB) is a chronic infectious disease caused by Mycobacterium tuberculosis (MTB) infection ^[1]; usually, Mycobacterium tuberculosis can be transmitted through the air, such as in droplets. Tuberculosis usually invades various organs of the body, but mainly invades the lungs, so it is called pulmonary tuberculosis ^[2]. Before the COVID-19 outbreak in 2019, TB was the leading cause of death from single-source infections and the 13th leading cause of death in the world. Despite significant progress in preventing and mitigating TB over the past two decades, there were 10.6 million new TB cases worldwide in 2021. Among them, the number of new TB cases in China was 639,548, accounting for 6% of the global incidence, and the number of cases ranked top five in the number of Class B infectious diseases in China. According to the World Health Organization (WHO), there are nearly 2 billion people with latent TB infections worldwide ^[3]. Therefore, the research of tuberculosis prevention and control is urgent.

In recent years, many scholars have developed mathematical models of tuberculosis. Waaler et al. ^[4] introduced the first mathematical model of tuberculosis. Dye et al. ^[5] and Song et al. ^[6] proposed a model showing fast and slow progress starting with two potential classes. Porco and Blower ^[7] proposed a mathematical model that takes into account the recurrence of TB. Bhunu et al. ^[8] first proposed a model for vaccinating susceptible people with the Bacillus Calmette - Guerin (BCG) vaccine and treating infected people. Beniers and Gomes ^[9] showed that the particulars of the relation between environmental mycobacteria (EM) abundance and vaccine effectiveness depend on the degree of protection conferred respectively by natural infection, vaccination, and EM. Yang et al. ^[10] formulated a mathematical model to explore the impact of vaccination and treatment on the transmission dynamics of tuberculosis. Ziv et al. ^[11] used a mathematical model to predict the impact of TB vaccines on disease prevention in countries with high incidence. As time has passed, more factors have been included in the study of tuberculosis transmission dynamics. Ren ^[12] considered the treatment and incomplete treatment of tuberculosis. Trauer et al. ^[13] added drug resistance to the model and analyzed tuberculosis in highly endemic regions in the Asia-Pacific region. Bhunu ed al. ^[14] considered the co infected population of HIV and TB. Guo et al. ^[15] formulated an age-structured TB model to study the effects of relapse and treatment on transmission dynamics of TB. Xue et al. ^[16] proposed a non-autonomous differential equation model with age structure and seasonal transmission rate. Bowong and Kurths ^[17] included the seasonality of tuberculosis in the study.

Vaccination is an effective way to prevent and control diseases. Li and Guo ^[18] studied the spread of a mutated COVID-19 virus in China at a certain percentage of vaccinations, and then modeled how the Omicron strain developed in competitive transmission with the Delta strain ^[19]. The vaccine to prevent tuberculosis is a BCG vaccine, which is a live vaccine made from the suspension of attenuated tuberculosis bacillus bovis. It can enhance the activity of macrophages, strengthen the ability of macrophages to kill tumor cells, activate T lymphocytes, and enhance the cellular immunity function of the body. It was first developed by French scientists Calmette and Guérin ^[20,21]. It has long been thought that vaccination with the BCG vaccine does not protect against acquiring mycobacterium tuberculosis infection, but does offer some protection against progression of the active disease ^[22]. In 2022, Li et al. ^[23] took vaccination as a single compartment and proposed a $SVEITR$ model with six compartments, including vaccination, fast and slow progression, incomplete treatment, and relapse. In ^[23], they also analyzed the epidemic law of tuberculosis in the United States and proposed the optimal control measures and cost-effectiveness.

As is well known, treating infected patients is one of the most effective ways to eliminate diseases. At present, the main treatment method for tuberculosis is drug treatment, which generally follows the five principles of "early treatment, appropriate drug use, combined treatment, regular treatment, and full course treatment". However, due to various limitations, the treatment of patients often fails to achieve the expected results. For example, the current level of screening and treatment for latent individuals in China is far from sufficient, to the extent that a large number of latent individuals develop into active infections. Due to incomplete treatment, some infected individuals develop into drug-resistant infections. Therefore, Jiang et al. ^[24] studied the optimal control of tuberculosis by using the screening of latent individuals as a control strategy. Xu et al. ^[25] divided tuberculosis patients into drug sensitive persons (DS) and drug resistant persons (DR) for research.

Studies have shown that tuberculosis bacteria can survive outside the host for a long time under the right conditions, and even 75% alcohol cannot eliminate them. There is also evidence that susceptible individuals can become infected through contact with pathogens that live freely in the environment. Thus, in recent years, the effects of contaminated environments (towels, handkerchief, toys, and doorknobs touched by infected people) have gradually been incorporated into modeling studies of TB. For example, Cai et al. ^[26] took the contaminated environment as a separate chamber to study tuberculosis in the contaminated environment in Jiangsu Province, China.

In conclusion, we can see that vaccination, treatment, and contaminated environments have a strong impact on the development trend of TB. However, most existing articles have considered a single influencing factor, and so far no studies have considered the combined effects of vaccination, treatment, and contaminated environments on TB. Therefore, with reference to the transmission mode of TB in existing models and based on Chinese TB data, this paper proposes a new model that takes into account vaccination, treatment, and contaminated environments, which extends the existing mathematical model. To analyze the effects of BCG vaccination, cleaning up contaminated environments and treatment on tuberculosis in China, we classified vaccination, treatment, and contaminated environment as different compartments, and used the Markov chain Monte Carlo (MCMC) method to estimate model parameters, and Latin Hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) to study the uncertainty and sensitivity of model parameters.

The organization of this paper is as follows: In Section 2, we present a new mathematical model of tuberculosis. In Section 3, we analyze the dynamic behavior of the model. In Section 4, we introduce some control strategies, and discuss the optimal control strategy. In Section 5, we take tuberculosis in China as an example, use this model to fit the real data, and conduct a cost-benefit analysis. In Section 6, we discuss and summarize our conclusions.

2.   Model formulation

In this section, we take vaccinated individuals, treated individuals, and contaminated environment as compartments $V$, $T$, and $W$, respectively, and establish a $SVEITRW$ model with seven compartments. In this model, the population is divided into six compartments: $S(t)$, $V(t)$, $E(t)$, $I(t)$, $T(t)$, and $R(t)$, where $S$ represents susceptible individuals, $V$ represents vaccinated individuals, $E$ represents exposed individuals, $I$ represents infected individuals, $T$ represents treated individuals, $R$ represents recovered individuals, and $W$ represents the density of pathogens in the contaminated environment (towels, handkerchieves, toys, bedding, etc.) at time $t$, assuming the population size at time $t$ is $N(t) = S(t)+V(t)+E(t)+I(t)+T(t)+R(t)$. It should be noted that although the BCG vaccine does not completely prevent mycobacterium tuberculosis infection, it does provide some protection for the vaccinated population. Therefore, in this paper we assume that inoculated individuals transition to the exposed state after infection. The structure of the infectious disease chamber of this model is shown in Figure 1. The biological significance of the parameters is shown in Table 1. The specific model is as follows:

where $k$ and $h$ are non-negative values, while other parameters have positive values.

3.   Stability analysis and persistence

Next, we discuss the dynamical behavior of system (2.1). First, we introduce some of the symbols that will be used in this article. Let $\mathbb{R}^{n}$ denote the $n$-dimensional Euclidean space, and

For a continuous and bounded function $f: [0, +\infty)\rightarrow R$, denote $f_{\infty} = \lim\limits_{t\to +\infty}\inf f(t)$, $f^{\infty} = \lim\limits_{t\to +\infty}\sup f(t)$. Next, we assume that

is the initial condition of system (2.1). Clearly, any solution of system (2.1) with non-negative initial conditions is non-negative, and the following lemma shows that the solutions are uniformly ultimately bounded.

3.1.   Positively invariant set

Lemma 3.1. The solutions of system (2.1) are uniformly and ultimately bounded, that is, the solutions of system (2.1) eventually enter

Proof. It is easy to get that $\mathbb{R}^{7}_{+}$ is an invariant set for system $(2.1)$, that is, any solution initiating from $\mathbb{R}^{7}_{+}$ remains in it. Since the total population $N(t) = S(t)+ V(t)+ E(t)+ I(t)+ T(t)+ R(t)\ge0$, and $W(t)\ge0$, then from the system $(2.1)$ we have

According to the comparison theorem, there exists $t_{1} > 0$ such that $N(t)\le\frac{\Lambda}{d}$, for $t\ge t_{1}$.

It follows from the last equation of (2.1) that for $t\ge t_{1}$,

Similarly, there exists $t_{2}\ge t_{1}$ such that $W(t)\le {\frac{\omega\Lambda}{dm}}$, for $t\ge t_{2}$. Therefore, the solutions of system (2.1) are uniformly ultimately bounded. This completes the proof.

□

3.2.   The basic reproduction number

Obviously, system (2.1) has a disease-free equilibrium point

The next-generation matrix method ^[27] is a common method to calculate the basic reproduction number $R_{0}$. We define $X = (E, I, T, R, W)$, and system (2.1) is written as $\dot X = \mathcal{F}-\mathcal{V}$, where

and

Jacobian matrices $F$ and $V$ of $\mathcal{F}$ and $\mathcal{V}$ at $E_{0}$ in (3.2) and (3.3) are obtained respectively, where

The basic reproduction number $R_{0}$ is the spectral radius of $FV^{-1}$, namely, $R_{0} = \rho(FV^{-1})$. Therefore,

where

and $\left(k_{1} = d+v+\gamma_{1}, k_{2} = d+\mu+h+\gamma_{2}, k_{3} = d+\gamma_{3}+\mu\right.$, $\left.k_{4} = \eta+d, B = k_{2}k_{3}k_{4}-\eta(h\gamma_{3}+k_{3}\gamma_{2}) > 0\right)$, where $R_{1}$ is the average number of secondary infections produced by infected individuals, and $R_{2}$ is the average number of secondary infected individuals generated by free-living viruses in the contaminated environment shed by the infected individuals.

3.3.   Stability of disease-free equilibrium

In this section, we discuss the stability of the disease-free equilibrium for system (2.1).

Theorem 3.1. The disease-free equilibrium $E_{0}$ is locally asymptotically stable if $R_{0} < 1$ and unstable if $R_{0} > 1$.

Proof. The Jacobian matrix $J(E_{0})$ of system (2.1) at $E_{0}$ can be expressed as

where $X_{1} = (1-p)\beta_{1}S^{0}+\beta_{2}V^{0}$, $X_{2} = p\beta_{1}S^{0}$, $X_{3} = (1-q)\sigma_{1}S^{0}+\sigma_{2}V^{0}$, $X_{4} = q\sigma_{1}S^{0}$, and we can obtain the eigenvalues for (3.5) as $-(k+d)$, $-d$ and the roots of

From (3.6), we have

Next, we prove that all eigenvalues of equation (3.7) have negative real parts. For this reason, similar to the proof of Theorem 2 in Jiang et al. ^[24], we assume that there is at least one eigenvalue $\bar{\lambda}$ such that $\rm{Re}\bar{\lambda} \ge0$. From (3.7), we can conclude that

Therefore, we have

From (3.8), we have

Next, we will divide both sides of the inequality in equation (3.9) by $k_{1}k_{2}k_{3}k_{4}-k_{1}\eta(h\gamma_{3}+k_{3}\gamma_{2})$, and we have

This contradicts $R_{0} < 1$, which means our assumption is not valid. Therefore, in (3.5) all eigenvalues have negative real parts. This completes the proof.

□

Theorem 3.2. If $R_{0}\le1$, then the disease-free equilibrium $E_{0}$ is globally asymptotically stable.

Proof. Construct the following Lyapunov function:

where

The derivative of $H$ along the solutions of model (2.1), together with $S\le S^{0}$ and $V\le V^{0}$ from the positive invariance $\Omega$, yields

Since all the parameters of model (2.1) are positive, we directly derive that $\dot H\le 0$ for $R_{0} < 1$. It follows that there exists a singleton $E^{0}$, as the maximal compact invariant set in $\big\{(S, V, E, I, T, R, W)\in\Omega: \dot H = 0\big\}$. Using LaSalle Invariance Principle ^[28], every solution of $(2.1)$ with initial condition in $\Omega$ approaches $E^{0}$ as $t\rightarrow +\infty$ whenever $R_{0}\le1$.

□

3.4.   Persistence

In this following, we discuss the persistence for system (2.1). To prove our conclusion, we first introduce the theory of uniform persistence ^[29].

Lemma 3.2. Assume that

[1] $f(X_{0})\subset X_{0}$ and f has a global attractor $A$;

[2] The maximal compact invariant set $A_{\partial} = A\cap M_{\partial}$ of $f$ in $\partial X_{0}$, possibly empty, admits a Morse decomposition $\left\{M_{1}, \cdots, M_{k}\right\}$ with the following properties:

(a) $M_{i}$ is isolated in $X$.

(b) $W^{s}(M_{i})\cap X_{0} = \varnothing$ for each $i\le l\le k$.

Then, there exists $\delta > 0$ such that for any compact internally chain transitive set $L$ with $L\nsubseteq M_{i}$ for all $1\le i\le k$, we have $inf_{x\in L}{\rm d}(x, \partial X_{0} > \delta)$, that is to say $f:X\rightarrow X$ is uniformly persistent with respect to $(X_{0}, \partial X_{0})$.

Theorem 3.3. If $R_{0} > 1$, then system (2.1) is uniformly persistent.

Proof. We prove the uniform persistence of system (2.1), which is equivalent to proving that if $R_{0} > 1$, there exists a small positive constant $\varepsilon > 0$ such that $I_{\infty} > \varepsilon$ for system (2.1) with initial value $\varphi(0)$ and $I(0)$. Using the persistence theory described above, we define

Next we prove that system (2.1) is uniformly persistent with respect to $(X, X_{0})$.

First, it can be seen that both $X$ and $X_{0}$ are positively invariant for system $(2.1)$, and $X_{0}$ is relatively closed in $X$. Moreover, $\partial X_{0}$ is a relatively closed set in $X$. It follows from Lemma 1 that solutions of system (2.1) are uniformly and ultimately bounded. Thus, system (2.1) is point dissipative, for system (2.1) there must exist a global attractor. Set

We first prove that

As it is obvious that $M^{'}_{\partial}\subseteq M_{\partial}$, we then only need to prove $M_{\partial}\subseteq M^{'}_{\partial}$. Using proof by contradiction, we make the assumption that $M_{\partial}\subseteq M^{'}_{\partial}$ is not true. Let $\varphi(t)$ be a solution of system (2.1) with initial condition $\varphi(0)$. Then, for any

and $\varphi(t)\neq M^{'}_{\partial}$, there must exist at least one of $E(t)$, $I(t)$, $T(t)$, $R(t)$, and $W(t)$ which is not zero. Without loss of generality, we assume that $E(t) = 0$, $I(t) = 0$, $T(t) = 0$, and $R(t) = 0$, but $W(t) > 0$. From system (2.1), one has the following equations $\forall t > 0$:

This means that $\varphi(t)\not\in \partial X_{0}$ for $t > 0$, or, in other words, $\varphi(t)\not\in M_{\partial}$, which contradicts the assumption that $\varphi(t)\in M_{\partial}$. Hence, we obtain $M_{\partial}\subseteq M^{'}_{\partial}$. Since $M_{\partial} = M^{'}_{\partial}$, we obtain that $M_{\partial}$ only has the disease-free equilibrium $E_{0} = (S^{0}, V^{0}, 0, 0, 0, 0, 0)$ and $E_{0}$ is a compact and isolate invariant for $\varphi(0)\in M_{\partial}$.

Next, we only need to prove that $W^{s}(E_{0})\cap X_{0} = \varnothing$, where $W^{s}(E_{0})$ denotes the stable manifold of $E_{0}$. That is, there exists a positive constant $\varepsilon$ such that for any solution $\Phi_{t}\big(\varphi(0)\big)$ of system (2.1) the initial condition $\varphi(0)\in X_{0}$, one has

where $D$ is a distance function in $X_{0}$. Again, by proof by contradiction, suppose not, and then $D\big(\Phi_{t}(\varphi(0)), E_{0}\big)^{\infty} < \varepsilon_{1}$, $\forall\varepsilon_{1} > 0$. Hence, there must exist $T > 0$ such that

for $t > T$. For $t > T$, we have

where

We consider an auxiliary system as follows:

where $u = (u_{1}, u_{2}, u_{3})^{T}$ and

and we have $S\big(A(\varepsilon_{1})\big) = \max\big\{\rm{Re} \lambda: \lambda \big.$ is an eigenvalue of $\big.A(\varepsilon_{1})\big\}$, where $S\big(A(\varepsilon_{1})\big)$ represents the stability modulus of matrix $A(\varepsilon_{1})$. Note that if $A(\varepsilon_{1})$ is irreducible and has non-negative off-diagonal elements, then $S\big(A(\varepsilon_{1})\big)$ is a simple eigenvalue of $A(\varepsilon_{1})$ with a positive eigenvector (see Theorem A. 5 in ^[30]). According to Lemma 2.1 in ^[31] and the proof of Theorem 2 in ^[27], the following equivalent inequalities can be obtained:

Since $S\big(A(\varepsilon_{1})\big)$ is continuous in $\varepsilon_{1}$, we choose $\varepsilon_{1} > 0$ small enough such that $S\big(A(\varepsilon_{1})\big) > 0$ as $R_{0} > 1$. Hence, let $u(t) = \big(u_{1}(t), u_{2}(t), u_{3}(t)\big)^{T}$ be a positive solution of the auxiliary system (3.12), which is strictly increasing with $u_{i}(t)\rightarrow +\infty$ as $t\rightarrow +\infty, i = 1, 2, 3$. According to the comparison principle, we have

From system (2.1), since

we have

This contradicts our assumption. Hence, $E_{0}$ is an isolated invariant set in $X$ and $W^{s}(E_{0})\cap X_{0} = \phi$. Therefore, the solution of system (2.1) is uniformly persistent if $R_{0} > 1$.

□

4.   The optimal control problems

In the previous section, we discuss that the basic reproduction number $R_{0}$ is an important threshold to measure whether tuberculosis is dying out, and analyze the dynamic behavior of system (2.1). Because of current development of tuberculosis, China cannot achieve the goal of eliminating tuberculosis by 2035 set by WHO ^[24]. According to the existing literature, the prevention and control of tuberculosis in China still faces the problems of high transmission rate among susceptible and infected people, insufficient vaccination coverage, large number of latent individuals ^[24], and low clearance rate of contaminated environments ^[26]. Therefore, in the following study, referring to the research of Guo and Li ^[32] on optimal control of the new online game addiction model, we will introduce four control variables $u_{1}(t)$, $u_{2}(t)$, $u_{3}(t)$, and $u_{4}(t)$ to reduce the number of tuberculosis patients in China.

Among them, the control variable $u_{1}(t)$ represents raising awareness of TB prevention in uninfected populations, such as avoiding contact with active TB patients by susceptible or vaccinated persons. The control variable $u_{2}(t)$ represents increasing the vaccination rate of susceptible populations. The control variable $u_{3}(t)$ represents screening exposed individuals to reduce the proportion of individuals who develop active TB. The control variable $u_{4}(t)$ represents increasing the clearance rate of pathogens in contaminated environments, such as timely sterilization of the environment (doorknobs, towels, handkerchief, toys, utensils, beds and toilet seats, bathroom sink faucets, etc.) that TB patients have been exposed to. Based on the above assumptions, we establish the following optimal control model:

Next, we consider the control variable set $u(t) = \big(u_{1}(t), u_{2}(t), u_{3}(t), u_{4}(t)\big)\in U$ depending on the state variables $S, V, E, I, T, R,$ and $W$. Control variables are bounded and Lebesgue measurable on $[0, 1]$ with $U = \big\{(u_{1}, u_{2}, u_{3}, u_{4}) |0\le u_{i}\le 1, t\in [0, T_{0}], i = 1, 2, 3, 4\big\}$. The objective functional to be minimized is represented by

where $Z_{1}, Z_{2}$, and $Z_{3}$ respectively represent the weight constants of latent individuals, infected individuals, and treated individuals. $\xi_{1}, \xi_{2}$, $\xi_{3}$, and $\xi_{4}$ respectively represent the weight constants of transmission rate, vaccination rate, detection rate, and clearance rate. Our goal is to find an optimal set of control measures that minimizes the cost of control and the number of people treated for latent, active infections. The control variables can be represented by quadratic linear combinations, and the optimal control set $\big(u_{1}^{*}(t), u_{2}^{*}(t), u_{3}^{*}(t), u_{4}^{*}(t)\big)$ satisfies $J\big(u_{1}^{*}(t), u_{2}^{*}(t), u_{3}^{*}(t), u_{4}^{*}(t)\big) = \mathop{\min}\limits_{U}J\big(u_{1}(t), u_{2}(t), u_{3}(t), u_{4}(t)\big)$.

In order to solve the above problems, we derive the Hamiltonian function corresponding to the necessary conditions of optimal control according to Pontryagin's maximum principle ^[33]

where $\lambda_{i}(t), i = 1, 2, ..., 7$ are adjoint functions.

The Pontryagin maximum principle ^[33] shows that the optimal control $u(t)$ and the state variable $x(t)$ satisfy

and

That is to say, $(4.4)$ correspond to this paper, and we have

If we take the partial derivative of $\mathcal{H}$ and substitute it into the above equations, we have the following system:

Since the state variables in the controlled system have no terminal value, the transversality conditions satisfy $\lambda_{i}(T)$ $(i = 1, 2, ..., 7)$ at the final time $T_{0}$. From Pontryagin's maximum principle, we can see that the optimal control variables $u_{1}^{*}(t), u_{2}^{*}(t), u_{3}^{*}(t), u_{4}^{*}(t)$ satisfy $\frac{\partial\mathcal{H}(t, x, u, \lambda)}{\partial u} = 0$, and combined with $0\le u_{i}\le 1, i = 1, 2, 3, 4$, we can get

where

5.   A case study

In this section, we will take tuberculosis in China as an example and use the MCMC method to estimate parameters. Then, we will calculate the basic reproduction number based on known data and estimate parameters, and use LHS and PRCC to analyze the sensitivity of each parameter of the basic reproduction number $R_{0}$.

5.1.   Sources of data and parameter estimation

In this section, we first introduce some parameters of model (2.1) based on annual tuberculosis data and death data due to illness from 2005 to 2021 obtained from the China Center for Disease Control and Prevention ^[34]. Then, we use the MCMC method to estimate the other parameters of model (2.1). The number of tuberculosis cases and deaths from 2005 to 2021 are shown in Table 2.

Based on the above data, current research results, and tuberculosis epidemiology, we have the following settings:

1) In 2005, the total population of China was $1,307.56$ million ^[35], so we assume that the initial total population $N(0) = 1,307,560,000$.

2) According to the report in China Demographic Yearbook ^[35], the average life expectancy in China was $75.51$ years in 2021. Since the natural mortality rate can be calculated as the inverse of life expectancy, we take $d = \frac{1}{75.51}$.

3) With an average of $2,943.8$ TB deaths between 2005 and 2021, we estimate the average TB death rate to be $\mu = \frac{2,943.8}{1,307,560,000}\approx 2.251\times10^{-6}$.

4) We assume that the total population will not change significantly from 2005 to 2021. Therefore, in the numerical simulation, we assume that the population $N(t)$ is a constant, that is, $N(t) = N(0)$, thus obtaining the susceptible recruitment rate $\Lambda = dN(0)\approx17,316,382$.

5) According to the reports by Liu et al. ^[36] and Jiang et al. ^[24], the coverage rate of BCG vaccines in China can reach approximately $70\%$, so we assume that the vaccination rate is $k = 0.7$.

6) According to the data from the China Center for Disease Control and Prevention ^[34], the number of new tuberculosis cases in China in 2005 was $1,259,308$, so we set the initial value of active TB patients $I(0) = 1,259,308$.

7) In addition, we used the MCMC algorithm to estimate the transmission rates $\beta_{1}$, $\beta_{2}$, $\sigma_{1}$, and $\sigma_{2}$, the recovery rate of latent individuals $\gamma_{1}$, the recovery rate of infected individuals $\gamma_{2}$, the treatment recovery rate $\gamma_{3}$, the progression rate from latent individuals to infected individuals $v$, and the initial values $V(0)$, $E(0)$, $T(0)$, and $R(0)$. Then, $S(0)$ can be calculated from $S(0) = N(0)-\big(V(0)+E(0)+I(0)+T(0)+R(0)\big)$.

8) We assume that an infected person sheds one unit of the virus, which means we assume $\omega = 1$. Then, we simulate the value of $W(0)$ using the MCMC method.

The values, units, standard deviations, errors, and data sources of each parameter are shown in Table 3.

The fitting results of model (2.1) were compared with the case data of China from 2005 to 2021, as shown in Figure 2. Among them, the red dots represent the actual reported cases of infection, the black curve is the patients fitted through MCMC, and the gray part is 50, 90, 95, and $99\%$ confidence intervals from deep to shallow. We can see that all points are within a $95\%$ confidence interval, indicating good fitting performance.

The basic reproduction number $R_{0}$ is a very important index to measure disease outbreak or extinction. In this paper, we used $15,000$ Markov chain samples to simulate the value of $R_{0}$, as shown in Figure 3(a). According to the parameters in Table 3 and the calculation method of the basic reproduction number, we determine that the average value of $R_{1}$ is $0.7046$, the average value of $R_{2}$ is $1.4056$, and the average value of $R_{0}$ is $2.1102$. In addition, $R_{0}$ satisfies a normal distribution, as shown in Figure 3(b), and all $R_{0}$ values are greater than 1, that is to say, tuberculosis will not disappear over time in China.

5.2.   Sensitivity analysis

In this part, we used the LHS method and PRCC ^[38] to analyze the uncertainty and sensitivity of $R_{0}$ relative to each parameter, as shown in Figure 4. As can be seen from the figure, the parameters that are positively correlated with the basic reproduction number $R_{0}$ are $\beta_{1}$, $\beta_{2}$, $\sigma_{1}$, $\sigma_{2}$, and $v$, respectively, and the parameters that are negatively correlated with the recovery rate are $\gamma_{1}$, $\gamma_{2}$, and $\gamma_{3}$. We believe that an absolute value of PRCC greater than $0.4$ indicates a high correlation between the input parameter and the output variable, a value between $0.2$ and $0.4$ indicates a moderate correlation, and a value between 0 and $0.2$ indicates a weak correlation. So, we can see that $\beta_{2}$, $\gamma_{1}$, $\gamma_{2}$, $\sigma_{2}$, and $v$ are highly correlated with $R_{0}$.

5.3.   Numerical simulation of optimal control

In this section, we first introduce the cost of the control strategy proposed in the fourth part. The research of Xue et al. ^[37] shows that the average cost of improving the awareness of uninfected people on tuberculosis prevention is 2 ＄/person/month, that is, the average cost of $u_{1}$ is 24 ＄/person/year. We assume that the average cost of clearing pathogens of contaminated environments is 60% of the average cost of $u_{1}$, which means that the average cost of $u_{4}$ is 14.4 ＄/person/year. According to the survey, the cost of receiving the BCG vaccine in China is approximately 200–400 ¥, with a median of 300. The protection period of the BCG vaccine is 10–20 years ^[3], and we assume that the protection period is 15 years. Therefore, the average cost of improving vaccine coverage is approximately 2.82 ＄/person/year. The routine examination items of tuberculosis are lung CT, sputum smear, etc. Pan et al. ^[39] show that the cost for these routine examinations is 420 to 2000 ¥. Due to the large number of tuberculosis latent individuals in China, considering the actual financial problems, we take the lowest cost here, that is, 420 ¥/person/year$\approx$59 ＄/person/year.

Based on the above settings, we will use the fourth-order Runge Kutta method to numerically simulate the system (2.1) and the accompanying system (4.7). We choose the weight constant values of the objective function as $Z_{1} = 10$, $Z_{2} = 10$, $Z_{3} = 10$, $\xi_{1} = 24$, $\xi_{2} = 2.82$, $\xi_{3} = 59$, and $\xi_{4} = 14.4$. To show more clearly the impact of different control strategies on TB in China, here we compare different control strategies, namely the single control strategy (implementing only one control), dual control strategy (implementing two control policies at the same time), triple control strategy (implementing three controls at the same time), and quadruple control strategy (four types of controls implemented simultaneously).

As shown in Figure 5, we compared the changing trends of infected individuals when they were not under control and when they were under four single control strategies. The dark blue solid line indicates that it is without control. The solid blue line represents $u_{1} = 0.5$. The red dashed line indicates $u_{2} = 0.5$. The solid black line represents $u_{3} = 0.5$. The solid green line represents $u_{4} = 0.5$. According to our prediction, there will still be 303,100 cases of tuberculosis in China by 2035. Although it has decreased by 65% compared to 2015, it is far from meeting the WHO target of a 90% decrease compared to the 2015 baseline. From the figure, we can see that the four control strategies can effectively reduce the number of tuberculosis cases. However, in the long run, increasing the vaccination rate ($u_{2}$) and clearing pathogens ($u_{4}$) from contaminated environments are more effective. Unfortunately, although the four single control strategies have effectively reduced the number of tuberculosis cases, they have not achieved the goals set by WHO.

As shown in Figure 6, we compared the changing trends of infected individuals when they were not under control and when they were under six dual control strategies. The dark blue solid line indicates uncontrolled. The blue dashed line indicates $u_{1} = u_{2} = 0.5$. The red dashed line indicates $u_{1} = u_{3} = 0.5$. The solid red line indicates $u_{1} = u_{4} = 0.5$. The solid green line represents the implementation strategy $u_{2} = u_{3} = 0.5$. The purple dashed line represents the implementation strategy $u_{2} = u_{4} = 0.5$. The solid black line represents the implementation strategy $u_{3} = u_{4} = 0.5$. We can see that the six dual control strategies can effectively reduce the number of tuberculosis cases, and the best effect is to implement strategies $u_{2}$ and $u_{4}$ at the same time. Compared to Figure 5, we can see that, overall, the effect of the dual control strategy is better than that of the single control strategy. Unfortunately, all dual control strategies still fail to achieve the goals proposed by the WHO.

Figure 7 is a comparison between uncontrolled and four triple control strategies. The dark blue solid line represents uncontrolled. The black dashed line represents $u_{1} = u_{2} = u_{3} = 0.5$. The red dashed line represents $u_{1} = u_{3} = u_{4} = 0.5$. The green dashed line represents $u_{1} = u_{2} = u_{4} = 0.5$. The solid black line represents $u_{2} = u_{3} = u_{4} = 0.5$. We can see that the changing trends of the four triple control strategies are similar. And, if $u_{1}$, $u_{2}$, $u_{4}$ or $u_{2}$, $u_{3}$, $u_{4}$ are controlled simultaneously, the goals proposed by the WHO will be achieved.

Figure 8 is a comparison diagram of the quadruple control strategy. The dark blue solid line indicates that it is not under control. The black dashed line represents $u_{1} = u_{2} = u_{3} = u_{4} = 0.5$. The green dashed line indicates the implementation of optimal control. We can see that when implementing the quadruple control strategy or optimal control, the number of infected individuals initially rapidly decreases, but gradually decreases after 2024 until it approaches zero. And, both the quadruple control strategy and optimal control can achieve the goals proposed by the WHO. Finally, compared with Figures 5–8, we suggest that the control strategy be implemented as soon as possible, and the four control strategies should be implemented at the same time.

5.4.   Cost-effectiveness analysis

In this section, we conduct a cost-benefit analysis of the control strategy proposed in the previous section. The analysis is mainly carried out from three aspects: total cost, unit cost (the ratio of total cost to total number of infections avoided), and incremental cost-benefit ratio (ICER).

The ICER can be expressed as

where the total costs are calculated by

Comparing Figures 5–8, it can be seen that although the single control strategy and the dual control strategy have reduced the number of tuberculosis cases, they have not reached the goals proposed by the WHO. Therefore, in this section we only conduct a cost-benefit analysis on the triple control strategy and the quadruple control strategy. We have developed control strategies in ascending order based on the number of infections avoided, as shown in Table 4. The total number of people who avoid infection refers to the difference between the number of uncontrolled infected individuals and the number of controlled infected individuals.

As shown in Table 4, we will record different strategy combinations as strategies A, B, C, D, E, and F. We can see that strategies C, D, E, and F can all achieve the goals proposed by the WHO, and strategy F has the highest total number of infection averted. From the total cost analysis and unit cost analysis, strategy D has the lowest total cost and unit cost, and can achieve the goal. The total cost of strategies A, B, C, and E exceeds 100 billion US dollars. According to the research of Li et al. ^[40], this is far more than China's actual investment in tuberculosis.

As for ICER, we can see from Table 4 that strategy E has the highest ICER value, indicating that, compared to strategy D, strategy E has higher costs and fewer outputs. Therefore, we delete strategy E. We then recalculate the ICER value of strategy F to $3.12\times10^{4}$ (deleting strategy E will only affect strategy F, therefore only recalculate strategy F). Due to the fact that ICER(F)$> \ $ICER(D), we deleted strategy F. Using a similar method, we delete strategies A, B, and C. Finally, we conclude that, starting from 2020, strategy D is the most cost-effective control method.

6.   Conclusions and discussion

TB is an ancient, chronic infectious disease that has affected humans for thousands of years, spread primarily through respiratory droplets. In 2014, the WHO proposed the end TB strategy, which aims to end the global TB epidemic by 2035. However, most of the existing mathematical models consider a single factor. Therefore, in order to eliminate TB as soon as possible, we set up a new $SVEITRW$ model with seven compartments based on TB data in China from 2005 to 2021, and analyze the combined effects of vaccination, treatment effect, and contaminated environments on TB. It is proved that when $R_{0} < 1$, the disease-free equilibrium point of the disease is globally stable, and when $R_{0} > 1$, system (2.1) is uniformly persistent. Then, the MCMC method is used to fit the model parameters, and the average value of the basic reproduction number of tuberculosis in China from 2005 to 2021 is calculated to be 2.1102. This means that TB will not disappear in China if we do not control it. In addition, we also discuss the distribution of $R_{0}$. As shown in Figure 3(b), $R_{0}$ is normally distributed, and the value of $R_{0}$ is almost greater than 1, that is, tuberculosis will spread in China. We perform sensitivity analysis of the parameters using PRCC, and find that $\beta_{1}$, $\beta_{2}$, $\sigma_{1}$, $\sigma_{2}$, and $v$ are positively correlated with $R_{0}$, and the recovery rate $\gamma_{1}$ and $\gamma_{2}$ are negatively correlated with $R_{0}$. Finally, we conduct an optimal control and cost-benefit analysis. Based on practical factors, we consider the improvement of the tuberculosis prevention awareness of uninfected people, the vaccination rate, the screening of potential individuals to reduce the proportion of individuals who develop active tuberculosis, and the clearance rate of pathogens in the contaminated environment as four control strategies. The analysis shows that implementing four control strategies simultaneously has the best effect, and implementing $u_{2}, u_{3}, u_{4}$ simultaneously is the most cost-effective. Therefore, it is recommended to focus on improving vaccine coverage, screening potential individuals, and improving the clearance rate of pathogens in contaminated environments, and to implement control strategies as soon as possible.

TB is a disease with extremely complex transmission patterns, and our model takes into account both vaccination, contaminated environments, treatment effectiveness, and fast and slow progression compared to existing models. However, due to the complexity of its transmission patterns, it is difficult to fully consider all aspects of its epidemiology. For example, our data is based on years and ignores the periodicity of TB. In fact, according to the study of Xue et al. ^[37], the incidence of tuberculosis has a strong periodicity, which is generally high in spring and low in winter. Therefore, based on this study, we can adopt a model with a periodic propagation rate, consider the epidemiological characteristics of tuberculosis such as seasonality, drug resistance, exogenous reinfection and endogenous reactivation, further analyze the epidemic trend of tuberculosis in China, formulate corresponding control strategies, and consider the cost-effectiveness of strategies. These are issues left for future work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the NNSF of China (12361101).

Conflict of interests

The authors declare that there is no conflict of interest in this work.