Transmission dynamics model of Tuberculosis with optimal control strategies in Haramaya district, Ethiopia

In this study, we use a compartmental nonlinear deterministic mathematical model to investigate the effect of different optimal control strategies in controlling Tuberculosis (TB) disease transmission in the community. We employ stability theory of differential equations to investigate the qualitative behavior of the model by obtaining the basic reproduction number and determining the local stability conditions for the disease-free and endemic equilibria. We consider three control strategies representing distancing, case finding, and treatment efforts and numerically compare the levels of exposed and infectious populations with and without control strategies. The results suggest that combination of all controls is the best strategy to eradicate TB disease from the community at an optimal level with minimum cost of interventions.


Introduction
TB is a chronic infection of the bacteria caused by Mycobacterium TB, which was discovered by Koch in 1882. Different scholars classify TB in various ways based on the symptoms and treatments of the disease. The medical community divides TB into two categories, pulmonary (smear-positive and smear-negative) and extra-pulmonary. Pulmonary TB is a common form of TB that affects the lungs, whereas extra-pulmonary TB affects other parts of the body such as pleura, lymph nodes, back, joints, urinary tract of the genitor, nervous system, and abdomen [14]. A susceptible person is infected with TB when he or she inhales the TB germs that are released into the air while the infectious person with pulmonary TB disease coughs, sneezes, cries, or sings [6,14]. The first stage of infection is the latency period. At this stage an individual does not exhibit any symptoms of the disease and does not infect others. Such individual is said to be infected with latent TB. The second stage is the time of active TB infection when a person begins to display symptoms such as general weakness or fatigue, fever, weight loss, lack of appetite, and night sweating and can also infect others [19].
TB is an ancient disease, which caused more suffering and death than any other infectious disease and remains a significant global health issue. Even today, after advanced screening, diagnosis, and treatment methods have been established, one third of the world's population has been exposed and is infected with TB [18]. Although TB cases are currently reducing in many developed countries, they are rising in Africa, Eastern parts of Europe, and Asia [17]. The government of Ethiopia, in collaboration with development partners, bilateral and multilateral organizations, and wider community, is working together to combat this deadly disease to meet the sustainable development targets by 2030 and end TB by 2035. However, according to the World Health Organization (WHO) report 2019 [23], Ethiopia is among 30 high TB and TB/HIV burden countries in which 70% of notified cases are within the age group of 15-54 years. This indicates that TB is affecting mostly the productive working group so that it can bring back the economy of the country which demands extensive work to decide the right control strategies.
Mathematical modeling of infectious disease is used as a powerful tool for studying the dynamics of disease transmission and the effect of different intervention strategies to inform public health policy makers on the implementation of effective intervention programs to combat infections. Several researchers have built different mathematical models to investigate the dynamics of TB disease transmission and optimal control strategies; for example, see [5,15,16,20,21,24,25] and some references therein. Gao and Huang [12] analyzed TB model that incorporates vaccination of susceptible individuals, identification for treatment of latently infected individuals, and treatment of individuals with active TB. Their analysis showed that the combined implementation of three controls is most effective and less expensive among different strategies. In [2] an optimal TB control model was analyzed by considering the case detection of TB infections as optimal control parameter. Investigating the effectiveness of optimal control by comparing the levels of exposed and infectious populations with and without optimal control, the authors suggested that the optimal control approach would deliver better results in terms of reducing the number of infectives. Abouelkheir et al. [1] proposed a mathematical model for the optimal control strategy of a drug-resistant TB by considering a treatment program. The numerical simulations of their study showed that the implementation of a treatment as a control has a positive impact on the reduction of infectious individuals. Whang et al. [21] developed a dynamic TB transmission model in South Korea using a modified Susceptible-Exposed-Infectious-Recovered model with the time-dependent parameters and proposed optimal treatment strategies. Using the data of active-TB incidence obtained from 1970 to 2009, they estimated the parameters by the least-squares curve fitting. They considered three control mechanisms (distancing, case finding, and case holding efforts) and used optimal control programs to minimize the number of exposed and infectious individuals and the cost of implementing the control treatment. The results showed that the most effective control factor for preventing TB transmission is a distancing control. All of the above studies showed substantial results for TB disease transmission dynamics by taking situations in different countries into consideration. In this paper, we propose a dynamic model for TB disease transmission in Hramaya district, Ethiopia. The model is a modified and extended version of the model presented in [9] with optimal control strategies for the control of the TB disease.
The rest of the paper is organized as follows. In Sect. 2, we describe the basic model formulation. In Sect. 3, we present qualitative analysis of the basic model. In Sect. 4, we formulate and analytically study an optimal control problem using Pontryagin's maximum principle. In Sect. 5, we present the numerical results and discussions. Finally, Sect. 6 concludes the paper.

Model description and formulation
We formulate a mathematical model that describes the dynamics of TB infection in a population. The model divides individuals in the population under study into four compartments according to their epidemiological status. These are the susceptible individuals (S), who are not yet infected by the TB bacterium but with a possibility to be infected; the exposed individuals (E), who recently infected but are not infectious; the infectious individuals (I), who have become infected with TB and are able to transmit the disease; and the recovered individuals (R), who have recovered from the disease. We assume that the total population N remains constant (that is, natural birth rate and death rate are equal, and there is no emigration or immigration) and there are no TB disease-related deaths. So, the natural birth and death rate is μ. We also assume that an individual is firstly exposed to the disease through contacts with infectious individuals and does not become infectious instantly. That is, after the initial infection, an individual stays in a latent period for some time and then either becomes infectious (by developing the disease) or recovered (when immunological defense kills the inhaled bacilli). Hence the rate at which people become susceptible starts with rate μN . From that rate, the normal death rate of people that have not yet been exposed to TB, μS, and the rate at which susceptibles by meeting with infectious become exposed, βSI/N , are subtracted. From those exposed people, those who die, μE, those who become infectious after being exposed, kE, and those who become recovered, αE, are reduced. To describe the rate at which people in the population become infectious, from kE, the percentage of death among infectious, μI, and the percentage of recovery from the TB infectious, γ I, are subtracted. The last equation describes the rate at which people in the population recover from a disease; that is, individuals who do not progress to TB infectious αE are moved from class E to the recovered class R, and the recovered individuals γ I are also moved from class I to R. From that amount, the percentage of death among recovered people μR is subtracted. This model description is shown in Fig. 1. Based on this model flow diagram, the dynamics of the model is given by the following system of differential equations: with initial conditions The total population N(t) = N at time t is the sum of the populations in the compartments S, E, I, and R. Since the variable R does not appear in the first three equations of model (1), it is possible to examine model (1) by studying the subsystem and determining R from R = N -(S + E + I) or dR dt = αE + γ I -μR.

Model analysis
Since model (1) tracks human populations, all parameters associated with it are nonnegative.

Invariant region
From model (1) the addition of all the equations gives which indicates that the total population N is constant. Also, since model (1) represents interaction between individuals in the human population, it is plausible to assume that all the state variables of the model are nonnegative for all t ≥ 0 (see the proof of Theorem 3.1). Further, it can be shown that the region = (S, E, I, R) ∈ R 4 + : is positively invariant. Thus the solution of model (1) with initial conditions in remains there for t ≥ 0. Furthermore, the usual existence, uniqueness, and continuation results hold in for the equations of model (1) so that model (1) is well posed both mathematically and epidemiologically. So it is sufficient to study the dynamics of the basic model (1) in . Proof Given that the initial data S 0 , E 0 , I 0 , R 0 are nonnegative, it is clear from the first equation of model (1) that

Positivity of the solutions
Integrating both sides of this inequality and applying the technique of separation of variables along with initial condition yield Also from the second equation of model (1) we have Integrating both sides, we have Then solving by technique of separation of variables and applying the initial condition, we get Applying similar steps to the third and fourth equations of model (1), we obtain, respectively, and This completes the proof of the theorem.
Therefore the solution of model (1) is positive. Since studying model (1) is similar to studying the reduced model (3), we focus the following discussions on model (3).

Equilibria and the basic reproduction number
The system of model (3) has two nonnegative equilibrium points called the disease-free equilibrium (DFE), where E = I = 0, and endemic equilibrium (EE), where E = 0 and I = 0. To obtain the DFE and EE points, the right-hand side of model (3) is set to zero and then solved for the values of S, E, and I. Then the DFE point of the respective states becomes and the EE point becomes To calculate the basic reproduction number R 0 of model (3), we use the method of the next generation matrix given in [10] and obtain

Local stability of the equilibria
For the stability analysis of the DFE and EE points, the linearization matrix (Jacobian matrix) of model (3) is evaluated at each equilibrium point, and then its eigenvalues are determined by solving its characteristic equation. Based on the eigenvalues, the equilibrium points will either be stable (if all the eigenvalues of the Jacobian evaluated at the equilibrium point contain negative real parts) or unstable (if at least one of the eigenvalues of the Jacobian evaluated at the equilibrium point has a positive real part). (3) is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.

Theorem 3.2 The DFE of model
Proof The Jacobian matrix of model (3) is For the DFE point E 1 of Eq. (10), The eigenvalues of the Jacobian matrix J are the solutions of the characteristic equation where I is the identity matrix, that is, This leads to the characteristic equation given by where Based on the Routh-Hurwitz stability criteria, without solving Eq. (17), we get that the sufficient and necessary conditions for stability are The first inequality is clearly satisfied since all model parameters are positive, and the second inequality a 3 > 0 holds if R 0 < 1. For the third condition, Therefore by the Routh-Hurwitz criteria all eigenvalues of J(E 1 ) are negative when R 0 < 1, so that E 1 is locally asymptotically stable.
Then its characteristic equation is given by Clearly, b 1 > 0, and if R 0 > 1, then b 3 > 0, and Therefore by the Routh-Hurwitz criteria, all eigenvalues of J(E 2 ) are negative when R 0 > 1, so that E 2 is locally asymptotically stable.

Sensitivity analysis of the model parameters
In the study of biological dynamics the transmission dynamics of infectious disease model sensitivity analysis plays an important role. Using sensitivity analysis, the role of each parameter used in the model can be investigated, and a strategy can easily be developed to control the spread of infection in the community.
In this section, the normalized forward sensitivity index of R 0 to a model parameter, which is a ratio of the relative change in R 0 to the relative change in the parameter defined by where p represents all the model parameters (μ, β, k, α, γ ), is carried out in the sense of [3]. Using (23), the sensitivity indices of R 0 with respect to the model parameters are presented in Table 1 using the values μ = 0.015267, β = 10, k = 0.0444, α = 0.2325, and γ = 1.25 taken from Sect. 5. The positive sign of sensitivity index of R 0 to the model parameters indicates that an increase (or decrease) in the value of each of the parameter leads to increase (or decrease) in R 0 , whereas the negative sign indicates that an increase (or decrease) in the value of parameter leads to a corresponding decrease (or increase) in R 0 of model (3). From Table 1, β and k have the impact of expanding the disease in the community if their values are increased, whereas μ, α, and γ have an influence of minimizing the burden of the disease in the community as their values increase.

Optimal control problem
Optimal control is the method of determining control and state trajectories for a dynamic system over a period of time to minimize a performance index or an objective functional.
In this section, three time-dependent controls (distancing, case finding, and treatment) are incorporated to the model of Eq. (1). The distancing control u 1 (t) is the effort of preventing susceptible individuals from becoming infectious using strategies such as early detection and isolation of infectious individuals, and health education campaign (prevention of careless spitting, coughing, and sneezing). It is incorporated by adding a control term that characterizes the contact between susceptible and infectious individuals so that the rate of infection will be reduced. The case finding control u 2 (t) is the effort of identifying TB exposed individuals through screening and put under treatment to reduce number of individuals that may progress to TB infectious. The treatment control u 3 (t) is the effort of ensuring for those that are infectious treatment and monitoring in taking their drugs to reduce the number of individuals developing and dying of TB. Then, after incorporating the controls into the model described by Eq. (1), we get the following optimal control model: with the initial conditions of Eq. (2). Here our purpose is finding the optimal values of u 1 , u 2 , and u 3 so that the state and control trajectories minimize the objective functional given by subject to the state system (24). In Eq. (25), t f is a fixed final time, whereas the coefficients B 1 , B 2 , and B 3 are positive weight constants, which balance the cost factors associated with control measures u 1 , u 2 , and u 3 , respectively. The cost of each control measure is assumed to be nonlinear and take quadratic form, that is, B 1 2 u 2 1 , B 2 2 u 2 2 , and B 3 2 u 2 3 are the costs of the control measure associated with u 1 , u 2 , and u 3 , respectively. Thus we seek the optimal controls, u * where U = {(u 1 , u 2 , u 3 ) : each u i , i = 1, 2, 3 is Lebesgue measurable with 0 ≤ u i ≤ u i max for 0 ≤ t ≤ t f } is the set of acceptable controls.
Proof The nontrivial criteria regarding the set of admissible controls U and the set of end conditions from the theorem of Fleming and Rishel [11] are described and checked below.
A. The set of all solutions to Eqs. (24)-(26) with corresponding control functions in U is nonempty.
B. The state system (24) can be written as a linear function of the control variables with coefficients dependent on time and the state variables.
C. The integrand L of (25), L(t, S, E, I, R, u 1 , u 2 , u 3 ) = E + I + B 1 2 u 2 1 + B 2 2 u 2 2 + B 3 2 u 2 3 is convex on U and additionally satisfies L(t, S, E, I, R, u 1 where d 1 > 0 and τ > 1. To establish condition A, Picard-Lindelöf 's theorem from [3,4] is referred. If the solutions of the state equations are bounded and the state equations are Lipschitz continuous in the state variables, then there is a unique solution corresponding to every admissible control U. With the boundedness of the solutions of model (1) already justified, it follows that the state system is continuous and bounded. It is also equally direct to show the boundedness of the partial derivatives with respect to the state variables in the state system, which confirms that the system is Lipschitz in relation to the state variables [7]. This completes the proof of condition A. Condition B is confirmed by observing linear dependence on controls u 1 , u 2 , and u 3 of the state equations. Finally, to verify condition C, by definition from [8] any constant, linear, and quadratic functions are convex. Therefore L(t, S, E, I, R, u 1 , u 2 , u 3 ) is convex on U. To prove the bound on L, note that by the definition of U we have Therefore L(t, S, E, I, R, u 1 ,

Characterization of the optimal controls
Necessary conditions that optimal solutions need to satisfy are obtained from Pontryagin's maximum principle. This principle converts Eqs.

Theorem 4.2
There exist optimal controls u * 1 , u * 2 , and u * 3 and the corresponding state solutions S * , E * , I * , and R * that minimize J(u 1 , u 2 , u 3 ) over U, and therefore there exist adjoint variables λ 1 , λ 2 , λ 3 , and λ 4 such that with transversality conditions Furthermore, the characterization of the optimal controls is given by Proof The adjoint system is obtained by taking the negative partial derivatives of H with respect to each state variable [13], evaluated at optimal controls and corresponding state variables as follows: Also, for characterization of the optimal controls, we first obtain the optimality conditions by taking the partial derivatives of H with respect to each control u i , i = 1, 2, 3, and setting them to zero, that is, Then, using the property of the control bounds 0 ≤ u i ≤ u i max , i = 1, 2, 3, the controls are given as This can be written in compact form as follows: We point out that the optimal controls and the corresponding state solutions are found by solving the following optimality system, which consists of the state system (24) with its initial conditions coupled with the adjoint system (28) and its transversality conditions (29) together with the characterization of the optimal controls (30): λ i (t f ) = 0, i = 1, 2, 3, 4. (31)

Uniqueness of the optimality system
Since the state and adjoint variables are bounded (because the adjoint system (28) is also linear in λ i , i = 1, 2, 3, 4, as the state system) and satisfy the Lipschitz condition rela-

Numerical results and discussions
In this section, we first estimate the model parameters. Second, we analyze the stability of the equilibrium points. Then we investigate the effect of the controls in reducing the infected population with TB under different strategies.

Estimation of model parameters
Since some parameter values such as demographic data are well known from the literature, we have estimated unknown parameters based on the data obtained from the health office of the Haramaya district and Haramaya hospital to keep the model more realistic. The natural mortality μ is postulated to be equal to the inverse of the life expectancy at birth, which is now about 65.5 years in Ethiopia [22]. Therefore the natural death rate μ in Ethiopia is 1 65.5 = 0.015267 per year. The mean infectious period 1 γ is estimated in the range of 0.5-2 years. Here we choose γ = 0.5+2 2 = 1.25 year -1 . The parameter β is in the range of 10-15 year -1 by the assumption that a person with active TB infects an average of 10-15 other people every year [18]. So, we have taken β = 10 year -1 . The parameters k and α, which cannot be obtained by the observed data or references, are estimated using the least-squares data fitting method with the help of a MATLAB tool fminsearch, which is part of the optimization toolbox. This is done by fitting model (1) to the observed data of active TB incidence in 2010-2019 ensuring that the minimum error occurred whilst getting the best fits. The fitted curve to the observed data is shown in Fig. 2 and the estimated values of k and α are presented in Table 2. To quantify parameter uncertainty and construct confidence intervals (CI), we rely on the parametric bootstrapping method. Bootstrapping is a statistical method for assigning accuracy measures to the sample estimates [6]. In this method, multiple observations are repeatedly sampled from the best-fit model to quantify parameter uncertainty by assuming that the time series follow a Poisson distribution centered on the mean at the time points.
In addition to the parameters, we take the following as an initial population of the district relating to the year 2010, and for numerical results, we divide the total population N = 361,000, as follows. We assume that 80% of the population are susceptible [18], whereas the number I 0 of infectious individuals in 2010 was estimated by the number of smear positive patients obtained from the health office of the Haramaya district and Haramaya hospital. So, with this assumption, we have the following initial data for the simulation purpose: (S 0 , E 0 , I 0 , R 0 ) = (288,800, 3455, 155, 68,590).

Stability analysis of equilibrium points
Based on the above real data, we found from Eq. (12) that R 0 = 1.2011 > 1. This in principle shows that TB does spread in the community of the district. Since R 0 > 1, the prevalence of TB will result in an epidemic. This is due to the fact that the rate of transmission is greater than the recovery rate. From Eq. (10), the DFE point is determined to be E 1 = (361,000, 0, 0), and from Eq. (11), EE point E 2 = ( N R 0 , μN(μ+γ ) βk (R 0 -1), μN β (R 0 -1)) = (300557. 8, 3158.4, 110.8).
The characteristic equation from (17) becomes λ 3 + 1.5727λ 2 -0.0506λ -0.0011 = 0, which shows that the DFE point is unstable by the Routh-Hurwitz criterion. This means that when an individual infected with Mycobacterium TB is present in a susceptible population, it will eventually result in an outbreak of the disease. The characteristic equation from (21) also becomes λ 3 + 1.5758λ 2 + 0.0286λ + 0.0011 = 0, which shows that the EE point is stable. This means that TB in the study area will persist.

Optimal control problem analysis
In this section, we numerically analyze the impact of control strategies on controlling transmission of the TB disease using the Matlab R2019a program. We solve the optimality system (31), which consists of eight ordinary differential equations from state and adjoint equations, using the forward-backward sweep method developed by Lenhart and Workman [13]. The numerical procedure starts with an initial guess on the control variables. Then, using fourth-order Runge-Kutta scheme, we solve the state equations simultaneously forward in time starting from the initial conditions, and using the state solutions, we iii. Strategy C: Applying treatment (u 3 ) only. iv. Strategy D: Combination of use of distancing (u 1 ) and case finding (u 2 ). v. Strategy E: Combination of use of distancing (u 1 ) and treatment (u 3 ). vi. Strategy F: Combination of use of case finding (u 2 ) and treatment (u 3 ). vii. Strategy G: Combination of use of distancing (u 1 ), case finding (u 2 ), and treatment (u 3 ).

Strategy A
Here we only use distancing control u 1 to optimize the objective functional J in (25). From Fig. 3 we observe that the strategy has a significant effect in reducing the number of exposed and infectious populations when compared to the numbers in the case without controls.

Strategy B
Under this strategy, we use the case finding control u 2 to optimize the objective functional J in (25) while setting other controls to zero. From Fig. 4 we observe that this strategy is effective in reducing the number of exposed and infectious populations.

Strategy C
In this strategy, we only use treatment control u 3 to optimize the objective functional J in (25) while setting other controls to zero. From Fig. 5 we observe that the strategy is effective in reducing the number of exposed and infectious populations. From the numerical computation we get that the numbers of infectives with strategies A, B, and C are 306, 701, and 1292, respectively, at final time. In Figs. 6(a), 6(b), and 6(c), the control profiles of strategies A, B, and C are displayed, which suggest that controls u 1 , u 2 , and u 3 of the respective strategy are at upper bound for period of 8.1, 8.3, and 7.7 years, respectively, before sharply drop to zero at final time. These clearly indicate that the number of infection is reduced by a larger amount in strategy A than in strategy B and C, whereas strategy C requires less effort than strategies A and B.

Strategy D
The distancing control u 1 and case finding control u 2 are applied under this strategy, whereas treatment control u 3 is set to zero. We observe from Fig. 7 that the number of exposed and infectious populations decrease when compared to the numbers in the case without controls.

Strategy E
Here both distancing control u 1 and treatment control u 3 are used, whereas case finding control u 2 is set to zero. We observe from Fig. 8 that the strategy is effective in reducing the number of exposed and infectious populations.

Strategy F
In this strategy, the objective functional J is optimized using case finding control u 2 and treatment control u 3 . We observe from Fig. 9 that this strategy shows a significant effect in reducing the number of exposed and infectious populations.

Strategy G
Here we used all the three control strategies u 1 , u 2 , and u 3 to optimize the objective functional J. We observe from Fig. 10 that the number of exposed and infectious populations highly decrease when compared to the numbers in the case without controls.
The results from the numerical computations show that the numbers of infectives with strategies D, E, F, and G are 89, 305, 259, and 88, respectively, at final time, which indicates that strategy G is more effective in reducing the number of infection than other strategies. In Figs. 11(a), 11(b), 11(c), and 11(d) the control profiles of strategy D, E, F, and G are shown, which indicate that strategy G requires less effort than other strategies.

Figure 12 Cost function of the intervention strategies
From the above results, all strategies are effective in reducing the number of exposed and infectious populations when compared with a situation without controls. The strategies with only one control are less effective in reducing the number of exposed and infectious populations when compared to strategies that incorporate more than one control. Also, even if the strategies with only one control are less effective, they are cheaper than the strategies that incorporate more than one control, as indicated in Fig. 12.

Conclusion
In this paper, we developed a dynamic model for TB transmission in the Haramaya district of Ethiopia based on active-TB incidence data recorded by the district health office and the Haramaya hospital. The model qualitative analysis demonstrates that model solutions are bounded and positive. The basic reproduction number and local stability conditions of the equilibria were determined, and the basic reproduction number was used to perform the sensitivity analysis that identified the role of each parameter used in the model to the spread of TB disease. Three time-dependent control measures such as distancing, case finding through screening of TB exposed individuals, and treatment of TB infectious population were applied to the model so that basic model was extended into an optimal control model. Theoretically, we proved the existence of an optimal control and obtained necessary conditions that optimal solutions need to satisfy using Pontryagin's maximum principle. Numerically, the levels of exposed and infectious populations with and without optimal controls were compared under different intervention strategies to investigate the roles of control measures independently and as combination in eradication of the TB disease. Implementing strategies with only one control is less effective in reducing the number of exposed and infectious populations when compared with strategies that incorporate more than one control. From strategies that incorporate more than one control measure, the strategy with combination of all three control measures is the best strategy in reducing the number of TB infection with minimum cost of interventions.