key: cord-268959-wh28s0ws
authors: Gao, Da-peng; Huang, Nan-jing
title: Optimal control analysis of a tuberculosis model()()
date: 2017-12-29
journal: Appl Math Model
DOI: 10.1016/j.apm.2017.12.027
sha: 
doc_id: 268959
cord_uid: wh28s0ws

In this paper, we extend the model of Liu and Zhang (Math Comput Model 54:836-845, 2011) by incorporating three control terms and apply optimal control theory to the resulting model. Optimal control strategies are proposed to minimize both the disease burden and the intervention cost. We prove the existence and uniqueness of optimal control paths and obtain these optimal paths analytically using Pontryagin’s Maximum Principle. We analyse our results numerically to compare various strategies of proposed controls. It is observed that implementation of three controls is most effective and less expensive among all the strategies. Thus, we conclude that in order to reduce tuberculosis threat all the three controls must be taken into consideration concurrently.

Tuberculosis (TB) is a major global health threat caused by Mycobacterium tuberculosis (Mtb) that most often affects the lungs (pulmonary TB) but can affect other sites as well (extrapulmonary TB). Fortunately, TB is a preventable and curable disease. Following the World Health Organization (WHO), more than 51 million patients have been successfully treated in countries that adopted the WHO strategy since 1995 [1] . However, the global burden of TB remains enormous. In 2011, there were an estimated 8.7 million new cases of TB (13% co-infected with HIV) and 1.4 million people died from TB [1] . In 2013, there were 9 million new TB cases and 1.5 million TB deaths (in which 0.4 million are HIV-positive patients) [2] . The increase of new cases lies in multiple factors such as the spread of HIV, the collapse of public health programs, the emergence of drug-resistant strains of Mtb [3] [4] [5] and exogenous re-infections, where a latently-infected individual acquires a new infection from another infectious individual (see, for example, [6, 7] and references therein). The TB control approach is the chemoprophylaxis treatment in the absence of an effective vaccine. Neglecting the compliance with drug treatments might result in a relapse and antibiotic resistant TB, i.e. multidrug-resistant TB (MDR-TB) [8] . MDR-TB is one of the most serious health problems and the progress in curing that is slow. Critical funding gaps for TB care and control make it difficult to sustain recent gains, promote further research and develop new drugs and vaccines [1] .

It is well known that the mathematical models play vital roles in the dynamics and control of many epidemics including malaria, severe acute respiratory syndromes (SARS) and TB (see, for example, [9] [10] [11] [12] [13] [14] [15] [16] [17] and references therein). Many mathematical dynamic models for TB have been studied extensively in the literature; for instance we refer the reader to [18] [19] [20] [21] [22] [23] . "Fast progressors" and "slow progressors" are two different ways considered by most models to progress to active disease after infection. The probability of latent and treated individuals who might be reinfected [24] is also taken into account by some latest models as the infection and/or disease do not confer full protection which is already well recognized [25] . Some previous models of TB, especially the predictive models endeavor to calculate a threshold called basic reproductive number. The dynamics of transmission of the disease has been analyzed in terms of the reduction of the basic reproductive number. Since the first paper published by Jung et al. [8] in 2002, the time dependent optimal control strategies have been employed in the study of dynamics of TB mathematical models by many authors (see, for example, [8, [26] [27] [28] [29] [30] [31] [32] [33] [34] ). Valuable theoretical results generated by both approaches of studying control strategies can be used to guide epidemic control programs. Various objective criteria may be adopted, according to a chosen goal (or goals).

In 2011, considering that the treated individuals are also infectious, Liu and Zhang [24] proposed a nonlinear TB model system as follows:

where a population of size N ( t ) is partitioned into five subclasses: susceptible ( S ( t )), vaccinated ( V ( t )), individuals infected with TB in latent stage ( L ( t )), individuals infected with TB in the active stage ( I ( t )) and treated individuals infected with TB ( T ( t )). The parameter represents the constant recruitment rate of susceptible individuals. Parameter p is the rate at which susceptible individuals are moved into the vaccination process. The natural death rate is μ. The disease-induced death rate coefficient for individuals in compartment I is α. Parameter ρ represents the rate that individuals successfully treated of TB return to the latent TB stage and parameter γ is the treatment rate in the infective class. It is assumed that the time before latent individuals become infectious obeys an exponential distribution, with mean waiting time 1 δ . Thus, δ is the rate at which an individual leaves the latent compartment to become infectious. Susceptible individuals acquire TB infection from individuals with active TB at rate βS(I + ρ 1 T ) , where β is the disease transmission coefficient, and the parameter ρ 1 < 1 accounts for the reduction in infectiousness among individuals with active TB who are treated (in comparison to those who are not treated). Fraction l of susceptible individuals who acquire TB infection is assumed to enter the latent TB class ( L ), at the rate βS(I + ρ 1 T ) , while the remainder moves to the active TB class ( I ). It is assumed that individuals in the latent class do not transmit infection. The vaccinated individuals are infected at rate ρ 2 βV (I + ρ 1 T ) , where ρ 2 < 1 represents the reduction in risk of infection owing to vaccination (for a vaccine that offers 100% protection, ρ 2 = 0 ; thus in reality 0 < ρ 2 < 1). For model (1.1) , Liu and Zhang [24] first obtained the basic reproduction number

.

After that they showed that if the basic reproduction number R 0 of the model system is below one then disease dies out otherwise it persists. This indicates that the basic reproduction number, R 0 , plays a crucial role in characterizing dynamics of the disease and could be used to suggest or design TB control strategies. However, they did not further investigate time dependent control strategies since their discussions had been concentrated on prevalence of TB at equilibria. Note that, ∂ p | (1 . 1) < 0 , also due to Liu and Zhang [24] , which implies that increasing the vaccination rate of susceptible individuals against infection of TB is conducive to disease control. In order to understand the dynamics under what circumstances TB can be controlled or curtailed, we implement the theory of optimal control. Three intervention strategies, called controls, are included in model system (1.1) . We write controls as functions of time and assign reasonable upper and lower bounds to them. First, a "case finding" (identification of latently infected individuals) control mechanism is incorporated in model (1.1) by replacing constant vaccination rate with u 1 ( t ). Second, a "case holding" (treatment of individuals with active TB) control mechanism is incorporated in model (1.1) by replacing constant successfully treatment rate with u 2 ( t ). Finally, another "case holding" is instituted in model (1.1) by replacing constant treatment rate with u 3 ( t ). It is required that u 1 ( t ), u 2 ( t ) and u 3 ( t ) are bounded Lebesgues integrable functions. We intend to determine optimal control strategies that minimize not only the infected individuals but also the cost of these three controls. As mentioned by Kumar and Srivastava [14] , the larger proportion of population is covered, the higher efforts and expenditure is required to provide vaccination. Hence, it is suitable to use the nonlinearity of order four to describe the high expensiveness concerned with the vaccination process. Following the idea of Kumar and Srivastava [14] , in this paper, the cost constructions that account for high nonlinear relationship (nonlinearity of order four i.e. u 4 1 (t) ) between cost and effort s during vaccination process has been employed in the application of optimal control to TB mathematical models.

Using the same parameters and class names as in model system (1.1) , the system of differential equations describing the controlled model is as follows:

(1. 2) In view of the fact that vaccination and treatment are easily available and implementable control strategies to stop TB, we obtain model system (1.2) by considering vaccination and treatment as control variables from model system (1.1) . The detailed analysis of optimal control problem is given in the following sections. The subsequent part of this paper is as follows: we present the formulation of the optimal control problems in Section 2 . In Section 3 , we investigate the existence of an optimal control function and derive an optimal system characterizing the optimal control. In Section 4 , we analyze the uniqueness of the optimal system, in which the state system coupled with co-states. It is worth pointing out that the proof of existence and uniqueness involving biquadratic form of the control variables has not been explored so far as we know. We perform numerical simulations to illustrate our theoretical results in Section 5 . The paper ends with a conclusion in Section 6 .

This section is devoted to the study of our model system (1.2) , when we administered vaccination and treatment policies over a fixed time window. For that we design an objective functional by keeping the biomedical goal in our mind. Our goal is to minimize the number of infected individuals (including latent, infectious, being treated individuals) with TB virus while at the same time keeping the cost of implementing these three control strategies very low. In mathematical perspective, for a fixed terminal time t f , the problem is to minimize the objective functional

Here, the total cost on a finite time horizon [0, t f ] consists of the cost induced by the disease itself and the cost induced by vaccination and treatment effort s. We split the disease cost into the cost induced by latent group,

proportional to the number of individuals infected with TB in latent stage, the cost induced by infected group,

proportional to the number of individuals infected with TB in the active stage, and the cost induced by treated group,

proportional to the number of treated individuals infected with TB. Modeling the control cost is a more delicate issue. The cost involved in vaccination process requires relatively higher efforts and expenditure because of covering larger proportion of population. For example, vaccination efforts target first quarter of population will be less than the next quarter and so on. Therefore, the cost for covering larger proportion of population involved in vaccination process will be not only increasing but also growth of cost will be faster and steeper. Hence, we consider the biquadratic form in the vaccination control to represent high expensiveness of vaccination process [16] . That is, the cost involved in vaccination process is taken as

It is also assumed that cost of treatment is nonlinear and takes a quadratic nature, which is found to be consistent with previous works in the literature (see, e.g. [8] ). In other words, the cost incurred in improving successfully treatment rate (the effort s such as regularity of drug intake) is taken as

Similarly, the cost incurred in providing treatment is taken as 6) are the weight constants of the infected TB individuals and control measures. They can be chosen to balance cost factors due to size and importance of the parts of the objective functional. We assume that there are practical limitations on the maximum rate at which individuals may be vaccinated or treated during a given period of time and the maximum successfully treated rate. In the rest of this paper, for simplicity of notation, we represent u 1 (t) = u 1 , u 2 (t) = u 2 and u 3 (t) = u 3 . We seek optimal controls u 1 , u 2 and u 3 in U ad such that

is the control set. In optimal control theory, there are several basic problems included such as to prove the existence of an optimal control, to characterize the optimal control, to prove the uniqueness of the control, to compute the optimal control numerically and to investigate how the optimal control depends on various parameters in the model.

In this section, we will study the sufficient condition for the existence of an optimal control of our system of ordinary differential equations (1.2) . We refer to the conditions in Theorem III 4.1 and its corresponding Corollary in [35] . After that, we characterize the optimal control functions by using Pontryagin's Maximum Principle [36, 37] , and then we derive necessary conditions for our control problem.

In this subsection, we will investigate the existence of an optimal control of our model (1.2) . The boundedness of solutions to system (1.2) for finite time interval is needed to establish the existence of an optimal control and the proof of the uniqueness of the optimality system. We begin by examining the priori boundedness of the state solutions of model (1.2) .

For model system (1.2) to be epidemiologically meaningful, it is important to prove that all its state variables are nonnegative for all time. In other words, solutions of model system (1.2) with nonnegative initial data remain nonnegative for all time t > 0. Suppose, for example, the variable S becomes zero for some time

To establish the upper bounds for the solutions, we consider an equation for the total population size N . The rate of change of the total population, obtained by adding all the equations in model (1.2) , is given by 

Proof. To prove this theorem, we follow the requirements from Theorem 4.1 and Corollary 4.1 in [35] and verify nontrivial requirements. Let r(t, x , u ) be the right-hand side of (1.2) . We need to show the following conditions are satisfied:

(1) r is of class C 1 and there exists a constant C such that

(2) The admissible set F of all solutions to system (1.2) with corresponding control in U ad is nonempty; In order to verify above conditions, we write

Then it is easy to see that r(t, x , u ) is of class C 1 and | r(t, 0 , 0) | = . Moreover, one has

Since S, V, L, I and T are bounded, there exists a constant C such that

This means that condition (1) holds. Thanks to condition (1), there exists a unique solution to system (1.2) for a constant control, which further implies that condition (2) holds.

In addition,

Condition (4) is obvious from the definition and the proof of this theorem can be completed by verifying condition (5) .

In order to show the convexity of the integrand in the objective functional g(t, x , u ) , we have to prove that

and so the proof is complete. 

Since there exists an optimal control for minimizing function (2.2) subject to system (1.2) , we drive the necessary conditions for this optimal control by using Pontryagin's Maximum Principle [36, 37] . We need to define Hamiltonian for deriving the necessary conditions as, Here, λ = (λ 1 , λ 2 , λ 3 , λ 4 , λ 5 ) ∈ R 5 is known as adjoint variable. The optimality system of equations is found by taking the appropriate partial derivatives of Hamiltonian (3.1) with respect to the associated state variables.

The following theorem is a consequence of the maximum principle.

Given an optimal control pair (u * 1 , u * 2 , u *

3 ) and corresponding solutions to the state system S * , V * , L * , I * , T * , that minimize objective functional (2.2) , there exists adjoint variables λ 1 , λ 2 , λ 3 , λ 4 and λ 5 , satisfying , u 1 max ,

Proof. The result follows from a direct application of a version of Pontryagin's Maximum Principle for bounded controls [36, 37] . The differential equations governing the adjoint variables are obtained by differentiation of Hamiltonian function (3.1) , evaluated at the optimal control. Then the adjoint system can be written as follows: (3.5) and = 1 , 2 , 3) .

We point out that the optimality system consists of the state system (1.2) with the initial conditions, adjoint system (3.2) with transversality conditions (3.3) , and optimality condition (3.4) . Thus, we have the following optimality system:

Due to the boundedness for the state system, the adjoint system has bounded coefficients which is linear in each of the adjoint variables. Therefore, the adjoint system have finite upper bounds. We will state the following proposition (without proof) needed for the proof of the uniqueness of the optimality system for the small time window.

Proposition 4.1. [13] The function u * (m ) = min { max { m, a } , b} which is Lipschitz continuous in m, with a < b for some fixed nonnegative constants. Now, we are in position to prove the uniqueness of the optimality system in the same way as shown in [12 , 13 , 17] . Proof. The proof is done by contradiction. Note that the boundedness for the state system and the adjoint system, the Lipschitz continuity of optimal control functions, and some elementary inequalities play key roles in obtaining a contradiction. , u 1 max ;

For ease of notation, we generally omit the dependence on time in the following except in the case that a specific time is intended. We consider the difference of the resulting equations for x i and x i , y i and y i , i = 1 , 2 , . . . , 5 . From the first equation of (3.7) , we get dx 1 dt

and d x 1 dt

It follows from the above two equations that

(4.1)

Multiplying both sides of (4.1) by (x 1 − x 1 ) and then integrating both sides from 0 to t f , one has

In order to simplify the right-hand expressions of (4.2) , we need some elementary inequalities. An elementary inequality

It follows that | a − b| ≤ 6 32(a 6 + b 6 ) and thus | a By the elementary inequality (a + b) 2 ≤ 2(a 2 + b 2 ) , we have

where C depends on bounds for x 1 , y 1 . Another inequality is also needed to obtain the specific estimate:

where M i (i = 1 , 2 , . . . , 8) and N j ( j = 1 , 2 , . . . , 9) depend on the coefficients and the bounds of the state variables and co-state variables. Table 1 . Table 1 .

We employed parameters values similar to those reported in previous studies whenever possible. Table 1 provides a summary of the parameter values used in numerical simulation of model (3.7) . Due to lack of data, we make hypothesis on some parameter values within realistic ranges. The initial values for the states are given by Table 1 . Table 1 . Table 1 . Table 1 . We start comparing the case of the population with optimal control actions and without control actions. In Figs. 1 -5 , the population with optimal control actions are shown in solid lines, compared with the population without control actions which are shown in dashed lines. We observe from Fig. 1 that the susceptible population is higher under control compared to the case without control. The results in Fig. 2 predict a higher level of vaccinated population for the aforementioned tuberculosis intervention strategies under control compared to the case without control. The results depicted in Figs. 3 and 4 clearly suggest that the optimal control results are very effective in control of individuals infected with TB in latent stage and in the active stage, as expected. In Fig. 5 , the population of treated humans declines sharply in the first 18 years and after that decrease further, while Fig. 4 shows slight increase in the last two years of the intervention strategy. This is due to the fact that the treatment control u 3 becomes zero after 18 years. The simulation results in Fig. 6 suggest that this strategy would require that the vaccination control u 1 should be at maximum for almost the entire period of intervention, while treatment controls u 2 and u 3 should maintain maximum effort s f or 13 years and 18 years respectively before decreasing to zero.

In what follows, we consider different control strategies in the reduction of infected individuals (including latent, infectious, being treated individuals). We observe in Fig. 7 that the infected individuals are lowest when three controls are considered. As the cost is one of the important parameter to decide for applicable strategy, we performed the cost design analysis and made comparative study of the cost for various strategies. The cost profiles for various control strategies are given in Fig. 8 . Note that the cost is least when three controls are considered, that is, the implementation of three control measures to restrict the epidemic while at the same time keeping the cost of vaccination and treatment very low are good policies for the achievement of our goal. In Fig. 9 we present the effect of varying the disease transmission coefficient β on the number of individuals infected with TB in the active stage ( I ( t )). We observe that after about 12 years the effect of β is negligible on individuals infected with TB in the active stage. Fig. 10 shows the effect of varying fraction coefficient l on the number of individuals infected with TB in the active stage ( I ( t )). It is easy to see that the individuals infected with TB in the active stage affected much by l .

In this work, we have studied a simple mathematical model of TB presented by Liu and Zhang [24] . In order to control the spread of TB, we applied a preventive control in the form of vaccination and two treatment controls to susceptible, latent, and individuals infected with TB in the active stage. Next, we established the objective functional by remembering the fundamental biomedical goals in our mind. Theoretically, we proved the existence of optimal control and studied the characterization of optimal control by Pontryagin's Maximum Principle. In order to assure that the optimal control is unique in the depictions, we investigated the uniqueness of the optimal control system.

Numerically, the simulations of the extended TB model (1.2) indicate that the vaccination and treatment strategies are effective in reducing TB disease transmission, especially we may get maximum disease control by using minimum cost as shown in Figs. 7 and 8 . It could therefore be concluded that the comprehensive effect of all the three controls not only minimizes cost burden but also keeps a tab on infective population.

Undoubtedly, dynamics of TB infection is far more complicated and varied than the one captured by this mathematical model. Therefore, we will incorporate the time delay for drug resistance into the present model so as to provide us a better assessment for designing TB dynamics. In addition, we note that the principle technique for optimal control problem is to solve a set of "necessary conditions" that an optimal control and corresponding state must satisfy. If the Hamiltonian is linear in the control variable u , it would be difficult to solve for u * from the optimality equation. We will leave these cases as our future work.

Tuberculosis control: past 10 years and future progress

Evolution of WHO, 1948-2001 policies for tuberculosis control

For the global surveillance and monitoring project: assessment of worldwide tuberculosis control

Exogenous reinfection in tuberculosis

A model for tuberculosis with exogenous reinfection, Theor

Optimal control of treatments in a two-strain tuberculosis model

Optimal control of vector-borne diseases: treatment and prevention

The intrinsic transmission dynamics of tuberculosis epidemic

Dynamical models of tuberculosis and their applications

Optimal control applied to vaccination and treatment strategies for various epidemiological models

The combined effects of optimal control in cancer remission

Vaccination and treatment as control interventions in an infectious disease model with their cost optimization

Optimal control of an influenza model with seasonal forcing and age-dependent transmission rates

Keeping options open:an optimal control model with trajectories that reach a DNSS point in positive time

An optimal strategy for HIV multitherapy

Control strategies for tuberculosis epidemics: new models for old problems

Mathematical models for the disease dynamics of tuberculosis

Modeling epidemics of multidrug-resistant m. tuberculosis of heterogeneous fitness

Prospects for worldwide tuberculosis control under the who dots strategy. directly observed shortcourse therapy

The reinfection threshold promotes variability in tuberculosis epidemiology and vaccine efficacy

The natural history of tuberculosis: the implications of age-dependent risks of disease and the role of reinfection

Global stability for a tuberculosis model

Rate of reinfection tuberculosis after successful treatment is higher than rate of new tuberculosis

Optimal interventions strategies for tuberculosis

Optimal intervention strategy for prevention tuberculosis using a smoking-tuberculosis model

Optimal control for a tuberculosis model with undected cases in cameroon

Modeling the impact of early therapy for latent tuberculosis patients and its optimal control analysis

Cost-effectiveness analysis of optimal control measures for tuberculosis

Optimal control applied to tuberculosis models. the IEA-EEF european congress of epidemiology 2012: epidemiology for a fair and healthy society

Optimal control for a tuberculosis model with reinfection and post-exposure interventions

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

Deterministic and Stochastic Optimal Control

Mathematical Theory of Optimal Processes

The authors are grateful to the editor and the referees for their valuable comments and suggestions.

3 (x 1 − x 1 ) 2 + (y 1 − y 1 ) 2 + (y 2 − y 2 ) 2 .(IV) case of x 1 (y 1 − y 2 ) < x 1 ( y 1 − y 2 ) < 0 . By aid of 3 An easy induction gives thatwhere K ≥ 1 + 2 d . In order to simplify the right-hand of (4.2) , another common expression can be used repeatedly,where C depends on bounds for x and y . Based on the above arguments, we findwhere M 1 , N 1 are appropriate upper-bounds. Similarly, we can obtain the following inequalities for (x i (t f ) , x i (t f )) and (y j (0) , y j (0)) with i = 2 , 3 , 4 , 5 and j = 1 , 2 , 3 , 4 , 5 :Now it follows from (4.5), (4.6), (4.7), (4.8), (4.9), (4.10), (4.11), (4.12), (4.13), (4.14) that 0 (x 1 − x 1 ) 2 dt will be nonnegative. Similar arguments apply to remaining integral terms, we can obtain all of the other λ s and t f s . Take the maximum of all of the λ s used as λ and the minimum of the t f s used as t f , then the coefficient of each integral term in (4.15) is nonnegative. This implies that Hence the solution of (3.7) is unique for small time. This ends the proof. Remark 4.1. We note that Theorem 4.1 shows that the solution of (3.7) is unique when t f is small. As pointed out by Lenhart and Workman [36] , we can usually obtain uniqueness of the solutions for (3.7) , but only for small time t f . This small time condition is due to reverse time orientations of the state equation and adjoint equation.

In this section, a fourth-order Runge-Kutta algorithm is used to solve optimal system (3.7) numerically. With initial conditions for the states, chosen an initial guess for the controls, the state differential equations (1.2) are solved forward in time by a fourth-order Runge-Kutta scheme. With the application of the current iteration solution of the state equations (1.2) and the transversality conditions (3.3) , the adjoint system (3.2) 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 (3.4) , and then the iterative process is repeated. This procedure is continued iteratively till current state, adjoint, and control values converge sufficiently.