Mathematical Modeling and Simulation to Control the Spread of Multidrug-Resistant Tuberculosis

—Tuberculosis is an infectious disease caused by Mycobacterium tuberculosis. Tuberculosis that fails treatment will develop into multidrug-resistant tuberculosis. Research on the TB epidemic continues, particularly in the ﬁeld of applied mathematics with modeling. In this study, we analyzed a suitable strategy in controlling the development of susceptible individuals to active tuberculosis and even multidrug-resistant tuberculosis. In this work, local stability analysis was carried out around the equilibrium point. Also, to see the most inﬂuential parameters in the epidemic, a sensitivity analysis was performed on basic reproductive factors. Besides, the ﬁnal work was to do numerical simulations with some cases, so that the model could describe the disease’s phenomena and characteristics.


I. INTRODUCTION
T UBERCULOSIS (TB) is an infectious disease caused by Mycobacterium tuberculosis (Mtb). This disease is spread through the patient's saliva splashes. Irregular TB treatment causes TB to develop into drug-resistant TB or multidrugresistant TB (MDR-TB) and even rifampicin-resistant (RR-TB). WHO has published a global TB report every year since 1997. TB causes death in about a quarter of the world's population and 5-10% of people with latent TB can develop active TB disease during their life if the immune system declines. However, the potential for the development of active TB is higher among people with HIV.
In 2018, WHO announced around 1.2 million (range, 1.1−1.3 million) TB deaths among HIV negative patients and around 251, 000 deaths (range, 223, 000 − 281, 000) among HIV positive patients. Meanwhile, the total case of MDR / RR-TB globally is 186, 772. Moreover, Indonesia is one of the countries with a high MDR-TB burden in the world, around 24 thousand MDR / RR-TB cases have increased compared to the previous year [1]. The vaccine used to fight TB is Bacille Calmette -Guérin (BCG), developed in 1921. At present, an experimental vaccine is designed to stop latent TB from becoming TB disease. Recently, a study indicates that administering vaccine with M 72/AS01 E to susceptible adults can provide immunity and protect against the progression of pulmonary TB disease for at least three years [2]. Therefore, the development of vaccines that can limit initial infection and reactivation of latent TB needs to be done to control TB [3].
S. Suddin is with the Faculty of Education Sciences, Mathematics Education, Universitas Timor, East Nusa Tenggara, Indonesia e-mail: sulasri.suddin@gmail.com. E. N. Bano is with the Faculty of Science and Technology, Mathematics, Universitas Timor, East Nusa Tenggara, Indonesia e-mail: iranaik-teas@gmail.com.
The research development of TB is carried out to achieve sustainable development goals and final TB strategies. Research on the dynamics of TB has been carried out by Liu and Zhang [4] concerning the effects of vaccination and treatment. Furthermore, about vaccination, Suddin [5] developed a research related to endogenous reactivation and exogenous reinfection of TB spread. On the other hand, Ronoh et al. [6] developed a deterministic model for MDR-TB, which subsequently by Nagar et al. [7] considered first and secondline treatment in populations infected with MDR-TB. Then, by paying attention to treatment, Mishra and Srivasta [8] constructed MDR-TB models in the presence of quarantine and vaccination in vulnerable individuals.
However, work on strategies for eradicating the MDR-TB outbreak using a combination of vaccine and treatment continues to be very limited. Therefore, the researchers intend to carry out research related to controlling MDR-TB outbreaks. One of them is that vaccination plays an important role in the spread of TB and also considers other things, such as the presence of endogenous reactivation in latent TB individuals. Based on the information described, the researcher is interested in conducting similar research by applying it to the research environment, which is the Kefamenanu City sub-district area.
This research is expected to provide a strong scientific basis for determining the rate of vaccination in controlling MDR-TB. This research is important as an information to support the development of multidisciplinary science, mathematics, and medical science regarding the spread of MDR-TB. In addition, it can be used as a reference and comparison for health services in Kefamenanu regarding MDR-TB transmission.
In this work, the paper is structured as follows. Section II presents a formulation of a mathematical model for controlling MDR-TB epidemics and equilibrium stability. Section III presents sensitivity analysis. Section IV presents numerical simulations. The conclusions section is presented in Section V.

A. Formulation of the Epidemic Model
In the epidemic model, the researchers assumed that the effect of vaccination for both the susceptible and exposed humans was reinforced by the presence of a new vaccine being developed for people. The compartment model is divided into six classes namely susceptible human S(t), exposed human E(t), infected human I(t), resistant infected human R md (t), recovered human R(t), and vaccinated human V (t) at the time t, respectively. The parameters used in this model are the recruitment in the population by births with rate δ, natural death rate of human with rate µ, disease-induced deaths in infected TB human and MDR-TB with rate d 1 and d 2 respectively, the susceptible became exposed with rate β, endogenous reactivation with rate , recovery rate of infected human from infected TB and MDR-TB after treatment with k and α respectively, resistance rate of treatment γ, the vaccination rate coefficient for the susceptible human and exposed human ρ and ψ, respectively.
Based on the facts and assumptions, the following system of differential equations was obtained which is described in Fig. 1.

C. Computation of Basic Reproduction Number
Reproduction number, denoted by R 0 , is a key quantity, i.e. very important role, in epidemiological models. It can represent the average number of new infections by an infected individual. It is useful to predict how strong a disease outbreak and can be used to guide and evaluate the control strategies for the disease.
The researchers use the next generation matrix [9] for the model to set the reproduction number. The next generation matrix is divided in two states, infected (F ) and non-infected (V ). F is a matrix which contains new infection rate and V is a matrix which contains the transfer of individuals inside and outside of the infectious compartment, not new infection). Then, the spectral radius of F V −1 is equal to R 0 . From system (1), we define . Therefore, where

D. Local Stability Analysis
In this subsection, we will determine the stability of the equilibrium points by linearizing the system of differential equations (1) by obtaining its Jacobian at the equilibrium point.
By Jacobian matrix of the system (1) at disease-free equilibrium e 0 , we have (1) is locally asymptotically stable at disease-free equilibrium e 0 . If R 0 = 1 then e 0 is stable. If R 0 > 1 then e 0 is unstable.
Proof: We have following the characteristic equation of the Jacobian Matrix at e 0 , . We observed form above that all the eigenvalues are negative except λ 5,6 . Clearly, λ 5,6 < 0 if only if R 0 < 1. Hence, the system (1) at the disease-free equilibrium is locally asymptotically stable.
The endemic equilibrium point e * can be expressed in the following R 0 form: The Jacobian matrix of the system (1) at endemic equilibrium e * is given by (1) is locally asymptotically stable at endemic equilibrium e * . If R 0 = 1 then e * is stable.
Proof: The characteristic equation of the Jacobian matrix is as follows From (4), we have three negative eigenvalues, i.e. λ 1 * ,2 * = −µ, λ 3 * = −c 3 . While, the other eigenvalues are obtained by solving the following equation

Parameter
Value Source  Based on the Routh-Hurwitz Criterion, all the eigenvalues have negative real parts. Hence, system (1) is locally asymptotically stable at endemic equilibrium.
We illustrate the analytical results of the model by carrying out numerical simulations using a set of estimates for parameter values obtained from the literature. Model simulation of the dynamic on MDR-TB spread by using Maple and Matlab software. The parameters used in the epidemic model are presented in Table I.
The average number of populations per year based on data from RSUD Kefamenanu. For situations that do not allow the data will be taken proportionally based on the existing data. The following is a numerical simulation in the presence of vaccination in susceptible and latent TB individuals, where we divide into four cases and presented in four graphs. This simulation uses preliminary data in Kefamenanu City District. for (ψ = ρ = 0.01). Both of these conditions indicate that the disease is endemic in the population. Initially, the curve of the infected population has increased (peak). Especially, vaccine administration rates were only given to latent TB  individuals, the peak of curve I is slightly higher than curve R md (Figure 2). Whereas the rate of vaccine administration in vulnerable individuals and individuals with latent TB is low, the peak of curve R md is higher than curve I (Figure 3). From the two treatments in Figures 2 and 3, we can see that low vaccination interventions in both populations, the period needed for a decrease is longer than administering a vaccine to exposed individuals. This means that the disease will disappear or no epidemic occurs in the population. In Figures 4 and 5, we can see that the curve of infected increases (peak) at an early time. Further, at a certain time, all curves continue to decrease and converge to a constant value. For the two conditions, we can describe that the higher vaccination intervention is given to only vulnerable, the time needed to eliminate the disease is slightly longer than high intervention on both populations (Compare Fig. 4 and Fig. 5).
To describe the relationship between giving vaccines to susceptible individuals (ρ) and latent TB sufferers (ψ), the following is presented in the numerical simulation. Figure 6 is showing that the infection will persist and lead to an epidemic if no interventions or giving higher vaccination on latent TB people. Conversely, if given a high vaccination in both populations, or giving only to vulnerable individuals, TB disease will disappear from the population.

III. SENSITIVITY ANALYSIS
In this section, we will conduct a sensitivity analysis to find out which parameters are the most influential in the disease epidemic. In this case, we focus on the sensitivity of the parameters β, δ, µ, d 1 , , γ, ω, ρ, and k to the basic reproduction number (3).
To this aim, p and V denote the generic parameters and variables to be sought from (1), we evaluate the normalized sensitivity index [10] which is defined as follows.
The sensitivity index of the parameters β, δ, µ, d 1 , , γ, ψ, ρ, and k with R 0 is obtained using (6) with considering the parameter values in Table I. The results of the sensitivity analysis are presented in Table II for the order of parameters the most influential in the spread of MDR-TB disease.
The result of the sensitivity analysis shows that the contact rate (β), the birth rate (δ), and the natural death rate (µ) were the three parameters that most clearly influenced the change of reproduction number (R 0 ). Then it is followed by the recovery rate (k) in individuals with TB disease and progression rate form TB disease to MDR-TB (γ). The parameters β, δ, and have a positive sensitivity index. It means that if these  parameters value increases, the value of R 0 will also increase. Conversely, if the value of these parameters decreases, then the value of R 0 will also decrease. while the parameters µ, k, γ, ρ, ψ, and d 1 have a negative sensitivity index. It means that if these parameters value increases, the value of R 0 will decrease and vice versa it will increase. An increase or decrease of R 0 is obtained based on the magnitude of the sensitivity index.
If the parameters β and δ are increased by 10%, then the basic reproduction number R 0 increased by approximately 10% starting at 2.091357470 to 2.300493216. Conversely, if the parameters β and δ are decreased, then the basic reproduction number R 0 decreased by 11.11111% or 1.882221723. Moreover, if the parameters µ, k, and γ are increased, then the basic reproduction number R 0 decreased by 7.626922%, 5%, and 4.70104%, respectively. Conversely, if these parameters are decreased, then the basic reproduction number is increased.
Based on the sensitivity analysis conducted, we present numerical simulation of the first three variables play an important role in the spread of the epidemic, i.e. β, δ, and µ. Figure 7 presents several conditions. The first and third graph (low natural mortality) show that the effect of increasing contact rate or birth rate cause the disease to persist in the population. Similarly, we can see in the second graph that an epidemic occurs if birth rates and transmission rates continue to increase. Otherwise, contact rates and birth rates are low, the disease will disappear. This provides the possibility for health policymakers to reduce birth rates and contact rates in areas with high rates of TB transmission.
Hereinafter, the following graph is given the interaction between the parameters of contact rate, successful treatment in patients with active TB, and patients who fail to take medication and develop into MDR-TB, i.e. β, k, and γ, respectively. Figure 8 describes that there are three conditions. Based on these simulations, high contact rates of TB disease cause epidemics even though the rate of progression of active TB to MDR-TB is small. In contrast, although more and more active TB develops into MDR-TB, the disease will disappear from the population if the contact rate is small (see the first graph). The second case describes that a high contact rate will contribute to the spread of the disease, even though the treatment rate in patients with active TB is also high. However, if the treatment rates is high and contact rates is small, then no epidemic occurs (see the second graph). The third condition indicates that even if the rate of progression from active TB to MDR-TB is low, there will still be an epidemic as a treatment in TB disease is getting smaller. Vice versa, there will be no epidemic (see the third graph). This is cause most people with active TB disease, it means the more vulnerable being infected.

IV. NUMERICAL SIMULATIONS
Numerical simulations will be discussed for the spreading of MDR-TB disease by reviewing some of the parameter conditions that are known to be most influential based on the sensitivity index. The following simulation is focused on endemic equilibrium points with ρ = ψ = 0.01 (other parameters can be seen in Table I) and the initial condition is chosen the same in the previous section.
To analyze the effect of parameter β transmission rates of susceptible individuals infected with active TB disease, then in this work we do variations of β that can be seen in the following numerical results.
In Figure 9. when we vary β from 35% to 60%, we see that the impact of the rate of contact is a bit larger, especially for vulnerable individuals.
To see the effect of parameter k level of treatment on patients with active TB on the spread of MDR-TB disease,    then in this work we make variations of k that can be seen in the following numerical results.
When we vary k from 25% to 60%, we see in Figure 10 that efficacy of treatment in active TB individuals is a bit larger, especially for the vulnerable and vaccinated individuals. Despite that the impact is a bit larger, in this case, the number of infected individuals always remain high, especially for the latent TB and patients with pulmonary MDR-TB.
The following is a numerical simulation of the effect of the variation of parameter γ the level of progression of Mtb germ resistance to MDR TB on each population.
When we vary γ from 25% to 60%, then in Figure 11. show that the effect of the development of active TB into MDR-TB is a bit larger, especially for the susceptible.
The following is a simulation of variations in the rate of vaccine intervention (ρ) to individuals susceptible to the spread of MDR-TB with high contact rate (β = 0.6).
When we vary ρ until 30%, we see now in Figure 12 that the impact of vaccination is a bit larger. Despite that the impact is a bit larger, the number of infected individuals always remains high. This implies that vaccination in vulnerable individuals alone is not sufficient to control the disease effectively.
In general, Figure 9-12. gives the trajectory plot of the model (1) for different initial conditions with R 0 > 1. From this figure, we can observe that infected individuals are always present in the population. This means that the trajectories converge to the endemic equilibrium point. Thus, whenever R 0 > 1, the disease persists in the host population  Fig. 12. Evaluation of susceptible S(t), exposed E(t), infectious I(t), individuals with pulmonary MDR-TB R md (t), recovered R(t), and vaccinated V (t) of the MDR-TB model (1) with varying of ρ.
as established in Theorem 2. In Figure 9. We report the result that the higher transmission rate, the higher the number of infected individuals by Mtb. This returns with the number of vulnerable, vaccinated, and cured individuals that are getting lower. Figure 10. interprets that increasing the rate of treatment in active TB patients has a positive effect on all populations, especially the number of infected TB individuals is decreasing. Further, in Figure 11, we can see that the stable progression of patients with pulmonary MDR-TB clearly causes the number of individuals to be infected and recovered to decrease. Moreover, Figure 12. informs that increasing the vaccination rate to vulnerable individuals gives a better impact on the population.

V. CONCLUSIONS
The results of the work show that giving vaccines to vulnerable individuals or knowing that latent individuals make the disease clear from the population. These findings can be used as information for researchers or organizations that produce vaccines for the prevention of TB disease in vulnerable adolescents and also adults as previously stated [2]. While waiting for the vaccine to be available, it is important to think of a strategy to reduce or limit contact between susceptible and individuals with TB disease. The other results show the need for appropriate health interventions in both vulnerable and active TB disease as a bridge causing the outbreak of TB in population. One of the recommendations in the presence of quarantine or isolation (severe) for individuals with TB disease which carried out in hospitals to reduce the likelihood of contact with family members especially those causing latent TB. This is in line with research that has been previously conducted [11]. Moreover, the need for health socialization to the public about awareness of the importance of maintaining a healthy lifestyle. Especially, TB sufferers must understand that TB is not cough, so it is crucial to treatment with early and regularly to prevent getting MDR-TB.
The future work in this model is required for a control class to limit contact between active TB and vulnerable as previously stated by distinguishing two class, i.e. infant and adolescents class. This is done in conjunction with the vaccine currently being developed. Also, further researchers can consider the other biological facts such as the development of active TB is also influenced by exogenous reinfection and a history of other illnesses suffered by patients [12]. One of them is a person with latent TB infection, HIV infection is a risk factor that plays a role in its development into an active TB disease [13].