A Mathematical Model of Tuberculosis with Drug Resistance Effects

Despite the enormous progress in prevention and treatment, tuberculosis disease remains a leading cause of death worldwide and one of the major sources of concern is the drug resistant strain, MDR-TB (multidrug resistant tuberculosis) and XDR-TB (extensively drug resistant tuberculosis). In this work, we extend the standard SEIRS epidemiology model of tuberculosis to include MDR-TB. For that, we considered compartments of susceptible, exposed, infected, resistant to a first line of treatment and recovered humans and we modeled the natural growth, the interactions between these populations and the effects of treatments. We calculate the basic reproduction number, 0  , using the next generation method. The DFE and the EE are established and their stability analysis done to show that they are locally and globally asymptotically stable. Numerical analysis for the model with and without delay is done and demonstrated that in the case of patients with both active tuberculosis and MDR tuberculosis, both strains will still persist due to lack of permanent immunity to tuberculosis while the recovered can still lose their immunity to become susceptible again.


Introduction
Tuberculosis is an airborne disease caused by Mycobacterium tuberculosis (Mtb) bacteria [1]. It is an ancient In addition to posing a major health concern to low and middle income countries, tuberculosis affects economic growth negatively [4]. Psycho-social distress that communities go through is enormous. This involves thinking about the loss of their loved ones and the economic impact of taking care of the sick ones especially among the low income individuals. This impacts not only the individuals, but also the economic progress of the country. Over the last twenty five years, the mortality rate of tuberculosis has greatly decreased by forty five percent since and this is largely due to effective diagnosis and treatment [5]. However, the world is still far from defeating the disease. About 8 billion US dollars per year is needed for a full response to the global tuberculosis epidemic in low and middle income countries by the year 2015 with a funding gap of 2 billion US dollars per year. This amount excluded resources required for research and development, which was estimated to be about 2 billion US dollars yearly [2]. Clearly, this reveals that the current investment in tuberculosis falls below the low and middle income country's needs.
(Mtb) bacteria spreads through inhaling droplets from the cough or sneeze of a person suffering from active tuberculosis. The bacteria enters the body causing an Mtb infection affecting majorly the lungs but it can also affect any other part of the body including the urinary tract, brain, lymph nodes, bones, joints and the ear. Person(s) with lowered immunity such as those with HIV, diabetes, immune disorders, end-stage renal disease, those on drugs that suppress immunity, young children, pregnant women among others are at a higher risk of contracting the disease [6]. Also people suffering from malnutrition due to lifestyle, drug abuse or poverty equally are at risk of contracting tuberculosis. Another category of people largely at risk of contracting tuberculosis are those who work closely or live close to a person with active tuberculosis and they could include health-care workers, people living in crowded living spaces or confined places such as schools or prisons [7].
The intensity of transmission depends on factors related to the bacteria, the human host, the environment and migration. Non-climatic factors such as environmental development and urbanization, population movement and migration affect the severity and incidence of tuberculosis. When it comes to environment and urbanization, the incidence of tuberculosis is generally lower in prime urban areas than in rural areas as there is difficulty in accessing proper medical care in rural villages compared to urban areas. However, rapid urbanization of areas within or on the outskirts of urban centers is commonly done in an uncontrolled fashion without thought or planning. The settlers are mainly migrant workers from rural villages and they tend to settle mostly in poor, overcrowded houses commonly known as slums with hardly any proper sanitation and this in turn leads to increased exposure of the population to Mtb hence a possibility of amplification of the disease to epidemic proportions through lack of effective treatment [6]. Population movements have significant implications for tuberculosis transmission as migration introduces tuberculosis problem to the areas to which the migrants migrate to. Temporary migrant workers often bring the bacteria to lower prevalence areas and local transmission can be readily established [5].
Tuberculosis is curable provided an early diagnosis is made and one follows the proper treatment regimen which could take six months upto two years for the active tuberculosis to clear [8]. There is an emerging form of tuberculosis commonly known as multi-drug resistant (MDR) tuberculosis which is defined as tuberculosis resistant to both of the two most effective first line of antibiotic treatment of active tuberculosis i.e., isoniazid (INH) and rifampin (RIF), and it is harder and more expensive to treat [6]. It is currently a major health concern to medical workers and researchers and one can get MDR tuberculosis by either spending time with an MDR patient and breathing in the MDR tuberculosis bacteria or those with active tuberculosis not following their prescribed treatment regimen or TB medicine not being readily available to them. MDR tuberculosis is much more difficult to treat and the mortality of persons with this form of tuberculosis is far much higher if the second line of antibiotic treatment is not effected promptly [6].
The first tuberculosis model for the transmission dynamics of tuberculosis was developed by Waaler and Anderson (1962) [9] who used a mathematical model to study the epidemiology of tuberculosis. With time other models have been developed to help prevent the risk of transmission of tuberculosis [10]. Recent global reports of multidrug resistant and extensively drug-resistant tuberculosis have renewed concerns that antibiotic resistance may undermine progress in tuberculosis control [11]. Including MDR tuberculosis in mathematical models is relatively new and there are very few models with this aspect. Agusto et al. [12] used a deterministic model with isolation where they studied the transmission dynamics of three strains of mycobacterium tuberculosis (TB), namely, the drug sensitive, multi-drug-resistance (MDR) and extensively-drug-resistance (XDR) tuberculosis strains. Their result of the global sensitivity analysis indicated that the dominant parameters are the disease progression rate, the recovery rate, the infectivity parameter, the isolation rate, the rate of cost to follow up and fraction of fast progression rates. They also found that an increase in isolation rate leads to a decrease in the total number of individuals who are cost to follow up. Trauer, Denholm and McBryde (2014) [8] used a mathematical model to simulate tuberculosis transmission in the highly endemic regions of the Asia-pacific. They found out that their model could not be calibrated to the estimated incidence rate without allowing for re-infection during latency and that even in the presence of a moderate fitness cost and a lower value of 0  , MDR tuberculosis becomes the dominant strain at equilibrium. Improved treatment of drug susceptible tuberculosis did not result in decreased rates of MDR tuberculosis through prevention of the new resistance but rather resulted in a modest increase in a modest increase in MDR tuberculosis. Cohen and Murray (2004) [13] modelled epidemics of multi-drug resistant tuberculosis of heterogeneous fitness. Their model suggested that the threat of multidrug resistance to tuberculosis control depends on relative fitness of MDR strains and this implied that the strains is considerably less than that of drug-sensitive strains and that the emergence of resistance would not jeopardize the success of tuberculosis control efforts. Their results implied that current epidemiological measures and short-term trends in the burden of MDR tuberculosis do not provide evidence that MDR tuberculosis strains can be contained in the absence of specific efforts to limit transmission from those with MDR disease.

The Model Equations
We will first extend the standard SEIRS mathematical model for the transmission of tuberculosis which will demonstrate the transmission of the Mycobacterium tuberculosis in human hosts taking into account the multidrug resistant (MDR) tuberculosis. Most classical models developed for studying tuberculosis dynamics often ignore the resistant class and if they include them, they end up with very complicated models. We seek to present a simple model that can easily be analysed so as to properly understand the dynamics of this disease.
Humans can contract Mtb tuberculosis through contact with individuals who are infected with the disease. After which they enter the exposed (latent) phase where a proportion of this class develop active tuberculosis thus moving into the infectious class. If treatment is administered promptly, those who recover from the disease will move to the recovered class and those who delay treatment and develop MDR tuberculosis will move to the resistant class. Those who recover from MDR tuberculosis will move to the recovered class. Given that there is no permanent immunity to tuberculosis, the recovered can lose their immunity and become susceptible again. Figure 2 represents the flow of individuals into the different compartments and it has been constructed with these assumptions: recruitment is by births only, a variable population, a constant mortality rate, no permanent immunity to tuberculosis, no immediate infectivity. Parameters: λ , recruitment by births, β , force of infection, µ , death rate of humans, α , disease induced deaths in humans due to active tuberculosis, 1 α , disease induced deaths in humans due to MDR tuberculosis, γ , recovery rate of infected humans from active tuberculosis, δ , recovery rate of infected humans from MDR tuberculosis, h λ ,  , rate at which the exposed become infectious, σ , rate at which the infected with active tuberculosis become resistant ρ , rate at which recovered humans become susceptible.
The human population is categorized into five classes such that at time 0 t ≥ there are S, susceptible humans, E, exposed humans to tuberculosis, I, infected humans with active tuberculosis, ES R , resistant humans to the first line of treatment, H R , recovered humans. Thus the size of the human population is given as N = S + E + I + R ES + R. In our model, the recruitment into the susceptible human population is by births ( λ ). The size of the human population is further increased by the partially immune humans in R after they lose their immunity at the rate ρ . The size of the human population is decreased by natural deaths ( µ ) and exposure to Mtb. The exposed susceptibles to Mtb move to the exposed classes E with the force of infection being β resulting in an increase in the exposed class. The exposed class is further decreased by natural deaths ( µ ) and the proportion who move to the infected class I after developing active tuberculosis. The the infected class I is also reduced by natural deaths ( µ ), disease induced deaths ( α ), those who recover ( γ ) and also by those resistant to the first line of treatment ( σ ). Thus both the infected class (I) and the resistant class ( ES R ) gain partial immunity at the rates ( γ ) and ( δ ) respectively thus moving to the recovered class R thus reducing their respective classes and also increasing the recovered class. The resistant class ( ES R ) is also reduced by natural deaths ( µ ) and disease induced deaths ( 1 α ) while the recovered class is reduced by natural deaths ( µ ) and those who lose their partial immunity at the rate ρ . See Table 1 and Table 2 for the description of the state variables and parameters respectively which will be followed by the resulting differential equations. (1)

Equilibrium Points
We obtain the equilibrium points for the system of differential Equation (1) above by equating each of the equations to 0 as shown below Calculation results in two equilibrium points, one being the disease free equilibrium while the other being the endemic equilibrium. With:

Condition of Existence/Positivity
The system will remain positive provided this condition holds: Substituting for x yields: The expression obtained above then is the condition of existence. Lemma Let the initial data be

Calculation of the Reproduction Number
The reproduction number 0  is defined as the average number of secondary cases arising from an average primary case in an entirely susceptible population over the period of infection [14]. The reproduction number is used to predict whether the epidemic will spread or die out. Any epidemiological model has a disease free equilibrium (DFE) at which the population remains in the absence of the disease. According to Diekmann and Heesterbeek (2000) [14], we call 1 FV − the next generation matrix for the model and set the reproduction number, Let us look at the following system of differential equations.
The above system can be represented in matrix form as shown below where  is the Jacobian of the matrix of infection rates and  is the Jacobian of the matrix of transition rates at , 0, 0, 0, 0 N λ µ Next we used a Maple Software to obtain 1 −  and the generated result is given below: We then obtain the spectral radius of

Stability Analysis
Local Stability of the DFE In this section we will determine the stability of the disease free equilibrium points. We will establish this by linearizing the system of differential Equations (1) by obtaining its Jacobian at the disease free equilibrium . The Jacobian of system of differential Equations (1) is as shown: Next, we will establish the eigenvalues for the linearized system as follows: Thus, the disease free equilibrium is locally asymptotically stable given that all the eigenvalues 0 Λ < .

Global Stability of the Disease Free Equilibrium
The local dynamics of a general SEIRS model is determined by the reproduction number 0  . If 0 1 ≤  , then each infected individual in its entire period of infectiousness will produce less than one infected individual on average. This means that the disease will be wiped out of the population. If 0 1 >  , then each infected individual in its entire infectious period having contact with susceptible individuals will produce more than one infected individual implying that the disease persists in the population. If 0 1 =  , and this is defined as the disease threshold, then one individual infects one more individual. For 0 1 ≤  the disease free equilibrium is locally asymptotically stable while for 0 1 >  the disease free equilibrium becomes unstable. By using the theory of Lasalle-Lyapunov function V, we will show the global asymptotic stability. The disease free equili- This can be written as ( ) then the derivative along the trajectories is given by Thus by Lasalle's invariance principle the disease free equilibrium is globally asymptotically stable on Ω .

Global Stability of the Endemic Equilibrium Theorem 2.2
The endemic equilibrium ( ) given by Equation (3) is globally asymptotically stable on Ω . Proof To establish the global stability of the endemic equilibrium, Φ , we construct the lyapunov function 1 : V as described by Ullah, Zaman and Islam (2013) [15] and it is given as ( ) where 1 2 3 , , L L L are positive constants to be later considered. Taking the derivative of the lyapunov function 1 V as given in Equation (6) yields Thus from Equation (9) We define the set ( ) . Therefore the largest compact invariant set is the singleton set Φ which is the endemic equilibrium. By Lasalle Invariance principle Φ is globally asymptotically stable on Ω .

Numerical Simulation
In this section we shall explore the behavior of tuberculosis disease when introduced into a naive environment and conduct numerical simulations of an isolated system (that is, no immigration or emigration). We will pay particular attention to the infected class and the resistance class and the effect of delay in intervention. We used two types of delays: a time dependent and a constant. The model uses a yearly time step and is solved by a fourth order Runge-Kutta scheme. For each simulation, we start with 5 15 10 × susceptible humans, 5 4 10 × exposed humans, 5 1 10 × infected humans, 5 1 10 × resistant humans and 5 1 10 × recovered humans. We will run simulations, in an interval of 100 years, to assess the effect of delaying treatment in the resistant class and also obtaining the numerical solution for the system (1). We start by defining the the values of the parameters in Table 3.

Without delay:
Here, we assume that the treatment has the same effect for both the infected I and the resistant ES R . To treat this case, we had to numerically solve system (1) as shown in Figure 3.
To test the sensitivity of our system to the disease induced death rate 1 α of the resistant individuals ES R and the disease induced death rate α of the infected I, we draw a 3-D curve (Figure 4) perturbing the system using this expression 1 α α η = + where η is a small positive value. Figure 4 shows that when η increases then 1 α the death rate of the resistant increases implying the decrease in R as the result of increasing deaths in I and ES R . Note that for 0.3 η = the ES R class is much greater than I implying much more deaths.
With delay: Now, we suppose that the treatment affects each of the individuals I and ES R in a different way. For that, we  add a delay function to the last equations of system (1) as shown below To better investigate on the drug efficiency we used a function ( ) Comparing to Figure 3, the I and ES R curves coincide at the beginning and then separate and this clearly shows that eventually the drugs ends up having a different effect on each of the populations. This effect is seen on the recovered class too.
After that, we used a low delay with a short period 0.01 ω = as shown in Figure 6 and a fast one with a longer period 0.1 ω = as shown in Figure 7.   What we see in Figure 6 is the positive effect of the drugs on the resistant class as there is reduction in the resistant class and an increase in the recovered class at exactly the same time whereas Figure 7 shows that a fast time dependent gives similar results as Figure 3. This implies that the drugs have a similar effect on both I and ES R .

Conclusion
This study presents a simple yet more realistic deterministic model for the transmission dynamics of tuberculosis. In contrast to many tuberculosis models in literature, we include the class resistant to the first line of treatment for tuberculosis. The resistant category is of particular importance in modeling tuberculosis in low and middle income countries mostly found in Africa where public health is under developed. In particular, the number of individuals who seek medical intervention and those who actually recover from both active tuberculosis and MDR tuberculosis is worth noting. MDR parameter is used to measure the success of treatment administered to MDR patients. The recovery parameters are also used to measure the implications of failure in treating both active tuberculosis and MDR tuberculosis. Our results and simulations demonstrate that when MDR tuberculosis patients delay treatment or fail to go for treatment altogether this strain will still persist over time. Also in the case of patients with both active tuberculosis and MDR tuberculosis, we establish that both strains will still persist at some equilibrium since there is no permanent immunity to tuberculosis and the recovered can still lose their immunity and become susceptible again. Of particular concern is the relation of tuberculosis and HIV as most tuberculosis patients are aware of their HIV status. We hope a lot of efforts, resources and research will be placed for a rigorous analysis of this co-infection model so as to fully understand this problem and the interventions for possible prevention of this predicament.