Global stability and optimal control for a tuberculosis model with vaccination and treatment

We formulate a mathematical model to explore the impact of 
vaccination and treatment on the transmission dynamics of 
tuberculosis (TB). We develop a technique to prove that the basic 
reproduction number is the threshold of global stability of the 
disease-free and endemic equilibria. We then incorporate a control 
term and evaluate the cost of control strategies, and then perform 
an optimal control analysis by Pontryagin's maximum principle. Our 
numerical simulations suggest that the maximum vaccination strategy 
should be enforced regardless of its efficacy.


1.
Introduction. Tuberculosis (TB) is a widespread and deadly infectious disease mainly caused by Mycobacterium tuberculosis (Mtb). TB is an airborne disease transmitted through the coughing, sneezing, and spitting of its infectious individuals. One third of the world's current population has been infected with Mtb, and new infections occur at a rate of one per second. Fortunately, most of the cases will not develop into the full-blown disease, only about one in ten of these latent infections will eventually progress into an active disease. It was estimated that, in 2013, there are 9.0 million new TB cases and 1.5 million TB deaths (in which 0.4 million are HIV-positive patients) [20].
There have been several studies about modeling TB transmission with different prevention and control strategies, such as vaccination, treatment, quarantine and isolation [5,6,7,8]. But none of these models has taken into consideration of a combination of control methods. Based on the incomplete treatment TB model [21], we formulate here a model with vaccination and treatment. Evidence has shown that the protective efficacy of Bacillus Calmette-Grin (BCG) in the prevention of some serious TB types (e.g. meningitis) in children is about 75% [2]. However, BCG's efficacy on preventing pulmonary TB in adults is highly variable [20]. In our model, we assume that BCG is only given to newborns, and we use the model to explore the effectiveness. Our optimal control analysis provides a balanced strategy which coordinates individual's willingness on reducing the infected cases and minimizing the cost of control strategies.
This paper consists of four sections. In section 2, the model about TB transmission with vaccination and treatment is first presented, then equilibria, basic reproduction number, and global stability proof are given. In section 3, we extend the model with specific control strategies, use the Pontryagin's maximum principle in the optimal control analysis, and provide the simulation results. We conclude our findings in section 4.

2.
A TB model with vaccination and treatment.
2.1. Model formulation. The total population is partitioned into five compartments: susceptible (S); latent (L); infectious (I); treatment (T ), and vaccination (V ), in which individuals had been vaccinated when they were newborns and still possess partially immunity against TB. The vaccinated individuals may get infected through contact with infectious individuals and then become latent for TB, meanwhile, vaccinated individuals may lose the immunity and become susceptible for TB. In this model, we also consider that some individuals in the infectious compartment I may reenter the latent compartment directly due to self-cure.
The flowchart among compartments is schematically depicted in Figure 1. Figure 1. The transmission diagram for model (1).
Here, Λ is the recruitment rate, p (0 ≤ p ≤ 1) is the fraction of the newborns vaccinated, µ is the natural death rate coefficient, α 1 and α 2 are the disease-induced death rate coefficients for individuals in compartments I and T respectively. β is the transmission coefficient, and then, the rate at which the susceptible individuals are infected is βIS; the rate at which the vaccinated individuals are infected is σβIV . Here 0 ≤ σ ≤ 1, σ = 0 indicates that the vaccination protection rate is 100%, σ = 1 indicates that the vaccination protection efficacy is 0, that is, 1 − σ represents the reduction in infection risk due to vaccination efficacy. c (c ≥ 0) is the immunity waning rate coefficient in the vaccinated host. ε is the rate coefficient at which an individual leaves the latent compartment to become infectious. γ is the rate at which an infective individual leaves the infectious compartment I. q (0 ≤ q < 1) is the fraction of infectious individuals reentering the latent compartment through self-cure, while 1 − q is the fraction of infectious individuals transferred to the treatment compartment. δ is the rate at which a treated individual who leaves the compartment T . k (0 ≤ k ≤ 1) is the fraction of the success treatment, that is, the treated individuals will enter latent compartment L after treated; correspondingly, 1 − k is the fraction of the drug-resistant individuals among the treated individuals, that is, the treated individuals will enter infectious compartment I after treated, which reflects the treatment failure.
Denote Figure 1 leads to the following model of ordinary differential equations: (1) From system (1), it follows that (S + L It implies that the region Ω = {(S, L, I, T, V ) ∈ R 5 + : S + L + I + T + V ≤ Λ/µ} is a positively invariant set for the system (1). The dynamics of system (1) on the set Ω is taken into account in the following part.

2.3.
Global stability of equilibria. With respect to the global stability of equilibria of system (1), we have the following results.
Proof. Firstly, we prove the global asymptotically stability of the disease-free equi- Then Since S 0 and V 0 satisfy the following equations we can rewrite the first and fifth equation of system (1) as Define the Lyapunov function Then the derivative of W 1 with respect to t along solution of system (1) is given by Substituting the expressions of S 0 and V 0 into function F(S,V) yields is the singleton {E 0 }. Therefore, it follows from the LaSalle's invariable principle [12] that E 0 is globally asymptotically stable on the set Ω if R 0 ≤ 1. Secondly, we prove the global asymptotically stability of the endemic equilibrium E * (S * , L * , I * , T * , V * ).
Define function where then the derivative of H 2 with respect to t along solution of system (17) is given by Direct calculations show that and Furthermore, to prove that function G(x, y, z, u, v) is negative definite with respect to x = 1, y = 1, z = 1, u = 1 and v = 1, we use equations (4), and consider the following three cases: If βI * < cV * S * , function G(x, y, z, u, v) will be expressed as If βI * = cV * S * , function G(x, y, z, u, v) will be expressed as If βI * > cV * S * , function G(x, y, z, u, v) will be expressed as Since the arithmetical mean is greater than, or equal to the geometrical mean, we have H 2 | (17) ≤ 0 for x, y, z, u, v > 0 and H 2 | (17) = 0 if and only if x = 1, v = 1 and y = z = u.
Theorem 2.2 indicates that the basic reproduction number, R 0 , plays a crucial role in determining dynamics of the disease and could be used for evaluating control strategies. Note that, for system (1), k reflects the fraction of successful treatment. Direct calculations show that ∂R 0 /∂k| (1) < 0. Therefore, increasing the fraction of successful treatment coefficient helps to reduce TB infection. It seems that the increase of k mainly depends on the decrease of drug-resistance. However, much higher levels of resistance and poor treatment outcomes are of major concern in some parts of the world, even for multidrug-resistant tuberculosis (MDR-TB) patients, the major treatment gaps remain and funding is insufficient [20].
Direct calculations show that ∂R 0 /∂p| (1) < 0. Therefore, increasing the value of p helps to reduce the number of infections, that is, increasing the vaccination fraction of infant susceptible individuals against infection of TB is helpful in controlling the spread of TB. But efficacy of BCG vaccine, even for all kinds of vaccines which exist against TB, is still a question. Successful TB vaccine research and development will require investments well into the next decade [9].

3.
Optimal control of extended model.

3.
1. An extended TB model. Despite substantial growth in funding for TB since 2002, an annual gap of around US$ 2 billion still needs to be filled [20]. In this section, an optimal control model for the transmission dynamics of TB is formulated by extending the autonomous system (1) to include control functions, and our goal here is to study the optimal control strategies to curtail the epidemic and keep cost minimum. The non-autonomous system of differential equations is as follows: In system (27), the control function u 1 (t) is the fraction of the newborns who enter the vaccination compartment. The control function u 2 (t) is the fraction of the registered treated individuals who enter the latent compartment L. Correspondingly, 1 − u 2 (t) denotes the fraction of the registered treated individuals who enter the infectious compartment I. The control functions u 1 (t), u 2 (t) are bounded (a 1 ≤ u 1 (t) ≤ b 1 , a 2 ≤ u 1 (t) ≤ b 2 ) and Lebesgue integrable.
Here, the objective functional is defined as dt, (28) which is subjected to the state system (27). The goal is to minimize the infected populations (including latent, infectious, being treated populations), and the costs of implementing these two control strategies. The cost results from a variety of sources, assuming that the costs of vaccination and treatment are nonlinear, and then, quadratic form is taken here [11,17]. The coefficients, B 1 , B 2 , B 3 and C 1 , C 2 represent, respectively, the corresponding weight constants, and these weights are balancing cost factors due to size and importance of the parts of the objective functional.
An optimal control function (u * 1 (t), u * 2 (t)), such that reaches its minimum, where and a i , b i (i = 1, 2) are fixed positive constants.

3.2.
Analysis of optimal controls. The necessary conditions is that an optimal solution must satisfy the Pontryagin's maximum principle [14]. This principle converts (27), (28) and (29) into a problem of minimizing pointwise a Hamiltonian H, with respect to (u 1 (t), u 2 (t)): where l i (i = 1, . . . , 5) are the adjoint variables. By the Pontryagin's maximum principle and the existence result for the optimal control pairs from [10], we obtain the following theorem.
Theorem 3.1. There exists an optimal control (u * 1 (t), u * 2 (t)) and corresponding solution, S, L, I, T and V such that minimizes J(u 1 (t), u 2 (t)) over Ω. Furthermore, there exists adjoint functions, l 1 (t), . . . , l 5 (t), such that With transversality conditions the following characterization holds Proof. The existence of an optimal control can be obtained due to the convexity of integrand of J(u 1 (t), u 2 (t)) with respect to (u 1 (t), u 2 (t)) [10], a priori boundedness of the state solutions, and the Lipschitz property of the state system with respect to the state variables. The following is obtained by the Pontryagin's maximum principle [14]: and evaluated at the optimal control and corresponding states, which results in the stated adjoint system (32) with transversality (33). Consider the optimality conditions, ∂H/∂u 1 = 0, ∂H/∂u 2 = 0 and work out (u 1 (t), u 2 (t)), subject to the constraints, the characterizations (34) can be derived.
4. Numerical simulations. Here, we use an iterative method for resolving the optimal system by a fourth-order Runge-Kutta scheme. The optimal system, which is consisted of 10 ordinary differential equations from the state equations (27) and adjoint equations (32), and coupled with the control characterization (34). With initial conditions for the states, given an initial assumption for the controls, the state differential equations (27) are solved forward in time by a fourth-order Runge-Kutta scheme. With the application of the current iteration solution of the state equations (27) and the transversality conditions (33), the adjoint system (32) is solved backward in time, again, by a fourth-order Runge-Kutta scheme. Both state and adjoint values are used to update the controls with the characterizations given by (34), and then the iterative process is repeated. This iterative process will terminate if current state, adjoint, and control values converge sufficiently. The bounds on the controls come from the estimation on the effectiveness of the intervention strategy. For the control u 1 (t), it is reasonable to assume that its upper bound is 1.0 and lower bound is 0; for the control u 2 (t), it is reasonable to assume that upper bound is 0.75 and lower bound is 0. The initial values for the states are given by S(0) = 5000, L(0) = 2000, I(0) = 2000, T (0) = 1000, V (0) = 95000, respectively. The weights in the objective function are B 1 = 1, B 2 = 1000 and B 3 = 2000, which means that the minimization of the number of latent, infectious, and treated has different importance. Now, for the optimal system, an optimal control strategy in the case of C 1 = 10 3 , C 2 = 10 6 is studied. In Figure 2, the population with optimal control actions are shown in dashed lines, compared with the populations without control actions which are shown in solid lines. Its parameter values are presented in the Table 1. Lack of data, some parameter values are assumed within realistic ranges, and its time unit here is year. As depicted in Figure 2, the number of infected individuals decrease sharply when the optimal control strategies are used, compared to those of infected individuals without control actions.
In Figure 3, we consider different self-cure fraction q. As q increases, the stress for the treatment strategy will drop. Since q depends on the immunity of TB infected individual, we should consider other co-infective diseases which may decrease the   Table 1.
immunity of TB infected individual if we want to control TB disease, such as HIV.
And it is well known that TB/HIV is one of the most deadly disease in the world. Direct calculation shows that ∂R 0 /∂σ| (1) > 0, then increasing the protective efficacy of BCG for preventing serious forms of TB, that is, decreasing the value of σ, is an effective strategy for controlling TB. Furthermore, from Figure 4, for different σ, despite of the level of its efficacy, the vaccination strategies will be  fully adopted almost during all the procedure while the treatment strategies will be effected by the different σ values. Meantime, the most effective treatment strategy should be taken into account at the onset of the control. As σ becomes smaller, that is, the vaccination protection is more effective, the treatment will relieve the stress. Only if the vaccination does not have any protection, we should not adopt it. It is in coincidence with the fact that it always is suggested to adopt the vaccination strategy in many TB high prevalence regions or groups.  Furthermore, we consider different costs for the vaccination and successful treatment. C 1 reflects the relative costs of the vaccination. Figure 5 shows how the control strategies depend on the various weights of C 1 : as C 1 increases, the amount of u 1 (t) decreases. C 2 reflects the relative costs of the successful treatment. Figure  6 illustrates how the control strategies depend on the various weights of C 2 : as C 2 increases, the amount of u 2 (t) decreases.

Conclusion.
A TB model has been formulated and its global dynamical behaviors have been studied. Since we know only the existence of endemic equilibrium if R 0 > 1, it is very tricky to determine the negative definiteness or semi-definiteness of the derivative of the Lyapunov function. So, we have regrouped the derivative into different forms according to various parameter combinations. Then we have obtained the conclusion that the disease-free equilibrium is globally asymptotically stable if R 0 ≤ 1; the endemic equilibrium is globally asymptotically stable if R 0 > 1. It indicates that the basic reproduction number, R 0 , plays a crucial role in determining dynamics of the disease and could be used for controlling the spread of TB. R 0 is induced by vaccination and treatment is also discussed. To control the spread of TB, the optimal regulations of an extended model with vaccination and treatment has been proposed, and optimal control strategies for several scenarios have been identified. Numerical simulations of the optimal system of extended TB model (27) indicate that the vaccination and treatment strategies are effective in reducing TB disease transmission. Especially, if the BCG vaccine works, we should adopt maximum strategy during the procedure. The simulation results of optimal control also show that a cost-effective control strategy depends on the cost of the implementation of the control.