Optimal control analysis of a tuberculosis model.

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.

Introduction

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 . 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 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 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 where ρ 2 < 1 represents the reduction in risk of infection owing to vaccination (for a vaccine that offers 100% protection, ; thus in reality 0 < ρ 2 < 1). For model (1.1), Liu and Zhang [24] first obtained the basic reproduction number where 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, for system (1.1), ρ represents successful treatment rate. Computing the partial derivative of R 0 with respect to ρ gives that if . Therefore, improving successful treatment rate helps to reduce TB infection when . Using the fact that as long as Liu and Zhang [24] pointed out that the treatment of TB-infected individuals in the active stage would have positive impact if . In other words, if ρ 1 lies between two certain thresholds, then increasing ρ and γ have positive meaning on controlling TB infection. Furthermore, 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. ) between cost and efforts 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: 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.

The control problem

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, the problem is to minimize the objective functional

Here, the total cost on a finite time horizon [0, t] consists of the cost induced by the disease itself and the cost induced by vaccination and treatment efforts. 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 efforts such as regularity of drug intake) is taken as . Similarly, the cost incurred in providing treatment is taken as . The coefficients, 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 and . We seek optimal controls u 1, u 2 and u 3 in U such that where 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.

The necessary and sufficient conditions of optimal control

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.

The existence of an optimal control

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 i.e. while all other variables are nonnegative. Then, from the S equation, we have . Thus, S(t) ≥ 0 for all t > 0. Similarly, it can be shown that V(t) ≥ 0, L(t) ≥ 0, I(t) ≥ 0 and T(t) ≥ 0 for all t > 0. 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 which implies that provided that otherwise, N is bounded above by N(0). Thus, N is bounded above by . The upper bound for N is also the upper bound for S, V, L, I and T since they are all nonnegative. The remainder of this section will be devoted to the proof of existence of optimal control on the premise that the solutions to each state variables is bounded. For convenience we let and denote the vector of model states and controls, respectively.

There exists an optimal control to problem (2.2) .

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

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

The admissible set of all solutions to system (1.2) with corresponding control in U is nonempty;

;

The control set is closed, convex and compact;

The integrand of the objective functional is convex in U.

In order to verify above conditions, we write Then it is easy to see that is of class C 1 and . Moreover, one has where and 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, Thus, condition (3) is satisfied.

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 we have to prove that where and are two control vectors with q ∈ [0, 1]. It follows that and Furthermore, we have and so the proof is complete. □

It is easy to see that Theorem 3.1 still holds if or are replaced by convex functions of .

Characterization of an optimal control

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, 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 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 with transversality conditions Furthermore, we may characterize the optimal pair by the piecewise continuous functions

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: and evaluated at the optimal control and corresponding states, which results in adjoint system (3.2) with transversality conditions (3.3).

Finally, the optimality conditions require that, for the optimal pair on the interior of the control set, and this condition is simplified in Eqs. (3.4) with special attention on control arguments involving the bounds on the controls as defined with U in Eq. (3.2). This completes the proof. □

It is easy to see that Theorem 3.2 still holds if or are replaced by differential convex functions of .

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: where

Uniqueness of the 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.

The function 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 , , .

For sufficiently small t (3.7) is unique.

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. Now we proceed to prove this theorem. Suppose (S, V, L, I, T, λ 1, λ 2, λ 3, λ 4, λ 5) and are two non-identical solutions to the optimal system (3.7). To show that the two solutions are equivalent, it is convenient to make a change of variables.

Let where λ is a positive constant to be chosen later. With the new variables the optimality conditions become 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 and y and . From the first equation of (3.7), we get and It follows from the above two equations that Multiplying both sides of (4.1) by and then integrating both sides from 0 to t, one has In order to simplify the right-hand expressions of (4.2), we need some elementary inequalities. An elementary inequality shows that It follows that and thus By the elementary inequality we have where C depends on bounds for .

Another inequality is also needed to obtain the specific estimate: If 0 ≤ b ≤ a, then . Using Proposition 4.1 and the above inequalities, we have the following estimation of . We divide our estimation in five cases as follows.

(I) case of . It follows from (4.3) that where M depends on bounds for .

(II) case of . By the inequality and (4.4), one has

(III) case of . By means of and (4.4), we get

(IV) case of . By aid of and (4.4), we have

(V) case of . Using the inequality and (4.4) again, we obtain An easy induction gives that Set and Then where .

In order to simplify the right-hand of (4.2), another common expression can be used repeatedly, where depends on bounds for and y.

Based on the above arguments, we find where M 1, N 1 are appropriate upper-bounds. Similarly, we can obtain the following inequalities for and with and : where and depend on the coefficients and the bounds of the state variables and co-state variables.

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 Notice that the coefficients of all of integrals in (4.15) are nonnegative if we choose a sufficiently large λ and a sufficiently small t. For example, if we fix and then the coefficient for the integral will be nonnegative. Similar arguments apply to remaining integral terms, we can obtain all of the other λ and t. Take the maximum of all of the λ used as λ and the minimum of the t used as t, then the coefficient of each integral term in (4.15) is nonnegative. This implies that and Hence the solution of (3.7) is unique for small time. This ends the proof. □

We note that Theorem 4.1 shows that the solution of (3.7) is unique when t 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. This small time condition is due to reverse time orientations of the state equation and adjoint equation.

Numerical results and discussion

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.

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 respectively. The weights in the objective function are and which means that minimization of the number of latent, infectious, and treated has different importance. We choose optimal control for a period of 20 years in keeping with what is in other literature on treatment of tuberculosis such as [34] with references.

Table 1

Parameter values used for model simulations.

Parameter	Value	Unit	Source
Λ	1428	Person year ^−1	[34]
β	0.003	Person ^−1 year ^−1	[34]
ρ1	0.25	None	[32]
µ	170	Year ^−1	[11]
ρ2	0.3	None	Assumed
l	0.9	None	[32]
δ	0.00368	Year ^−1	[10]
α	0.17	Year ^−1	[2]
u1max , u3max	1	None	[14]
u2max	1.125	Year ^−1	[34]

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 efforts for 13 years and 18 years respectively before decreasing to zero.

Fig. 1

Dashed line: susceptible population without optimal control. Solid line: susceptible population with optimal control. Weight constants are and other parameter values are presented in Table 1.

Fig. 5

Dashed lines: treated population without optimal control. Solid lines: treated population with optimal control. Weight constants are and other parameter values are presented in Table 1.

Fig. 2

Dashed line: vaccinated population without optimal control. Solid line: vaccinated population with optimal control. Weight constants are and other parameter values are presented in Table 1.

Fig. 3

Dashed line: latent population without optimal control. Solid line: latent population with optimal control. Weight constants are and other parameter values are presented in Table 1.

Fig. 4

Dashed line: infectious population without optimal control. Solid line: infectious population with optimal control. Weight constants are and other parameter values are presented in Table 1.

Fig. 6

Three optimal control strategies profiles during 20 years. Weight constants are and other parameter values are presented in Table 1.

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.

Fig. 7

Infected individuals (including latent, infectious, being treated individuals) profiles under various control strategies during 20 years. Other parameter values are the same as those in Figs. 1–6.

Fig. 8

Cost profiles under various control strategies during 20 years. Other parameter values are the same as those in Figs. 1–6.

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.

Fig. 9

Simulations of TB model (1.2) showing the effect of varying transmission rates. Other parameter values are the same as those in Figs. 1–6.

Fig. 10

Simulations of TB model (1.2) showing the effect of varying fraction coefficients. Other parameter values are the same as those in Figs. 1–6.

Conclusion

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.