A Mathematical Model of Tuberculosis with Drug Resistance Effects ()
Received 11 February 2016; accepted 22 July 2016; published 25 July 2016
1. Introduction
Tuberculosis is an airborne disease caused by Mycobacterium tuberculosis (Mtb) bacteria [1] . It is an ancient disease with evidence of its existence being found in relics from ancient Egypt, India and China [1] . In the eighteenth century, Western Europe suffered terribly from this disease with a prevalence as high as 900 deaths per 100,000. This was largely due to poor ventilation, overcrowded housing, primitive sanitation, malnutrition among other risk factors that led to the epidemic [1] . Today, this disease ranks as the second leading cause of morbidity and mortality in the world from a single infectious agent, after the human immunodeficiency virus (HIV) [2] . Interestingly, about one-third of the world’s population is infected with Mtb with approximately nine million people developing active tuberculosis and up to nearly two million people worldwide die from the disease every year. In 2013, approximately nine million people contracted active tuberculosis and this included 1.1 million cases among people living with HIV and 550,000 children. Out of these nine million cases 1.5 million people succumbed to the disease and this included 360,000 among people who were HIV-positive, 510,000 were women out of which 180,000 were HIV-positive. Africa recorded the highest tuberculosis/HIV burden with three out of four Tuberculosis patients knowing their HIV status. Approximately 480,000 people developed multidrug-resistant (MDR) tuberculosis globally with 210,000 of those who developed MDR tuberculosis succumbing to it [2] . Figure 1 showed the global incidence of tuberculosis in 2008 [3] .
![]()
Figure 1. Global distribution of tuberculosis [3] .
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 tubercu- losis 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
, 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.
2. Model Analysis
2.1. 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 multi- drug 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 pre- sent 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.
![]()
Figure 2. Schematics of the compartmental model. State variables: S, susceptible human, E, exposed human with Mtb, I, Infected human with active tuberculosis, 
, resistant humans to first line of treatment, R, recovered humans from both active tuberculosis and MDR tuberculosis. Parameters:
, recruitment by births, 
, force of infection, 
, death rate of humans, 
, disease induced deaths in humans due to active tuberculosis, 
, 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, 
, 
, 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 ![]()
there are S, susceptible humans, E, exposed humans to tuberculosis, I, infected humans with active tuberculosis, ![]()
, resistant humans to the first line of treatment, ![]()
, recovered humans. Thus the size of the human population is given as N = S + E + I + RES + 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 (![]()
) 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 (![]()
) is also reduced by natural deaths (![]()
) and disease induced deaths (![]()
) 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.
![]()
Table 1. Description of state variables.
![]()
Table 2. Description of parameters.
2.2. Differential Equations
![]()
(1)
2.3. 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
![]()
(2)
Calculation results in two equilibrium points, one being the disease free equilibrium while the other being the endemic equilibrium.
Disease Free Equilibrium Point ![]()
is expressed as follows:
![]()
Endemic Equilibrium Point ![]()
is expressed as follows:
![]()
(3)
With: ![]()
2.4. 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![]()
. Then the solution set
![]()
of the system (1) is positive for all![]()
.
2.5. Calculation of the Reproduction Number
The reproduction number ![]()
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 ![]()
the next generation matrix for the model and set the reproduction num-
ber, ![]()
where ![]()
and ![]()
for ![]()
for the number of compartments,
and ![]()
for the infected compartments only. ![]()
denotes the spectral radius of a matrix A. F and V are ![]()
matrices, where m is the number of infected classes. Consider an infected individual introduced into compartment k of a disease-free population. The ![]()
entry of F is the rate at which an infected indi- vidual in compartment j produces new infections in compartment i, and the ![]()
entry of ![]()
is the average time an infected individual spends in compartment j during its lifetime, assuming that the population remains near the DFE and barring reinfection Hence, the ![]()
entry of the product ![]()
is the expected number of new infections in compartment i produced by the infected individual originally introduced into compartment k. Let us look at the following system of differential equations.
![]()
(4)
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 ![]()
so that
![]()
Next we used a Maple Software to obtain ![]()
and the generated result is given below:
![]()
We then obtain the spectral radius of![]()
, ![]()
, which is defined as the largest eigenvalue of ![]()
and the spectral radius for the above system is the basic reproduction number, ![]()
, and its expression is given in equation (5) below
![]()
(5)
2.6. 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:
![]()
Evaluating the determinant above yields these eigenvalues:
![]()
, ![]()
, ![]()
, ![]()
, ![]()
Thus, the disease free equilibrium is locally asymptotically stable given that all the eigenvalues![]()
.
Global Stability of the Disease Free Equilibrium
The local dynamics of a general SEIRS model is determined by the reproduction number![]()
. If![]()
, 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![]()
, 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![]()
, and this is defined as the disease threshold, then one individual infects one more individual. For ![]()
the disease free equilibrium is locally asymptotically stable while for ![]()
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- brium point is![]()
.
Theorem 2.1
If![]()
, then the disease-free equilibrium ![]()
of the system is globally asymptotically stable on![]()
.
Proof
We construct the following Lasalle-Lyapunov function ![]()
on the positively invariant compact set![]()
. Thus on![]()
, ![]()
is continuous and non negative. We define![]()
.
The system of ordinary differential equations given by Equation (4) can be written as
![]()
This can be written as ![]()
where ![]()
and![]()
.
If we define ![]()
then the derivative along the trajectories is given by ![]()
as
![]()
which is strictly decreasing when ![]()
thus![]()
.
We define the set![]()
. The largest invariant set is contained in the set E
for which ![]()
or ![]()
or![]()
. Thus ![]()
when![]()
. If ![]()
or ![]()
then![]()
. 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 ![]()
![]()
where ![]()
as described by Ullah, Zaman and Islam (2013) [15] and it is given as
![]()
(6)
where ![]()
are positive constants to be later considered.
Taking the derivative of the lyapunov function ![]()
as given in Equation (6) yields
![]()
(7)
![]()
(8)
where ![]()
and![]()
.
Choosing![]()
, Equation (8) becomes
![]()
(9)
Thus from Equation (9), ![]()
if and only if ![]()
and ![]()
We have that ![]()
if and
only if![]()
, ![]()
and ![]()
or when ![]()
and![]()
. 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![]()
.
3. 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 ![]()
susceptible humans, ![]()
exposed humans, ![]()
infected humans, ![]()
resistant humans and ![]()
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![]()
. 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 ![]()
of the resistant individuals ![]()
and the disease induced death rate ![]()
of the infected I, we draw a 3-D curve (Figure 4) perturbing the system using this expression ![]()
where ![]()
is a small positive value.
Figure 4 shows that when ![]()
increases then ![]()
the death rate of the resistant increases implying the decrease in R as the result of increasing deaths in I and![]()
. Note that for ![]()
the ![]()
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 ![]()
in a different way. For that, we
![]()
Figure 3. Evolution of each population without delay.
add a delay function to the last equations of system (1) as shown below
![]()
(10)
To better investigate on the drug efficiency we used a function ![]()
for two different time delays.
![]()
Figure 5 illustrates the dynamics of the populations under the effect of constant 5 years treatment delay.
Comparing to Figure 3, the I and ![]()
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 ![]()
as shown in Figure 6 and a fast one with a longer period ![]()
as shown in Figure 7.
![]()
Figure 5. Evolution of each population with constant delay.
![]()
Figure 6. Evolution of each population with the low time dependent delay.
![]()
Figure 7. Evolution of each population with the fast time dependent delay.
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![]()
.
4. 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 inter-
ventions for possible prevention of this predicament.
Acknowledgements
We would like to start by thanking all the organizers and sponsors of the CIMPA Kenya School 2015 making special mention of Prof. Sylvain Sylva in 2011, CIMPA coordinator for Sub Saharan Africa, for granting us an opportunity to participate in this school and for financially supporting our participation which gave us the opportunity to meet great mathematicians from all parts of the world. Our special thanks go to Prof. Dorothy Musekwa, Dr. Josephine Kagunda and Dr. Nelson Owuor who inspired us through an exercise that resulted in this work. We would wish to thank Prof. Wilson Lamb, Prof. Jacek Banasiak, Prof. Luboobi, Prof. Wandera Ogana, Dr. Winnie Mutuku, Dr. Irene Wattanga, Dr. Samuel Mwalili and Prof. David Malonza who offered us great insights as we formulated this problem. Lastly, we would like to thank our fellow participants from all over the world who offered us encouragement and sincere advice to better this work. May God bless each one of you richly.