Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Stability of a Two-Strain Tuberculosis Model with General Contact Rate

Stability of a Two-Strain Tuberculosis Model with General Contact Rate Stability of a Two-Strain Tuberculosis Model with General Contact Rate //// Hindawi Publishing Corporation Home Journals About Us About this Journal Submit a Manuscript Table of Contents Journal Menu Abstracting and Indexing Aims and Scope Annual Issues Article Processing Charges Articles in Press Author Guidelines Bibliographic Information Contact Information Editorial Board Editorial Workflow Free eTOC Alerts Reviewers Acknowledgment Subscription Information Open Special Issues Published Special Issues Special Issue Guidelines Abstract Full-Text PDF Full-Text HTML Linked References How to Cite this Article Abstract and Applied Analysis Volume 2010 (2010), Article ID 293747, 31 pages doi:10.1155/2010/293747 Research Article <h2>Stability of a Two-Strain Tuberculosis Model with General Contact Rate</h2> Hai-Feng Huo , 1 Shuai-Jun Dang , 1 and Yu-Ning Li 2 1 Institute of Applied Mathematics, Lanzhou University of Technology, Lanzhou, Gansu 730050, China 2 Department of Pediatrics, First Hospital of Lanzhou University, Lanzhou, Gansu 730000, China Received 1 November 2010; Accepted 31 December 2010 Academic Editor: D. Anderson Copyright © 2010 Hai-Feng Huo et al. This is an open access article distributed under the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract A two-strain tuberculosis model with general contact rate which allows tuberculosis patients with the drug-sensitive Mycobacterium tuberculosis strain to be treated is presented. The model includes both drug-sensitive and drug-resistant strains. A detailed qualitative analysis about positivity, boundedness, existence, uniqueness and global stability of the equilibria of this model is carried out. Analytical results of the model show that the quantities 𝑅 1 and 𝑅 2 , which represent the basic reproduction numbers of the sensitive and resistant strains, respectively, provide the threshold conditions which determine the competitive outcomes of the two strains. Numerical simulations are also conducted to confirm and extend the analytic results. 1. Introduction Tuberculosis (TB), an infectious disease caused by Mycobacterium tuberculosis ( M. tuberculosis ), is one of the world's leading causes of infectious mortality. According to the World Health Organization(WHO), one third of the world's population is infected with M. tuberculosis , leading to between two and three millions death each year. At present, about 95% of the estimated 8 million new cases of TB occurring each year are in developing countries, where 80% occur among people between the ages of 15 and 59 years [ 1 ]. For the time being, TB is becoming a world-wide problem. War, famine, homelessness, and a lack of medical care all contribute to the increasing incidence of tuberculosis among disadvantaged persons. Since TB is easily transmissible between persons, the increase in TB in any segment of the population represents a threat to all segments of the population. Sub-Saharan Africa remains the epicenter of the epidemic, but India, China, Indonesia, Bangladesh and Pakistan together account for more than half of the cases in the world [ 2 ]. TB was assumed to be on its way out in developed countries until the number of TB cases began to increase in the 1980s. With this return, we face the paradox of a well-known bacteria, fully treatable with efficient and affordable drugs according to internationally recommended guidelines, which yet causes increasing human suffering and death. Active TB cases may be pulmonary or extrapulmonary, but pulmonary cases are more infectious and form the bulk of most cases of active TB. The usual symptoms of active TB include fatigue, high fever, night sweats, loss of appetite, and a cough, but confirmation of active TB requires a positive sputum culture. Extrapulmonary accounts for between 5% and 30% of the total cases and may affect any part of the body. M. tuberculosis droplets are released into the air by sneezing and coughing of infectious individuals. Tubercle bacillus carried by such droplets live in the air for a short period of time [ 3 ]; it is believed that occasional contacts with an infectious individual rarely lead to transmission. TB is described as a slow disease because of its long and variable latency period distribution and its short and relatively narrow infectious period distribution. Individuals who are latently infected are neither clinically ill nor capable of transmitting TB [ 4 ]. Most latently infected individuals do not become infectious (active TB). About 5%–10% of the latently infected individuals develop active TB, that is, about 90%–95% remain latently infected. Antituberculosis drugs are a two-edged sword. While they destroy pathogenic M. tuberculosis , they also select for drug-resistant bacteria against with those drugs which are then ineffective. Global surveillance has shown that drug-resistant TB is widespread and is now a threat to TB control programs in many countries. The WHO distinguishes between two types of resistance: acquired resistance—resistance among previously treated patients; and primary resistance—resistance among new cases. In all regions studied, prevalence of acquired resistance is higher than prevalence of primary resistance, but the size of this difference varies between regions [ 5 ]. Treatment of TB consists of a combination of different drugs to avoid acquisition of resistance. Despite these precautions, drug resistance continues to emerge being favoured by the long duration of treatment and improper use of the antibiotics. Drug resistant TB has higher rates of treatment failure and longer periods of infectiousness in part due to the time lapse between TB diagnosis and obtaining drug-sensitivity test results [ 6 ]. In general, tuberculosis can be treated with antibiotics. However, most worrisome is resistance to the two first-line drugs, isoniazid and rifampicin, defined as multidrug resistance(MDR)-TB, which is an emerging problem of great importance to public health, with higher mortality rates than drug-sensitive TB, particularly in immunocompromised patients. MDR-TB patients require treatment with more toxic second line drugs and remain infectious for longer than patients infected with drug-sensitive strains, incurring higher costs due to prolonged hospitalization. Mathematical models can provide useful tools to analyze the spread and control of infectious diseases [ 7 ]. Different mathematical models for TB have been formulated and studied [ 3 , 5 , 8 – 16 ]. In most of the epidemics models for TB discussed in the literature, the question of contact rate has not been a central one. Nevertheless, the mode of transmission is crucially important for two reasons. First, it determines the probable response of the disease to control, second, the objective in many models of disease in animals is to predict what will happen when a pathogen is introduced into a system in which it does not currently exist [ 17 ]. Even more important, several laboratory studies have found that the mass action incidence rate is inadequate for describing pathogen transmission and standard incidence rate is considered, for human disease, more accurate than mass action incidence rate [ 17 ]. Moreover, epidemic models with nonlinear incidence rate have recently attracted a great deal of attention from mathematical models [ 18 – 22 ]. In fact, the incidence of a disease which is the number of new cases per unit time plays an important role in the study of mathematical epidemiology. Thieme and Castillo-Chavez [ 23 ] argued that the general form of a population-size-dependent incidence should be written as 𝛽 𝐶 ( 𝑁 ) 𝑆 𝐼 / 𝑁 , where 𝑆 and 𝐼 are the numbers of susceptible and infective at time 𝑡 , respectively, 𝛽 is the probability of transmitting the infection between two individuals taking part in a contact per unit time, 𝐶 ( 𝑁 ) is the probability for an individual to take part in a contact and 𝑁 is the size of the total population. In [ 24 ], Zhang and Ma think the above incidence frequently takes two forms in most of the literature. When 𝐶 ( 𝑁 ) = 𝑁 , the corresponding incidence is the bilinear form 𝛽 𝑁 𝑆 𝐼 / 𝑁 = 𝛽 𝑆 𝐼 . When 𝐶 ( 𝑁 ) = 𝜆 , the corresponding incidence 𝜆 𝑆 𝐼 / 𝑁 is called standard form. When the total population size 𝑁 is not quite large, since the number of contacts made by an individual per unit time should increase as the total population size 𝑁 increases, the bilinear form would be suitable, but when the total population size is too large, since the number of contacts made by an infective per unit time should be limited, or should grow less rapidly as the total population size 𝑁 increases, the bilinear form is not suitable and the standard form may be more realistic [ 24 ]. Therefore, the two forms of incidence mentioned above are actually two extreme cases for the total population size 𝑁 being very small and very large, respectively. Moreover, Heesterbeek and Metz [ 25 ] derived the expression for the saturating contact rate of individuals contacts in a population that mixes randomly, that is, 𝐶 ( 𝑁 ) = 𝑏 𝑁 √ 1 + 𝑏 𝑁 + , 1 + 2 𝑏 𝑁 ( 1 . 1 ) where 𝐶 ( 𝑁 ) is nondecreasing and 𝐶 ( 𝑁 ) / 𝑁 is nonincreasing. In Bowong and Tewa [ 26 ], a SEI type of tuberculosis model with a general contact rate is proposed to study the global asymptotically stability of disease-free equilibrium and endemic equilibrium. The stability of equilibria is derived through the use of Lyapunov stability theory and LaSalle's invariant set theorem. Nevertheless, the effect of drug-resistant strain of tuberculosis is not taken into account in their article. In this paper, we study a tuberculosis model with general contact rate that includes explicitly both drug-sensitive strain and drug-resistant strain. We adapt the approach of Bowong and Tewa [ 26 ] for modeling the effective chemoprophylaxis (given to latently infected individuals) and therapeutic treatments (given to infectious). We extend their model by including the latently infected with resistant strain class and infectious with resistant strain class. The new model allows us to examine the effects of drug treatment on the prevalence of both drug-sensitive and drug-resistant strains. Mathematical properties of the model system are studied both analytically and numerically. It is shown that the system has three possible equilibrium points including an endemic equilibrium at which both strains are present. A detailed analysis of global asymptotically stability is conducted, which shows that the dynamic behaviors of the system are determined by two quantities, 𝑅 1 and 𝑅 2 , which represent the basic reproduction numbers of the sensitive and resistant strains. Without regard to disease-induced death rate, we prove that the disease-free equilibrium is globally asymptotically stable when 𝑅 1 < 1 , 𝑅 2 < 1 , the unique drug-resistant-TB-strain-only equilibrium exists and is globally asymptotically stable when 𝑅 1 < 1 < 𝑅 2 . As for the unique interior equilibrium, existence, uniqueness and global stability can be proved under certain conditions. The paper is organized as follows. In Section 2 , we formulate a two-strain tuberculosis model with general contact rate. The threshold conditions for the existence and uniqueness equilibria are derived and the global asymptotically stability of equilibria are proved in Sections 3 and 4 is devoted to numerical simulations. In Section 5 , we summarize the findings and conclusions. 2. Model Description According to the transmitted features of TB, we subdivide the population into susceptible individuals ( 𝑆 ) , those exposed to drug-sensitive TB ( 𝐸 1 ) , individuals with symptoms of TB and drug-sensitive ( 𝐼 1 ) , those exposed to drug-resistant TB ( 𝐸 2 ) and those displaying symptoms of TB and drug-resistant ( 𝐼 2 ) , the total number of population at time 𝑡 is given by 𝑁 ( 𝑡 ) = 𝑆 ( 𝑡 ) + 𝐸 1 ( 𝑡 ) + 𝐼 1 ( 𝑡 ) + 𝐸 2 ( 𝑡 ) + 𝐼 2 ( 𝑡 ) . ( 2 . 1 ) The model is represented by the transfer diagram in Figure 1 . Here, Λ is the recruitment rate, 𝜇 is nondisease related death rate. Since exposed individuals are not capable of transmitting the disease, we assume that susceptible may become infected only through contacts with active infectious individual at a rate 𝛽 𝑖 ( 𝑁 ) 𝐼 𝑖 , 𝑖 = 1 , 2 , where 𝑖 = 1 and 𝑖 = 2 represent rates of infection by drug-sensitive strain and drug-resistant strain, and the transmission coefficient 𝛽 𝑖 ( 𝑁 ) is a nonnegative 𝐶 2 function of the total population 𝑁 . Susceptible are infected with drug-sensitive strain and drug-resistant strain entering classes 𝐸 1 and 𝐸 2 , respectively. For simplicity, we only account for treatment of sensitive strain, we assume that chemoprophylaxis of individuals in 𝐸 1 reduces their reactivation at a constant rate 𝑟 1 and that the initiation of therapeutics immediately removes individuals in 𝐼 1 from active status and places them into exposed class at a constant rate 𝑟 2 . The time before individuals in 𝐸 1 who does not received effective chemoprophylaxis become infectious is assumed to satisfy an exponential distribution, with mean waiting time 1 / 𝑘 1 . Thus, individuals leave the class 𝐸 1 to the class 𝐼 1 at a constant rate 𝑘 1 ( 1 − 𝑟 1 ) . The time before individuals in 𝐸 2 become infectious is also assumed to satisfy an exponential distribution, with mean waiting time 1 / 𝑘 2 . Thus, individuals leave the class 𝐸 2 to the class 𝐼 2 at a constant rate 𝑘 2 . Individuals sick with drug-sensitive strain receive treatment at rate 𝑟 2 and a proportion 𝑝 𝑟 2 respond to treatment and move into exposed class 𝐸 1 , in the remaining proportion ( 1 − 𝑝 ) 𝑟 2 , inappropriate treatment results in the development of drug-resistant strain and the individuals move into class 𝐸 2 . Individuals in the infectious class 𝐼 1 and 𝐼 2 suffer additional disease-induced death at rate 𝑑 1 and 𝑑 2 , where 𝑑 2 > 𝑑 1 . Then the following two strains tuberculosis model is formulated 𝑆  = Λ − 𝛽 1 ( 𝑁 ) 𝑆 𝐼 1 − 𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 𝐸 − 𝜇 𝑆 ,  1 = 𝛽 1 ( 𝑁 ) 𝑆 𝐼 1 + 𝑝 𝑟 2 𝐼 1 − 𝑘 1  1 − 𝑟 1  𝐸 1 − 𝜇 𝐸 1 , 𝐼  1 = 𝑘 1  1 − 𝑟 1  𝐸 1 −  𝑟 2 + 𝜇 + 𝑑 1  𝐼 1 , 𝐸  2 = 𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 1 −  𝜇 + 𝑘 2  𝐸 2 , 𝐼  2 = 𝑘 2 𝐸 2 −  𝜇 + 𝑑 2  𝐼 2 , ( 2 . 2 ) where Λ , 𝑘 1 , 𝑘 2 , 𝑟 1 , and 𝑟 2 are assumed to be positive and 𝑝 ∈ [ 0 , 1 ] . It is natural to assume that the transmission coefficient 𝛽 1 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) satisfies the following conditions: 𝛽 1 ( 𝑁 ) > 0 , 𝛽  1  ( 𝑁 ) ≤ 0 , 𝑁 𝛽 1  ( 𝑁 )  𝛽 ≥ 0 , 2 ( 𝑁 ) > 0 , 𝛽  2  ( 𝑁 ) ≤ 0 , 𝑁 𝛽 2  ( 𝑁 )  ≥ 0 . ( 2 . 3 ) Figure 1: Transfer diagram for the dynamics of tuberculosis with general contact rate. Obviously, we can see that 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 corresponds to the standard incidence rate, that 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 corresponds to the mass action incidence rate, and that 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) corresponds to the saturating contact rate, where 𝐶 ( 𝑁 ) = 𝑏 𝑁 √ 1 + 𝑏 𝑁 + . 1 + 2 𝑏 𝑁 ( 2 . 4 ) 2.1. Positivity and Boundedness System ( 2.2 ) describes human population and therefore it is necessary to prove that all the variables 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 are nonnegative for all time. Solutions of system ( 2.2 ) with positive initial data remains positive for all time 𝑡 ≥ 0 and are bounded. Theorem 2.1. If 𝑆 ( 0 ) ≥ 0 , 𝐸 1 ( 0 ) ≥ 0 , 𝐼 1 ( 0 ) ≥ 0 , 𝐸 2 ( 0 ) ≥ 0 , 𝐼 2 ( 0 ) ≥ 0 , the solutions 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2 ( 𝑡 ) of system ( 2.2 ) are positive for 𝑡 ≥ 0 . For system ( 2.2 ), the region Γ is positively invariant and all solutions starting in Γ approach, enter, or stay in Γ , where   Γ = 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2  ( 𝑡 ) ∈ ℝ 5 + Λ ∶ 𝑁 ( 𝑡 ) ≤ 𝜇  . ( 2 . 5 ) Proof. Under the given initial conditions, it is easy to prove that the solutions of the system ( 2.2 ) are positive; if not, we assume a contradiction: that there exists a first time 𝑡 1 such that 𝑆  𝑡 1  = 0 , 𝑆   𝑡 1  < 0 , 𝐸 1 ( 𝑡 ) ≥ 0 , 𝐼 1 ( 𝑡 ) ≥ 0 , 𝐸 2 ( 𝑡 ) ≥ 0 , 𝐼 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 1 , ( 2 . 6 ) there exists a 𝑡 2 𝐸 1  𝑡 2  = 0 , 𝐸  1  𝑡 2  < 0 , 𝑆 ( 𝑡 ) ≥ 0 , 𝐼 1 ( 𝑡 ) ≥ 0 , 𝐸 2 ( 𝑡 ) ≥ 0 , 𝐼 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 2 , ( 2 . 7 ) there exists a 𝑡 3 𝐼 1  𝑡 3  = 0 , 𝐼  1  𝑡 3  < 0 , 𝑆 ( 𝑡 ) ≥ 0 , 𝐸 1 ( 𝑡 ) ≥ 0 , 𝐸 2 ( 𝑡 ) ≥ 0 , 𝐼 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 3 , ( 2 . 8 ) there exists a 𝑡 4 𝐸 2  𝑡 4  = 0 , 𝐸  2  𝑡 4  < 0 , 𝑆 ( 𝑡 ) ≥ 0 , 𝐸 1 ( 𝑡 ) ≥ 0 , 𝐼 1 ( 𝑡 ) ≥ 0 , 𝐼 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 4 , ( 2 . 9 ) or there exists a 𝑡 5 𝐼 2  𝑡 5  = 0 , 𝐼  2  𝑡 5  < 0 , 𝑆 ( 𝑡 ) ≥ 0 , 𝐸 1 ( 𝑡 ) ≥ 0 , 𝐼 1 ( 𝑡 ) ≥ 0 , 𝐸 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 5 . ( 2 . 1 0 ) In the first case, we have 𝑆   𝑡 1  = Λ > 0 , ( 2 . 1 1 ) which is a contradiction meaning that 𝑆 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . In the second case, we have 𝐸  1  𝑡 2  = 𝐼 1  𝑡 2 𝛽   1  𝑡 ( 𝑁 ) 𝑆 2  + 𝑝 𝑟 2  ≥ 0 , ( 2 . 1 2 ) which is a contradiction meaning that 𝐸 1 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . In the third case, we have 𝐼  1  𝑡 3  = 𝑘 1  1 − 𝑟 1  𝐸 1  𝑡 3  ≥ 0 , ( 2 . 1 3 ) which is a contradiction meaning that 𝐼 1 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . In the fourth case, we have 𝐸  2  𝑡 4  = 𝛽 2  𝑡 ( 𝑁 ) 𝑆 4  𝐼 2  𝑡 4  + ( 1 − 𝑝 ) 𝑟 2 𝐼 1  𝑡 4  ≥ 0 , ( 2 . 1 4 ) which is a contradiction meaning that 𝐸 2 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . In the fifth case, we have 𝐼  2  𝑡 5  = 𝑘 2 𝐸 2  𝑡 5  ≥ 0 , ( 2 . 1 5 ) which is a contradiction meaning that 𝐼 2 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . Thus, in all cases, 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2 ( 𝑡 ) remain positive for 𝑡 ≥ 0 . Let ( 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2 ( 𝑡 ) ) ∈ ℝ 5 + be any solution with nonnegative initial condition, adding all equations in system ( 2.2 ) gives 𝑁  ( 𝑡 ) = Λ − 𝜇 𝑁 − 𝑑 1 𝐼 1 − 𝑑 2 𝐼 2 , ≤ Λ − 𝜇 𝑁 . ( 2 . 1 6 ) It follows that Λ 0 ≤ 𝑁 ( 𝑡 ) ≤ 𝜇 + 𝑁 ( 0 ) 𝑒 − 𝜇 𝑡 , ( 2 . 1 7 ) where 𝑁 ( 0 ) represents initial values of the total population. Thus 0 ≤ 𝑁 ( 𝑡 ) ≤ Λ / 𝜇 , as 𝑡 → ∞ . Therefore all feasible solutions of system ( 2.2 ) enter the region   Γ = 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2  ( 𝑡 ) ∈ ℝ 5 + Λ ∶ 𝑁 ( 𝑡 ) ≤ 𝜇  . ( 2 . 1 8 ) Hence, Γ is positively invariant and it is sufficient to consider solutions of system ( 2.2 ) in Γ . Existence, uniqueness and continuation results for system ( 2.2 ) hold in this region. It can be shown that 𝑁 ( 𝑡 ) is bounded and all the solutions starting in Γ approach, enter or stay in Γ . 3. Model Analysis There are one disease-free equilibrium 𝑃 0 and two possible endemic equilibria for system ( 2.2 ), the endemic equilibria include boundary equilibrium 𝑃 1 (when only the drug-resistant TB-strain is present) and the interior equilibrium 𝑃 2 (when both strains exist). 3.1. Global Stability of the Disease-Free Equilibrium The disease-free-equilibrium is given as 𝑃 0 =  𝑆 0  , 0 , 0 , 0 , 0 , 𝑆 0 = Λ 𝜇 . ( 3 . 1 ) The basic reproduction number is defined as the number of secondary infections produced by a single infectious individual during the entire infectious period. In our case the reproduction number defines the number of secondary TB infections produced by a single active TB individual during the entire infectious period. Mathematically it is defined as the spectral radius of the next generation matrix [ 27 ]. To determine the stability of the disease-free steady state 𝑃 0 , we use the results in van den Driessche and Watmough [ 27 ]. Reorder the components of 𝑃 0 as 𝐸 1 = 0 , 𝐼 1 = 0 , 𝐸 2 = 0 , 𝐼 2 = 0 , 𝑆 = Λ / 𝜇 . Set ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐹 ℱ = 1 𝐹 2 𝐹 3 𝐹 4 𝐹 5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝛽 1 ( 𝑁 ) 𝑆 𝐼 1 0 𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑘 𝒱 = 1  1 − 𝑟 1  𝐸 1 + 𝜇 𝐸 1 − 𝑝 𝑟 2 𝐼 1 𝑟 2 𝐼 1 +  𝜇 + 𝑑 1  𝐼 1 − 𝑘 1  1 − 𝑟 1  𝐸 1  𝜇 + 𝑘 2  𝐸 2 − ( 1 − 𝑝 ) 𝑟 2 𝐼 1  𝜇 + 𝑑 2  𝐼 2 − 𝑘 2 𝐸 2 𝜇 𝑆 + 𝛽 1 ( 𝑁 ) 𝑆 𝐼 1 + 𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . − Λ ( 3 . 2 ) The infected compartments are 𝐸 1 , 𝐼 1 , 𝐸 2 , and 𝐼 2 . Thus ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐹 = 0 𝛽 1  𝑆 0  𝑆 0 0 0 0 0 0 0 0 0 0 𝛽 1  𝑆 0  𝑆 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 𝑉 = 𝜇 + 𝑘 1  1 − 𝑟 1  − 𝑝 𝑟 2 0 0 − 𝑘 1  1 − 𝑟 1  𝜇 + 𝑑 1 + 𝑟 2 0 0 0 − ( 1 − 𝑝 ) 𝑟 2 𝜇 + 𝑘 2 0 0 0 − 𝑘 2 𝜇 + 𝑑 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 3 . 3 ) The dominant eigenvalues of 𝐹 𝑉 − 1 are given by 𝑅 1 = 𝑘 1  1 − 𝑟 1  𝛽 1  𝑆 0  𝑆 0  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  , 𝑅 2 = 𝑘 2 𝛽 2  𝑆 0  𝑆 0  𝜇 + 𝑘 2   𝜇 + 𝑑 2  , ( 3 . 4 ) where 𝑅 1 and 𝑅 2 are reproduction numbers for drug-sensitive TB strain only and drug-resistant TB strain only, respectively. It implies that the spectral radius of the matrix 𝐹 𝑉 − 1 is 𝜌  𝐹 𝑉 − 1   𝑅 = m a x 1 , 𝑅 2  . ( 3 . 5 ) If 𝑅 1 < 1 , and 𝑅 2 < 1 , then 𝜌 ( 𝐹 𝑉 − 1 ) < 1 . By Theorem 3.2 in van den Driessche and Watmough [ 27 ], we know that the disease-free steady state 𝑃 0 is locally asymptotically stable, 𝑃 0 is unstable if 𝑅 1 > 1 or 𝑅 2 > 1 . The following theorems provide the global stability of the disease-free equilibrium. To conduct an analytical analysis of global asymptotical behaviors of the disease-free equilibrium point, first of all, we assume that there is no disease-induced death rate, that is, 𝑑 1 = 0 , 𝑑 2 = 0 . Consequently, the total population size 𝑁 ( 𝑡 ) satisfies the equation 𝑑 𝑁 / 𝑑 𝑡 = Λ − 𝜇 𝑁 and 𝑁 ( 𝑡 ) → Λ / 𝜇 as 𝑡 → ∞ . Using results from Castillo-Chavez and Thieme [ 28 ] and Mischaikow et al. [ 29 ], we can obtain analytical results by considering the following limiting system of ( 2.2 ) in which the total population is assumed to be constant 𝑁 = 𝑆 0 = Λ / 𝜇 : 𝐸  1 = 𝛽 1  𝑆 0  𝐼 1  𝑆 0 − 𝐸 1 − 𝐼 1 − 𝐸 2 − 𝐼 2  + 𝑝 𝑟 2 𝐼 1 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝐸   1 , 𝐼  1 = 𝑘 1  1 − 𝑟 1  𝐸 1 −  𝑟 2  𝐼 + 𝜇 1 , 𝐸  2 = 𝛽 2  𝑆 0  𝐼 2  𝑆 0 − 𝐸 1 − 𝐼 1 − 𝐸 2 − 𝐼 2  + ( 1 − 𝑝 ) 𝑟 2 𝐼 1 −  𝜇 + 𝑘 2  𝐸 2 , 𝐼  2 = 𝑘 2 𝐸 2 − 𝜇 𝐼 2 . ( 3 . 6 ) Notice that the 𝑆 equation is eliminated from ( 3.6 ) and the variable 𝑆 is replaced by 𝑆 0 − 𝐸 1 − 𝐼 1 − 𝐸 2 − 𝐼 2 . By introducing the fractions 𝑥 1 = 𝐸 1 𝑁 , 𝑥 2 = 𝐼 1 𝑁 , 𝑦 1 = 𝐸 2 𝑁 , 𝑦 2 = 𝐼 2 𝑁 , ( 3 . 7 ) we obtain the equivalent limiting system from ( 3.6 ) 𝑥  1 = 𝑆 0 𝛽 1  𝑆 0   1 − 𝑥 1 − 𝑥 2 − 𝑦 1 − 𝑦 2  𝑥 2 + 𝑝 𝑟 2 𝑥 2 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   1 , 𝑥  2 = 𝑘 1  1 − 𝑟 1  𝑥 1 −  𝜇 + 𝑟 2  𝑥 2 , 𝑦  1 = 𝑆 0 𝛽 2  𝑆 0   1 − 𝑥 1 − 𝑥 2 − 𝑦 1 − 𝑦 2  𝑦 2 + ( 1 − 𝑝 ) 𝑟 2 𝑥 2 −  𝜇 + 𝑘 2  𝑦 1 , 𝑦  2 = 𝑘 2 𝑦 1 − 𝜇 𝑦 2 . ( 3 . 8 ) Obviously 0 ≤ 𝑥 1 + 𝑥 2 + 𝑦 1 + 𝑦 2 ≤ 1 , ( 3 . 9 ) for all time 𝑡 ≥ 0 . For a bounded real-valued function 𝑓 on [ 0 , ∞ ) , we define 𝑓 ∞ = l i m i n f 𝑡 → ∞ 𝑓 ( 𝑡 ) , 𝑓 ∞ = l i m s u p 𝑡 → ∞ 𝑓 ( 𝑡 ) . ( 3 . 1 0 ) Lemma 3.1 (see [ 30 ]). Let 𝑓 ∶ [ 0 , ∞ ) → ∞ be bounded and continuously differentiable. Then there are sequences 𝑡 𝑛 , 𝑠 𝑛 → ∞ with the following properties: 𝑓  𝑡 𝑛  ⟶ 𝑓 ∞ , 𝑓   𝑡 𝑛  𝑓  𝑠 ⟶ 0 , 𝑛  ⟶ 𝑓 ∞ , 𝑓   𝑠 𝑛  ⟶ 0 , ( 3 . 1 1 ) for 𝑛 → ∞ . Theorem 3.2. In the absence of disease-induced death rate, the disease-free equilibrium 𝑃 0 of system ( 2.2 ) is globally asymptotically stable if 𝑅 1 < 1 , 𝑅 2 < 1 . Proof. By Lemma 3.1 and the 𝑥  2 , 𝑦  2 equations in ( 3.8 ) we have 𝑥 ∞ 2 ≤ 𝑘 1  1 − 𝑟 1  𝜇 + 𝑟 2 𝑥 ∞ 1 , 𝑦 ∞ 2 ≤ 𝑘 2 𝜇 𝑦 ∞ 1 . ( 3 . 1 2 ) Using the 𝑥  1 equation in ( 3.8 ) and choosing 𝑡 𝑛 → ∞ such that 𝑥 1  𝑡 𝑛  ⟶ 𝑥 ∞ 1 , 𝑥  1  𝑡 𝑛  ⟶ 0 , 𝑡 ⟶ ∞ , ( 3 . 1 3 ) we get 0 ≤ 𝑆 0 𝛽 1  𝑆 0   1 − 𝑥 1 − 𝑥 2 − 𝑦 1 − 𝑦 2  ∞ 𝑥 ∞ 2 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1 ≤ 𝑆 0 𝛽 1  𝑆 0  𝑥 ∞ 2 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1 ≤ 𝑆 0 𝛽 1  𝑆 0  𝑘 1  1 − 𝑟 1  𝜇 + 𝑟 2 𝑥 ∞ 1 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1 ≤  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1  𝑘 1  1 − 𝑟 1  𝑆 0 𝛽 1  𝑆 0   𝜇 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1  ≤    − 1 𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1  𝑘 1  1 − 𝑟 1  𝑆 0 𝛽 1  𝑆 0   𝜇 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1   =  − 1 𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1  𝑅 1  . − 1 ( 3 . 1 4 ) It is shown that 𝑥 ∞ 1 ≤ 0 since 𝑅 1 < 1 , as 𝑥 ∞ 1 ≥ 0 , we have that 𝑥 ∞ 1 = 0 . The inequalities in ( 3.12 ) also imply that 𝑥 ∞ 2 = 0 . Using the 𝑦  1 equation in ( 3.8 ) and choosing 𝑠 𝑛 → ∞ such that 𝑦 1  𝑠 𝑛  ⟶ 𝑦 ∞ 1 , 𝑦  1  𝑠 𝑛  ⟶ 0 , 𝑡 → ∞ , ( 3 . 1 5 ) we get 0 ≤ 𝑆 0 𝛽 2  𝑆 0   1 − 𝑥 1 − 𝑥 2 − 𝑦 1 − 𝑦 2  ∞ 𝑦 ∞ 2 −  𝜇 + 𝑘 2  𝑦 ∞ 1 ≤ 𝑆 0 𝛽 2  𝑆 0  𝑦 ∞ 2 −  𝜇 + 𝑘 2  𝑦 ∞ 1 ≤ 𝑆 0 𝛽 2  𝑆 0  𝑘 2 𝜇 𝑦 ∞ 1 −  𝜇 + 𝑘 2  𝑦 ∞ 1 ≤  𝜇 + 𝑘 2  𝑦 ∞ 1  𝑘 2 𝑆 0 𝛽 2  𝑆 0  𝜇  𝜇 + 𝑘 2   =  − 1 𝜇 + 𝑘 2  𝑦 ∞ 1  𝑅 2  . − 1 ( 3 . 1 6 ) It is shown that 𝑦 ∞ 1 ≤ 0 since 𝑅 2 < 1 , as 𝑦 ∞ 1 ≥ 0 , we have that 𝑦 ∞ 1 = 0 . The inequalities in ( 3.12 ) also imply that 𝑦 ∞ 2 = 0 . Hence l i m 𝑡 → ∞ 𝑥 1 ( 𝑡 ) = l i m 𝑡 → ∞ 𝑥 2 ( 𝑡 ) = l i m 𝑡 → ∞ 𝑦 1 ( 𝑡 ) = l i m 𝑡 → ∞ 𝑦 2 ( 𝑡 ) = 0 . ( 3 . 1 7 ) Therefore, the disease-free equilibrium 𝑃 0 is globally asymptotically stable. Next, we consider that there is disease-induced death rate, which means that 𝑑 1 > 0 , 𝑑 2 > 0 . Then, we can get the following theorem. Theorem 3.3. For 𝑑 1 > 0 , 𝑑 2 > 0 , the disease-free equilibrium 𝑃 0 of system ( 2.2 ) is globally asymptotically stable if 𝑅 1 < 1 − 𝑘 2 𝑟 2 ( 1 − 𝑝 ) / ( ( 𝜇 + 𝑑 1 + 𝑟 2 ) ( 𝜇 + 𝑘 1 ( 1 − 𝑟 1 ) ) − 𝑘 1 𝑝 𝑟 2 ( 1 − 𝑟 1 ) ) , 𝑅 2 < 1 . Proof. Let us consider the following Lyapunov function: 𝑉  𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2  = 𝑘 1  1 − 𝑟 1  𝐸 1 +  𝜇 + 𝑘 1  1 − 𝑟 1 𝐼   1 + 𝑘 2 𝐸 2 +  𝜇 + 𝑘 2  𝐼 2 . ( 3 . 1 8 ) Its time derivative along the solutions of system ( 2.2 ) satisfies ̇ 𝑉 = 𝑘 1  1 − 𝑟 1  ̇ 𝐸 1 +  𝜇 + 𝑘 1  1 − 𝑟 1 ̇ 𝐼   1 + 𝑘 2 ̇ 𝐸 2 +  𝜇 + 𝑘 2  ̇ 𝐼 2 = 𝑘 1  1 − 𝑟 1 𝛽   1 ( 𝑁 ) 𝑆 𝐼 1 + 𝑝 𝑟 2 𝐼 1 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝐸   1  +  𝜇 + 𝑘 1  1 − 𝑟 1 𝑘    1  1 − 𝑟 1  𝐸 1 −  𝜇 + 𝑑 1 + 𝑟 2  𝐼 1  + 𝑘 2  𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 1 −  𝜇 + 𝑘 2  𝐸 2  +  𝜇 + 𝑘 2 𝑘   2 𝐸 2 −  𝜇 + 𝑑 2  𝐼 2  , =  𝑘 1  1 − 𝑟 1  𝛽 1 ( 𝑁 ) 𝑆 + 𝑘 1  1 − 𝑟 1  𝑝 𝑟 2 −  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1 𝐼    1 +  𝑘 2 𝛽 2  ( 𝑁 ) 𝑆 − 𝜇 + 𝑑 2   𝜇 + 𝑘 2 𝐼   2 + 𝑘 2 𝑟 2 ( 1 − 𝑝 ) 𝐼 1 . ( 3 . 1 9 ) Now, using ( 2.3 ), one has 𝛽 1 ( 𝑁 ) 𝑆 ≤ 𝛽 1 ( 𝑆 0 ) 𝑆 0 , 𝛽 2 ( 𝑁 ) 𝑆 ≤ 𝛽 2 ( 𝑆 0 ) 𝑆 0 , then ̇ 𝑉 ≤   𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1  1 − 𝑟 1  𝑝 𝑟 2 𝑅   1  𝐼 − 1 1 + 𝑘 2 𝑟 2 ( 1 − 𝑝 ) 𝐼 1 +  𝜇 + 𝑑 2   𝜇 + 𝑘 2 𝑅   2  𝐼 − 1 2 , =    𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1  1 − 𝑟 1  𝑝 𝑟 2 𝑅   1  − 1 + 𝑘 2 𝑟 2  𝐼 ( 1 − 𝑝 ) 1 +  𝜇 + 𝑑 2   𝜇 + 𝑘 2 𝑅   2  𝐼 − 1 2 . ( 3 . 2 0 ) Since  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1   = 𝜇 𝜇 + 𝑑 1 + 𝑟 2  + 𝑘 1  1 − 𝑟 1   𝜇 + 𝑑 1  + 𝑘 1 𝑟 2  1 − 𝑟 1  − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1   = 𝜇 𝜇 + 𝑑 1 + 𝑟 2  + 𝑘 1  1 − 𝑟 1   𝜇 + 𝑑 1  + 𝑘 1 𝑟 2  1 − 𝑟 1  ( 1 − 𝑝 ) ≥ 0 , ( 3 . 2 1 ) and 𝑅 1 < 1 − 𝑘 2 𝑟 2 ( 1 − 𝑝 ) / ( ( 𝜇 + 𝑑 1 + 𝑟 2 ) ( 𝜇 + 𝑘 1 ( 1 − 𝑟 1 ) ) − 𝑘 1 𝑝 𝑟 2 ( 1 − 𝑟 1 ) ) ≤ 1 , we have   𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1  1 − 𝑟 1  𝑝 𝑟 2 𝑅   1  − 1 + 𝑘 2 𝑟 2 ( 1 − 𝑝 ) ≤ 0 . ( 3 . 2 2 ) Thus, ̇ 𝑉 ≤ 0 if 𝑅 1 < 1 − 𝑘 2 𝑟 2 ( 1 − 𝑝 ) / ( ( 𝜇 + 𝑑 1 + 𝑟 2 ) ( 𝜇 + 𝑘 1 ( 1 − 𝑟 1 ) ) − 𝑘 1 𝑝 𝑟 2 ( 1 − 𝑟 1 ) ) , 𝑅 2 < 1 . It is obvious that ̇ 𝑉 = 0 if and only if 𝐼 1 = 0 , 𝐼 2 = 0 . Then, the largest invariant set of system ( 2.2 ) on the set { ( 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) ∈ ℝ 5 + , ̇ 𝑉 ( 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) = 0 } is the disease-free equilibrium 𝑃 0 . Therefore, it follows from LaSalle's invariance principle [ 31 , 32 ] that 𝑃 0 is globally stable if 𝑅 1 < 1 − 𝑘 2 𝑟 2 ( 1 − 𝑝 ) / ( ( 𝜇 + 𝑑 1 + 𝑟 2 ) ( 𝜇 + 𝑘 1 ( 1 − 𝑟 1 ) ) − 𝑘 1 𝑝 𝑟 2 ( 1 − 𝑟 1 ) ) , 𝑅 2 < 1 . 3.2. Existence and Uniqueness of the Drug-Resistant-TB-Strain-Only Equilibrium Here, we present a result concerning the existence and uniqueness of the drug-resistant-TB-strain-only equilibrium for the system ( 2.2 ). To do this, we shall make use of the basic reproduction ratio 𝑅 1 and 𝑅 2 . Let 𝑃 1   𝐸 = ( 𝑆 , 0 , 0 , 2 ,  𝐼 2 ) be boundary equilibrium of system ( 2.2 ). Then, the boundary equilibrium can be obtained by setting the right-hand side of each differential equations in system ( 2.2 ) equal to zero with the exception of second and third equation, giving Λ − 𝛽 2   𝑁   𝑆  𝐼 2  𝛽 − 𝜇 𝑆 = 0 , 2   𝑁   𝑆  𝐼 2 −  𝜇 + 𝑘 2   𝐸 2 𝑘 = 0 , 2  𝐸 2 −  𝜇 + 𝑑 2   𝐼 2 = 0 , ( 3 . 2 3 ) adding the above three equations, we have  Λ − 𝜇 𝑁 − 𝑑 2  𝐼 2 = 0 . ( 3 . 2 4 ) Using ( 3.24 ), the first and second equation of ( 3.23 ), we can easily express  𝑆 ,  𝐸 2 and  𝐼 2 in terms of  𝑁 in the form:  𝑆 = Λ 𝑑 2 𝜇 𝑑 2 + 𝛽 2   𝑁  𝑁  ,  𝐸   Λ − 𝜇 2 = Λ 𝛽 2   𝑁  𝑁    Λ − 𝜇  𝜇 + 𝑘 2   𝜇 𝑑 2 + 𝛽 2   𝑁  𝑁 ,  𝐼   Λ − 𝜇   2 =  𝑁 Λ − 𝜇 𝑑 2 . ( 3 . 2 5 ) Substituting  𝑆 ,  𝐸 2 ,  𝐼 2 in the third equations of ( 3.23 ) yields   𝑁  𝐹   𝑁  Λ − 𝜇 = 0 , ( 3 . 2 6 ) where 𝐹   𝑁  = Λ 𝑘 2 𝑑 2 𝛽 2   𝑁  −  𝜇 + 𝑘 2   𝜇 + 𝑑 2   𝜇 𝑑 2 + 𝛽 2   𝑁  𝑁   Λ − 𝜇   . ( 3 . 2 7 ) Clearly,  Λ − 𝜇 𝑁 = 0 is a fixed point of ( 3.23 ), which corresponds to the disease-free equilibrium 𝑃 0 . Since  𝑁 ∈ [ 0 , 𝑆 0 ] , we get 𝐹 ( 0 ) = Λ 𝑘 2 𝑑 2 𝛽 2  ( 0 ) − 𝜇 + 𝑘 2   𝜇 + 𝑑 2   𝜇 𝑑 2 + Λ 𝛽 2  ( 0 ) = Λ 𝑘 2 𝑑 2 𝛽 2  𝜇 ( 0 ) − 2  𝑑 + 𝜇 2 + 𝑘 2  + 𝑑 2 𝑘 2   𝜇 𝑑 2 + Λ 𝛽 2  ( 0 ) = Λ 𝑘 2 𝑑 2 𝛽 2  𝜇 ( 0 ) − 2  𝑑 + 𝜇 2 + 𝑘 2    𝜇 𝑑 2 + Λ 𝛽 2  ( 0 ) − Λ 𝑑 2 𝑘 2 𝛽 2 ( 0 ) − 𝜇 𝑘 2 𝑑 2 2  𝜇 = − 2  𝑘 + 𝜇 2 + 𝑑 2    𝜇 𝑑 2 + Λ 𝛽 2  ( 0 ) − 𝜇 𝑘 2 𝑑 2 2 , 𝐹  𝑆 0  = Λ 𝑘 2 𝑑 2 𝛽 2  𝑆 0  − 𝜇 𝑑 2  𝜇 + 𝑘 2   𝜇 + 𝑑 2   𝑘 = 𝜇 2 𝑑 2 𝛽 2  𝑆 0  𝑆 0 − 𝑑 2  𝜇 + 𝑘 2   𝜇 + 𝑑 2   = 𝜇 𝑑 2  𝜇 + 𝑘 2   𝜇 + 𝑑 2   𝑘 2 𝛽 2  𝑆 0  𝑆 0  𝜇 + 𝑘 2   𝜇 + 𝑑 2   − 1 = 𝜇 𝑑 2  𝜇 + 𝑘 2   𝜇 + 𝑑 2 𝑅   2  . − 1 ( 3 . 2 8 ) It appears that 𝐹 ( 0 ) < 0 and 𝐹 ( 𝑆 0 ) > 0 when 𝑅 2 > 1 . The existence of fixed point follows from the intermediate value-theorem. Now,  𝐹 ( 𝑁 ) is monotone increasing, so that  𝐹 ( 𝑁 ) = 0 has only one positive root in the interval [ 0 , 𝑆 0 ] . Thus, we have established the following result Lemma 3.4. When 𝑅 2 > 1 , the system ( 2.2 ) has a unique drug-resistant-TB-strain-only equilibrium 𝑃 1   𝐸 = ( 𝑆 , 0 , 0 , 2 ,  𝐼 2 ) with  𝑆 ,  𝐸 2 , and  𝐼 2 all nonnegative. 3.3. Global Stability of the Drug-Resistant-TB-Strain-Only Equilibrium Theorem 3.5. If 𝑅 1 < 1 < 𝑅 2 , the drug-resistant-TB-strain-only equilibrium 𝑃 1 exists and is locally asymptotically stable. Proof. Reorder the equilibrium 𝑃 1 as  𝐸 ( 0 , 0 , 2 ,  𝐼 2 ,  𝑆 ) . Similarly as in the proof of locally stability of the disease-free steady state 𝑃 0 , we have 𝐹 1 =  0 𝛽 1   𝑁   𝑆  , 𝑉 0 0 1 =  𝜇 + 𝑘 1  1 − 𝑟 1  − 𝑝 𝑟 2 − 𝑘 1  1 − 𝑟 1  𝜇 + 𝑑 1 + 𝑟 2  , 𝐹 1 𝑉 1 − 1  𝑘 = 𝐴 1  1 − 𝑟 1  𝛽 1   𝑁   𝑆  𝜇 + 𝑘 1  1 − 𝑟 1 𝛽   1   𝑁   𝑆  , 0 0 ( 3 . 2 9 ) where 1 𝐴 =  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  . ( 3 . 3 0 ) The spectral radius of 𝐹 1 𝑉 1 − 1 is given by 𝜌  𝐹 1 𝑉 1 − 1  = 𝑘 1  1 − 𝑟 1  𝛽 1   𝑁   𝑆  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  , ≤ 𝑘 1  1 − 𝑟 1  𝛽 1   𝑁   𝑁  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  . ( 3 . 3 1 ) Since ( 𝛽 1 ( 𝑁 ) 𝑁 )  ≥ 0 , and  𝑁 ≤ 𝑆 0 , then 𝛽 1   𝑁   𝑁 ≤ 𝛽 1  𝑆 0  𝑆 0 , ( 3 . 3 2 ) Thus 𝜌  𝐹 1 𝑉 1 − 1  ≤ 𝑘 1  1 − 𝑟 1  𝛽 1  𝑆 0  𝑆 0  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  = 𝑅 1 . ( 3 . 3 3 ) As 𝑅 1 < 1 , we have 𝜌 ( 𝐹 1 𝑉 1 − 1 ) < 1 , which implies that the drug-resistant-TB-strain-only equilibrium 𝑃 1 is locally asymptotically stable by Theorem 3.2 in van den Driessche and Watmough [ 27 ]. The following theorem provides the global stability of drug-resistant-TB-strain-only equilibrium. Theorem 3.6. If there is no disease-induced death rate, the drug-resistant-TB-strain-only equilibrium 𝑃 1 is global asymptotically stable when 𝑅 1 < 1 < 𝑅 2 . Proof. If 𝑅 1 < 1 < 𝑅 2 , from Theorem 3.5 , 𝑃 1 is locally asymptotically stable. In the following, we only need to prove that the 𝑃 1 is a global attractor. There is no disease-induced death rate, which means 𝑑 1 = 𝑑 2 = 0 , using a similar argument as in the proof of Theorem 3.2 , we can show that if 𝑅 1 < 1 , l i m 𝑡 → ∞ 𝑥 1 ( 𝑡 ) = l i m 𝑡 → ∞ 𝑥 2 ( 𝑡 ) = 0 . ( 3 . 3 4 ) Then the equivalent limiting system of ( 3.6 ) is 𝐸  2 = 𝛽 2  𝑆 0  𝐼 2  𝑆 0 − 𝐸 2 − 𝐼 2  −  𝜇 + 𝑘 2  𝐸 2 , 𝐼  2 = 𝑘 2 𝐸 2 − 𝜇 𝐼 2 . ( 3 . 3 5 ) Let ( 𝑃 , 𝑄 ) be the vector field defined by system ( 3.35 ). Then for the Dulac function 𝐷 ( 𝐸 2 , 𝐼 2 ) = 1 / 𝐸 2 𝐼 2 , there holds 𝜕 𝐸 2 𝜕 ( 𝐷 𝑃 ) + 𝐼 2 𝛽 ( 𝐷 𝑄 ) = − 2  𝑆 0 𝑆   0 − 𝐼 2  𝐸 2 2 − 𝑘 2 𝐼 2 2 < 0 . ( 3 . 3 6 ) From Dulac's criterion, we can see equivalent limiting system ( 3.35 ) does not have a limit cycle. Therefore, the local stability of 𝑃 1 implies the global stability in Int Γ . This completes the proof of Theorem 3.6 . 3.4. Existence and Uniqueness of the Interior Equilibrium First of all, we assume that (W 1 ) l i m 𝑁 → 0 ( 1 / 𝛽 1 ( 𝑁 ) ) = + ∞ , (W 2 ) 𝛽 1 ( 𝑁 ) / 𝛽 2 ( 𝑁 ) = 𝛽 1 / 𝛽 2 , then, existence, uniqueness and global stability of interior equilibrium can be obtained. Let 𝑃 2 = ( 𝑆 ∗ , 𝐸 ∗ 1 , 𝐼 ∗ 1 , 𝐸 ∗ 2 , 𝐼 ∗ 2 ) be interior equilibrium of system ( 2.2 ). Then, the interior equilibrium (steady state with 𝐼 1 > 0 , 𝐼 2 > 0 ) can be obtained by setting the right-hand side of each of the five differential equations in system ( 2.2 ) equal to zero, giving Λ − 𝛽 1  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1 − 𝛽 2  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 2 − 𝜇 𝑆 ∗ 𝛽 = 0 , 1  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1 + 𝑝 𝑟 2 𝐼 ∗ 1 − 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 − 𝜇 𝐸 ∗ 1 𝑘 = 0 , 1  1 − 𝑟 1  𝐸 ∗ 1 −  𝑟 2 + 𝜇 + 𝑑 1  𝐼 ∗ 1 𝛽 = 0 , 2  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 −  𝜇 + 𝑘 2  𝐸 ∗ 2 𝑘 = 0 , 2 𝐸 ∗ 2 −  𝜇 + 𝑑 2  𝐼 ∗ 2 = 0 , ( 3 . 3 7 ) adding the above five equations, we have Λ − 𝜇 𝑁 ∗ − 𝑑 1 𝐼 ∗ 1 − 𝑑 2 𝐼 ∗ 2 = 0 . ( 3 . 3 8 ) Using ( 3.37 ) and ( 3.38 ), we get the following equivalent equations set Λ − 𝛽 1  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1 − 𝛽 2  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 2 − 𝜇 𝑆 ∗ 𝛽 = 0 , 1  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1 + 𝑝 𝑟 2 𝐼 ∗ 1 − 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 − 𝜇 𝐸 ∗ 1 𝑘 = 0 , 1  1 − 𝑟 1  𝐸 ∗ 1 −  𝑟 2 + 𝜇 + 𝑑 1  𝐼 ∗ 1 𝑘 = 0 , 2 𝐸 ∗ 2 −  𝜇 + 𝑑 2  𝐼 ∗ 2 = 0 , Λ − 𝜇 𝑁 ∗ − 𝑑 1 𝐼 ∗ 1 − 𝑑 2 𝐼 ∗ 2 = 0 . ( 3 . 3 9 ) From ( 3.39 ), we have 𝑆 ∗ = 𝛽  𝑆 0  𝑆 0 𝑅 1 𝛽 1 ( 𝑁 ∗ ) , 𝐼 ∗ 1 = 𝑑 2  Λ − 𝜇 𝑆 ∗  − 𝑆 ∗ 𝛽 2  𝑁 ∗   Λ − 𝜇 𝑁 ∗  𝑑 2 𝑆 ∗ 𝛽 1 ( 𝑁 ∗ ) − 𝑑 1 𝑆 ∗ 𝛽 2 ( 𝑁 ∗ ) , 𝐸 ∗ 1 = 𝜇 + 𝑑 1 + 𝑟 2 𝑘 1  1 − 𝑟 1  𝐼 ∗ 1 , 𝐼 ∗ 2 =  Λ − 𝜇 𝑁 ∗  𝑑 2 − 𝑑 1 𝑑 2 𝐼 ∗ 1 , 𝐸 ∗ 2 =  𝜇 + 𝑑 2   Λ − 𝜇 𝑁 ∗  𝑘 2 𝑑 2 − 𝑑 1  𝜇 + 𝑑 2  𝑘 2 𝑑 2 𝐼 ∗ 1 . ( 3 . 4 0 ) Substituting ( 3.40 ) into 𝛽 2 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 − ( 𝜇 + 𝑘 2 ) 𝐸 ∗ 2 = 0 , 𝑁 ∗ satisfies the following equation: 𝐺  𝑁 ∗  = 0 , ( 3 . 4 1 ) where 𝐺  𝑁 ∗  = 𝛽 1  𝑆 0  𝑆 0  Λ − 𝜇 𝑁 ∗  𝑅 1  𝑘 2 𝑑 2 𝑆 0 𝛽 1  𝑆 0  𝛽 2  𝑁 ∗  𝑅 1 𝛽 1 ( 𝑁 ∗ ) − 𝑑 2 𝜇 1 − ( 1 − 𝑝 ) 𝑘 2 𝑑 2 𝑟 2 𝛽 2  𝑁 ∗  𝛽 1 ( 𝑁 ∗ )  + 𝑑 2 Λ  𝛽 1 − 1  𝑆 0  𝛽 1 ( 𝑁 ∗ ) 𝑅 1 𝑑   1 𝜇 1 + ( 1 − 𝑝 ) 𝑘 2 𝑑 2 𝑟 2 − 𝑑 1 𝑘 2 𝑆 0 𝛽 1  𝑆 0  𝛽 2  𝑁 ∗  𝑅 1 𝛽 1 ( 𝑁 ∗ )  , 𝜇 1 =  𝜇 + 𝑘 2   𝜇 + 𝑑 2  . ( 3 . 4 2 ) Since 𝛽 1 ( 𝑁 ) / 𝛽 2 ( 𝑁 ) = 𝛽 1 / 𝛽 2 , l i m 𝑁 → 0 ( 1 / 𝛽 1 ( 𝑁 ) ) = + ∞ , if 𝑅 1 > 𝑑 1 𝑘 2 𝛽 2 𝑆 0 𝛽 1 ( 𝑆 0 ) / 𝛽 1 [ 𝑑 1 𝜇 1 + ( 1 − 𝑝 ) 𝑘 2 𝑑 2 𝑟 2 ] , then l i m 𝑁 ∗ → 0 𝐺  𝑁 ∗  = − ∞ . ( 3 . 4 3 ) It means that 𝐺 ( 𝑁 ∗ ) = 0 has only one positive root in the interval [ 0 , 𝑆 0 ] if 𝐺 ( 𝑁 ∗ ) is monotone increasing and 𝐺 ( 𝑆 0 ) ≥ 0 , which is equivalent to 𝑅 1 > 𝐻 1 + 𝐻 2 𝐻 3 + 𝐻 4 , 𝑅 2 <  1 + ( 1 − 𝑝 ) 𝑑 2 𝑘 2 𝑟 2 𝑑 1 𝜇 1  𝑅 1 , ( 3 . 4 4 ) where 𝐻 1 = 𝑘 2 𝑑 2 𝛽 2 𝑆 0 𝛽 1  𝑆 0  𝛽 2 1  𝑁 ∗  , 𝐻 2 = 𝑑 1 𝑑 2 𝑘 2 𝛽 2 𝑆 0 𝛽 1  𝑆 0  𝛽  1  𝑁 ∗  , 𝐻 3 = 𝛽 2 1  𝑁 ∗ 𝑑   2 𝛽 1 𝜇 1 + 𝑑 2 𝑘 2 𝛽 2 𝑟 2  , 𝐻 ( 1 − 𝑝 ) 4 = 𝑑 2 𝛽  1  𝑁 ∗ 𝑑   1 𝛽 1 𝜇 1 + 𝑑 2 𝑘 2 𝛽 1 𝑟 2  . ( 1 − 𝑝 ) ( 3 . 4 5 ) Let  𝐻 𝐻 = m a x 1 , 1 + 𝐻 2 𝐻 3 + 𝐻 4 , 𝑑 1 𝑘 2 𝛽 2 𝑆 0 𝛽 1  𝑆 0  𝛽 1  𝑑 1 𝜇 1 + 𝑘 2 𝑑 2 𝑟 2   , 𝑓  𝑅 ( 1 − 𝑝 ) 1  =  1 + ( 1 − 𝑝 ) 𝑑 2 𝑘 2 𝑟 2 𝑑 1 𝜇 1  𝑅 1 . ( 3 . 4 6 ) When 𝑅 1 > 𝐻 , 𝑅 2  𝑅 < 𝑓 1  , ( 3 . 4 7 ) it shows that 𝐺 ( 𝑆 0 ) ≥ 0 , and 𝐺 ( 𝑁 ∗ ) is monotone increasing, by the intermediate value-theorem, 𝐺 ( 𝑁 ∗ ) = 0 has only one positive root in the interval [ 0 , 𝑆 0 ] , which signify that 𝑃 2 is unique. Obviously, 𝑆 ∗ > 0 . In order to make certain that 𝐸 ∗ 1 > 0 , 𝐼 ∗ 1 > 0 , 𝐸 ∗ 2 > 0 , 𝐼 ∗ 2 > 0 , we obtain that 𝐸 ∗ 1 > 0 , 𝐼 ∗ 1 > 0 , 𝐸 ∗ 2 > 0 , 𝐼 ∗ 2 > 0 if and only if 0 < 𝐼 ∗ 1 < Λ − 𝜇 𝑁 ∗ 𝑑 1 . ( 3 . 4 8 ) Due to 𝑑 2 > 𝑑 1 > 0 , introducing ( 3.40 ) into ( 3.48 ), we acquire 𝐵 1 < 𝑅 1 < 𝐵 2 , ( 3 . 4 9 ) where 𝐵 1 = 𝑆 0 𝛽 1  𝑆 0 𝛽   2  𝑁 ∗   Λ − 𝜇 𝑁 ∗  + 𝜇 𝑑 2  Λ 𝑑 2 𝛽 1 ( 𝑁 ∗ ) , 𝐵 2 = 𝑆 0 𝛽 1  𝑆 0 𝛽   1  𝑁 ∗   Λ − 𝜇 𝑁 ∗  + 𝜇 𝑑 1  Λ 𝑑 1 𝛽 1 ( 𝑁 ∗ ) , ( 3 . 5 0 ) and 𝐵 1 < 𝐵 2 . Now, we can get the following conclusion. Lemma 3.7. When conditions (W 1 ), (W 2 ), ( 3.47 ) and ( 3.49 ) are satisfied, the system ( 2.2 ) has a unique interior equilibrium 𝑃 2 = ( 𝑆 ∗ , 𝐸 ∗ 1 , 𝐼 ∗ 1 , 𝐸 ∗ 2 , 𝐼 ∗ 2 ) with 𝑆 ∗ , 𝐸 ∗ 1 , 𝐼 ∗ 1 , 𝐸 ∗ 2 , and 𝐼 ∗ 2 all nonnegative. Remark 3.8. For 𝛽 1 ( 𝑁 ) = 𝛽 1 and 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , condition (W 1 ) is not satisfied. However, at least when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , condition (W 1 ) is satisfied, we also can find our conclusions are still correct from Figures 7 and 10 . 3.5. Global Stability of the Interior Equilibrium Herein, we study the global stability of the interior equilibrium 𝑃 2 of system ( 2.2 ). For simplicity, we only consider the special case when the transmission coefficients of drug-sensitive TB-strain and drug-resistant TB-strain are equal, that is, 𝛽 1 ( 𝑁 ) = 𝛽 2 ( 𝑁 ) = 𝛽 ( 𝑁 ) , then, we get the following result. Theorem 3.9. If conditions (W 1 ), (W 2 ), ( 3.47 ) and ( 3.49 ) are satisfied, interior equilibrium 𝑃 2 of system ( 2.2 ) is is globally asymptotically stable when 𝑆 𝑆 ∗ = 𝐸 1 𝐸 ∗ 1 = 𝐸 2 𝐸 ∗ 2 𝐼 ≥ 1 , 2 𝐼 ∗ 2 ≥ 𝐼 1 𝐼 ∗ 1 𝑆 o r 𝑆 ∗ = 𝐸 1 𝐸 ∗ 1 = 𝐸 2 𝐸 ∗ 2 𝐼 ≤ 1 , 2 𝐼 ∗ 2 ≤ 𝐼 1 𝐼 ∗ 1 . ( 3 . 5 1 ) Proof. Following [ 33 ], we consider the Lyapunov function: 𝑈  𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2  =  𝑆 − 𝑆 ∗  +  𝐸 l n 𝑆 1 − 𝐸 ∗ 1 l n 𝐸 ∗ 1   𝐼 + 𝐶 1 − 𝐼 ∗ 1 l n 𝐼 ∗ 1  +  𝐸 2 − 𝐸 ∗ 2 l n 𝐸 ∗ 2   𝐼 + 𝐷 2 − 𝐼 ∗ 2 l n 𝐼 ∗ 2  , ( 3 . 5 2 ) where 𝐶 and 𝐷 are positive constants to be determined later. Differentiating this function with respect to time along the solutions of system ( 2.2 ) yields ̇  𝑆 𝑈 = 1 − ∗ 𝑆  ̇  𝐸 𝑆 + 1 − ∗ 1 𝐸 1  ̇ 𝐸 1  𝐼 + 𝐶 1 − ∗ 1 𝐼 1  ̇ 𝐼 1 +  𝐸 1 − ∗ 2 𝐸 2  ̇ 𝐸 2  𝐼 + 𝐷 1 − ∗ 2 𝐼 2  ̇ 𝐼 2 =  𝑆 1 − ∗ 𝑆   Λ − 𝛽 ( 𝑁 ) 𝑆 𝐼 1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 2  +  𝐸 − 𝜇 𝑆 1 − ∗ 1 𝐸 1   𝛽 ( 𝑁 ) 𝑆 𝐼 1 + 𝑝 𝑟 2 𝐼 1 − 𝑘 1  1 − 𝑟 1  𝐸 1 − 𝜇 𝐸 1   𝐼 + 𝐶 1 − ∗ 1 𝐼 1   𝑘 1  1 − 𝑟 1  𝐸 1 −  𝑟 2 + 𝜇 + 𝑑 1  𝐼 1  +  𝐸 1 − ∗ 2 𝐸 2   𝛽 ( 𝑁 ) 𝑆 𝐼 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 1 −  𝜇 + 𝑘 2  𝐸 2   𝐼 + 𝐷 1 − ∗ 2 𝐼 2   𝑘 2 𝐸 2 −  𝜇 + 𝑑 2  𝐼 2  . ( 3 . 5 3 ) Considering ( 3.37 ), we can deduce that  𝑁 Λ = 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 + 𝜇 𝑆 ∗ , 𝜇 + 𝑘 1  1 − 𝑟 1   𝑁 = 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1 𝐸 ∗ 1 + 𝑝 𝑟 2 𝐼 ∗ 1 𝐸 ∗ 1 , 𝜇 + 𝑑 1 + 𝑟 2 = 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 𝐼 ∗ 1 , 𝜇 + 𝑘 2  𝑁 = 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 𝐸 ∗ 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 𝐸 ∗ 2 , 𝜇 + 𝑑 2 = 𝑘 2 𝐸 ∗ 2 𝐼 ∗ 2 . ( 3 . 5 4 ) Then, ( 3.53 ) becomes ̇  𝑆 𝑈 = 1 − ∗ 𝑆   𝛽  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 + 𝜇 𝑆 ∗ − 𝛽 ( 𝑁 ) 𝑆 𝐼 1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 2  +  𝐸 − 𝜇 𝑆 1 − ∗ 1 𝐸 1   𝛽 ( 𝑁 ) 𝑆 𝐼 1 + 𝑝 𝑟 2 𝐼 1  𝑁 − 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1 𝐸 1 𝐸 ∗ 1 − 𝑝 𝑟 2 𝐼 ∗ 1 𝐸 1 𝐸 ∗ 1   𝐼 + 𝐶 1 − ∗ 1 𝐼 1 𝑘   1  1 − 𝑟 1  𝐸 1 − 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 𝐼 1 𝐼 ∗ 1  +  𝐸 1 − ∗ 2 𝐸 2   𝛽 ( 𝑁 ) 𝑆 𝐼 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 1  𝑁 − 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 𝐸 2 𝐸 ∗ 2 − ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 𝐸 2 𝐸 ∗ 2   𝐼 + 𝐷 1 − ∗ 2 𝐼 2 𝑘   2 𝐸 2 − 𝑘 2 𝐸 ∗ 2 𝐼 2 𝐼 ∗ 2  𝜇  = − 𝑆 − 𝑆 ∗  2 𝑆  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1  1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 1 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 1 𝑆   1 − ∗ 𝑆   𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2  1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 2 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 𝑆   1 − ∗ 𝑆  +  𝐸 1 − ∗ 1 𝐸 1 𝛽  𝑁   ∗  𝑆 ∗ 𝐼 ∗ 1  𝛽 ( 𝑁 ) 𝑆 𝐼 1 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 1 − 𝐸 1 𝐸 ∗ 1  + 𝑝 𝑟 2 𝐼 ∗ 1  𝐼 1 𝐼 ∗ 1 − 𝐸 1 𝐸 ∗ 1  𝐼   + 𝐶 1 − ∗ 1 𝐼 1 𝑘   1  1 − 𝑟 1  𝐸 ∗ 1  𝐸 1 𝐸 ∗ 1 − 𝐼 1 𝐼 ∗ 1 +  𝐸   1 − ∗ 2 𝐸 2 𝛽  𝑁   ∗  𝑆 ∗ 𝐼 ∗ 2  𝛽 ( 𝑁 ) 𝑆 𝐼 2 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 − 𝐸 2 𝐸 ∗ 2  + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  𝐼 1 𝐼 ∗ 1 − 𝐸 2 𝐸 ∗ 2  𝐼   + 𝐷 1 − ∗ 2 𝐼 2 𝑘   2 𝐸 ∗ 2  𝐸 2 𝐸 ∗ 2 − 𝐼 2 𝐼 ∗ 2 .   ( 3 . 5 5 ) Now, using ( 3.37 ), we have 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 =  𝑟 2 + 𝜇 + 𝑑 1  𝐼 ∗ 1 , 𝑘 2 𝐸 ∗ 2 =  𝜇 + 𝑑 2  𝐼 ∗ 2 . ( 3 . 5 6 ) Then, ( 3.55 ) can be rewritten as follows: ̇ 𝜇  𝑈 = − 𝑆 − 𝑆 ∗  2 𝑆  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1   1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 1 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 1 𝑆   1 − ∗ 𝑆  +  𝐸 1 − ∗ 1 𝐸 1   𝛽 ( 𝑁 ) 𝑆 𝐼 1 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 1 − 𝐸 1 𝐸 ∗ 1   + 𝑝 𝑟 2 𝐼 ∗ 1   𝐸 1 − ∗ 1 𝐸 1 𝐼   1 𝐼 ∗ 1 − 𝐸 1 𝐸 ∗ 1  + 𝐶  𝜇 + 𝑟 2 + 𝑑 1  𝑝 𝑟 2  𝐼 1 − ∗ 1 𝐼 1 𝐸   1 𝐸 ∗ 1 − 𝐼 1 𝐼 ∗ 1    𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2   1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 2 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 𝑆   1 − ∗ 𝑆  +  𝛽 ( 𝑁 ) 𝑆 𝐼 2 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 − 𝐸 2 𝐸 ∗ 2 𝐸   1 − ∗ 2 𝐸 2   + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1   𝐼 1 𝐼 ∗ 1 − 𝐸 2 𝐸 ∗ 2 𝐸   1 − ∗ 2 𝐸 2  + 𝐷  𝜇 + 𝑑 2  𝐼 ∗ 2 ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  𝐼 1 − ∗ 2 𝐼 2 𝐸   2 𝐸 ∗ 2 − 𝐼 2 𝐼 ∗ 2   . ( 3 . 5 7 ) Now, let ( 𝑥 , 𝑤 1 , 𝑧 1 , 𝑤 2 , 𝑧 2 ) = ( 𝑆 / 𝑆 ∗ , 𝐸 1 / 𝐸 ∗ 1 , 𝐼 1 / 𝐼 ∗ 1 , 𝐸 2 / 𝐸 ∗ 2 , 𝐼 2 / 𝐼 ∗ 2 ) and 𝑔 ( 𝑁 ) = 𝛽 ( 𝑁 ) / 𝛽 ( 𝑁 ∗ ) , then we get ̇ 𝜇  𝑈 = − 𝑆 − 𝑆 ∗  2 𝑆  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1   1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 1   1 1 − 𝑥  +  1 1 − 𝑤 1   𝑔 ( 𝑁 ) 𝑥 𝑧 1 − 𝑤 1   + 𝑝 𝑟 2 𝐼 ∗ 1   1 1 − 𝑤 1   𝑧 1 − 𝑤 1  + 𝐶  𝜇 + 𝑟 2 + 𝑑 1  𝑝 𝑟 2  1 1 − 𝑧 1   𝑤 1 − 𝑧 1    𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2   1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 2   1 1 − 𝑥  +  𝑔 ( 𝑁 ) 𝑥 𝑧 2 − 𝑤 2   1 1 − 𝑤 2   + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1   𝑧 1 − 𝑤 2   1 1 − 𝑤 2  + 𝐷  𝜇 + 𝑑 2  𝐼 ∗ 2 ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  1 1 − 𝑧 2   𝑤 2 − 𝑧 2   , 𝜇  = − 𝑆 − 𝑆 ∗  2 𝑆  + 𝑓 𝑥 , 𝑤 1 , 𝑧 1 , 𝑤 2 , 𝑧 2  , ( 3 . 5 8 ) where 𝑓  𝑥 , 𝑤 1 , 𝑧 1 , 𝑤 2 , 𝑧 2   𝑁 = 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1 𝑓 1  𝑥 , 𝑤 1 , 𝑧 1  + 𝑝 𝑟 2 𝐼 ∗ 1 𝑓 2  𝑤 1 , 𝑧 1   𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 𝑓 3  𝑥 , 𝑤 2 , 𝑧 2  + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 𝑓 4  𝑧 1 , 𝑤 2 , 𝑧 2  , 𝑓 1  𝑥 , 𝑤 1 , 𝑧 1  =  1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 1   1 1 − 𝑥  +  1 1 − 𝑤 1   𝑔 ( 𝑁 ) 𝑥 𝑧 1 − 𝑤 1  , 𝑓 2  𝑤 1 , 𝑧 1  =  1 1 − 𝑤 1   𝑧 1 − 𝑤 1  + 𝐶  𝜇 + 𝑟 2 + 𝑑 1  𝑝 𝑟 2  1 1 − 𝑧 1   𝑤 1 − 𝑧 1  , 𝑓 3  𝑥 , 𝑤 2 , 𝑧 2  =  1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 2   1 1 − 𝑥  +  𝑔 ( 𝑁 ) 𝑥 𝑧 2 − 𝑤 2   1 1 − 𝑤 2  , 𝑓 4  𝑧 1 , 𝑤 2 , 𝑧 2  =  𝑧 1 − 𝑤 2   1 1 − 𝑤 2  + 𝐷  𝜇 + 𝑑 2  𝐼 ∗ 2 ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  1 1 − 𝑧 2   𝑤 2 − 𝑧 2  . ( 3 . 5 9 ) We can choose 𝐶 = 𝑝 𝑟 2  𝜇 + 𝑟 2 + 𝑑 1  , 𝐷 = ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  𝜇 + 𝑑 2  𝐼 ∗ 2 . ( 3 . 6 0 ) For 𝑥 = 𝑤 1 = 𝑤 2 , we have 𝑓 1  𝑥 , 𝑤 1 , 𝑧 1  =  1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 1   1 1 − 𝑥  +  1 1 − 𝑤 1   𝑔 ( 𝑁 ) 𝑥 𝑧 1 − 𝑤 1  = 𝑔 ( 𝑁 ) 𝑥 𝑧 1  1 𝑥 1 − 1 + 1 − 𝑤 1  1 + 1 − 𝑥 − 𝑤 1  1 1 − 𝑤 1  1 = 1 − 𝑥 − 𝑤 1 1 + 1 = 2 − 𝑥 − 𝑥 , ( 3 . 6 1 ) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑥 = 1 𝑓 2  𝑤 1 , 𝑧 1  =  1 1 − 𝑤 1   𝑧 1 − 𝑤 1  +  1 1 − 𝑧 1   𝑤 1 − 𝑧 1  =  𝑧 1 − 𝑤 1   1 1 − 𝑤 1 1 − 1 + 𝑧 1  =  𝑧 1 − 𝑤 1   1 𝑧 1 − 1 𝑤 1  𝑤 = 2 − 1 𝑧 1 − 𝑧 1 𝑤 1 , ( 3 . 6 2 ) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑤 1 = 𝑧 1 𝑓 3  𝑥 , 𝑤 2 , 𝑧 2  =  1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 2   1 1 − 𝑥  +  𝑔 ( 𝑁 ) 𝑥 𝑧 2 − 𝑤 2   1 1 − 𝑤 2  = 𝑔 ( 𝑁 ) 𝑥 𝑧 2  1 1 − 𝑤 2 + 1 𝑥  1 − 1 + 1 − 𝑥 − 𝑤 2 1 + 1 = 1 − 𝑥 − 𝑤 2 1 + 1 = 2 − 𝑥 − 𝑥 , ( 3 . 6 3 ) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑥 = 1 . For 𝑤 2 ≥ 1 , 𝑧 2 ≥ 𝑧 1 or 𝑤 2 ≤ 1 , 𝑧 2 ≤ 𝑧 1 , we have 𝑓 4  𝑧 1 , 𝑤 2 , 𝑧 2  =  𝑧 1 − 𝑤 2   1 1 − 𝑤 2  +  1 1 − 𝑧 2   𝑤 2 − 𝑧 2  ≤  𝑧 2 − 𝑤 2   1 1 − 𝑤 2 1 − 1 + 𝑧 2  =  𝑧 2 − 𝑤 2   1 𝑧 2 − 1 𝑤 2  𝑤 = 2 − 2 𝑧 2 − 𝑧 2 𝑤 2 , ( 3 . 6 4 ) which is less than or equal to zero by the arithmetic-geometric mean inequality, with equality if and only if 𝑤 2 = 𝑧 2 . Thus, ̇ 𝑈 ( 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) is less or equal to zero with equality only if 𝑥 = 1 , 𝑤 1 = 𝑧 1 , 𝑤 2 = 𝑧 2 . Since 𝑥 = 𝑤 1 = 𝑤 2 , the largest invariant set of system ( 2.2 ) on the set { ( 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) ∈ ℝ 5 + , ̇ 𝑈 ( 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) = 0 } is the endemic equilibrium point 𝑃 2 . By LaSalle's invariance principle [ 31 , 32 ], we can conclude that the interior equilibrium 𝑃 2 is globally asymptotically stable when condition ( 3.51 ) is satisfied. Remark 3.10. It is possible for condition ( 3.51 ) to fail, in which case the global stability of the interior equilibrium of system ( 2.2 ) has not been established. Figures 7 , 10 and 13 , however, seem to support the idea that the interior equilibrium of system ( 2.2 ) is still global asymptotically stable even in this case. 4. Numerical Simulations To illustrate the theoretical results contained in this paper, we give some simulations using the parameter value in Table 1 . Numerical results are displayed in the following figures. Table 1: Description and estimation of parameters. Taking no account of disease-induced death rate, we assume that 𝑑 1 = 𝑑 2 = 0 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 0 . 0 0 0 1 (so that 𝑅 1 < 1 , 𝑅 2 < 1 ). For 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 and 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑏 𝑁 / ( 1 + 𝑏 𝑁 + 1 + 2 𝑏 𝑁 ) , 𝑏 = 1 , Figures 2 , 3 , and 4 show the trajectories plots and their plane figures, respectively, which are in agreement with Theorem 3.2 . Figure 2: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 0 . 0 0 0 1 , 𝑑 1 = 𝑑 2 = 0 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 3: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 0 . 0 0 0 1 , 𝑑 1 = 𝑑 2 = 0 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 4: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑁 / ( 1 + 𝑁 + 1 + 2 𝑁 ) , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 0 . 0 0 0 1 , 𝑑 1 = 𝑑 2 = 0 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. When there is disease-induced death rate, which means that 𝑑 1 > 0 , 𝑑 2 > 0 , we take 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 . Firstly, we consider system ( 2.2 ) with 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 and choose 𝛽 1 = 𝛽 2 = 0 . 0 0 1 , so that 𝑅 1 < 1 , 𝑅 2 < 1 , Figure 5 presents the trajectories plot and plane figure. From this figure, we can see that the trajectories of system ( 2.2 ) converge to the disease-free equilibrium. We also choose 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 , so that 𝑅 1 < 1 , 𝑅 2 > 1 , Figure 6 presents the trajectories plot and plane figure. It shows that the drug-resistant-TB-strain-only equilibrium is global asymptotically stable. Then, we choose 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 , so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) , the trajectories plot and its plane figure are depicted in Figure 7 , we can observe that the trajectories converge to the unique interior equilibrium. This implies that both the sensitive and resistant strains coexist and prevail in the host population as shown in Theorem 3.9 . Figure 5: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 = 𝛽 2 = 0 . 0 0 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 6: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 > 1 . It shows that the boundary equilibrium 𝑃 1 is stable. Figure 7: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) . It shows that the interior equilibrium 𝑃 2 is stable. Next, we need to think the case when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 . Figures 8 , 9 and 10 show the trajectories plot and its plane figure of system ( 2.2 ) when 𝛽 1 = 𝛽 2 = 0 . 0 0 1 (so that 𝑅 1 < 1 , 𝑅 2 < 1 ), 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 (so that 𝑅 1 < 1 , 𝑅 2 > 1 ), and 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 (so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) ), respectively. As the previous cases, we can find that the trajectories of system ( 2.2 ) converge to the disease-free equilibrium when 𝑅 1 < 1 , 𝑅 2 < 1 (see Figure 8 ), the drug-resistant-TB-strain-only equilibrium is global asymptotically stable when 𝑅 1 < 1 , 𝑅 2 > 1 (see Figure 9 ), while the trajectories of system ( 2.2 ) converge to the unique interior equilibrium when 𝑅 1 > 1 , 𝑅 2 < 𝑓 ( 𝑅 1 ) (see Figure 10 ). Figure 8: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 , 𝛽 1 = 𝛽 2 = 0 . 0 0 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 9: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 > 1 . It shows that the boundary equilibrium 𝑃 1 is stable. Figure 10: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 , 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) . It shows that the interior equilibrium 𝑃 2 is stable. Eventually, we assume that 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑏 𝑁 / ( 1 + 𝑏 𝑁 + 1 + 2 𝑏 𝑁 ) in system ( 2.2 ), of which 𝑏 = 1 . Numerical results are depicted in Figures 11 , 12 , and 13 . Figure 11 presents the trajectories and its plane figure when 𝛽 1 = 𝛽 2 = 0 . 0 0 1 (so that 𝑅 1 < 1 , 𝑅 2 < 1 ). From this figure, it clearly appears that the disease disappears in the host population. Figure 12 presents the trajectories and its plane figure when 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 (so that 𝑅 1 < 1 , 𝑅 2 > 1 ), we can see that the drug-resistant-TB-strain-only equilibrium is global asymptotically stable. Figure 13 gives a description of the trajectories and its plane figure when 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 (so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) ). By observation, we conclude that both the sensitive and resistant strains will coexist and prevail in the host population. Figure 11: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑁 / ( 1 + 𝑁 + 1 + 2 𝑁 ) , 𝛽 1 = 𝛽 2 = 0 . 0 0 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 12: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑁 / ( 1 + 𝑁 + 1 + 2 𝑁 ) , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 > 1 . It shows that the boundary equilibrium 𝑃 1 is stable. Figure 13: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑁 / ( 1 + 𝑁 + 1 + 2 𝑁 ) , 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) . It shows that the interior equilibrium 𝑃 2 is stable. 5. Discussion In this paper, we studied a two-strain tuberculosis model that includes both drug-sensitive and drug-resistant strains. The novelty of our model is that we incorporate general contact rate 𝛽 𝑖 ( 𝑁 ) , 𝑖 = 1 , 2 into a two-strain tuberculosis model. By using the Lyapunov stability theory, LaSalle's invariant set theorem and Dulac's criterion, the global stability of equilibria of the proposed model is proved. The basic reproduction number 𝑅 1 , 𝑅 2 provide the threshold conditions which determine the competitive outcomes of the two strains. When 𝑅 1 < 1 , 𝑅 2 < 1 , the disease-free equilibrium 𝑃 0 is globally stable and the disease will die out. When 𝑅 1 < 1 , 𝑅 2 > 1 , a unique boundary equilibrium 𝑃 1 exists and is globally stable. While when 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) , there exists a unique interior equilibrium 𝑃 2 which is globally asymptotically stable under certain conditions. A fairly good agreement is obtained between the analytical and numerical results. With the appearance of numerous second-line antituberculous drugs, such as capreomycin, amikacin, and kanamycin, the treatment of drug-resistant-TB has become possible. We will leave this for future study. Acknowledgments This work was supported by the NNSF of China (10961018), the Key Project of Chinese Ministry of Education (209131), the Project Sponsored by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, the NSF of Gansu Province of China (3ZS042-B25-013), the NSF of Bureau of Education of Gansu Province of China (0803-01) and the Development Program for Outstanding Young Teachers in Lanzhou University of Technology (Q200703 ) and the Doctor's Foundation of Lanzhou University of Technology. <h4>References</h4> C. Dye, S. Scheele, P. Dolin, V. Pathania, and M. C. Raviglione, “ Global burden of tuberculosis: estimated incidence, prevalence, and mortality by country ,” Journal of the American Medical Association , vol. 282, no. 7, pp. 677–686, 1999. T. R. Frieden, T. R. Sterling, S. S. Munsiff, C. J. Watt, and C. Dye, “ Tuberculosis ,” Lancet , vol. 362, no. 9387, pp. 887–899, 2003. B. Song, C. Castillo-Chavez, and J. P. Aparicio, “ Tuberculosis models with fast and slow dynamics: the role of close and casual contacts ,” Mathematical Biosciences , vol. 180, pp. 187–205, 2002. B. Miller, “Preventive therapy for tuberculosis,” Medical Clinics of North America , vol. 77, no. 6, pp. 1263–1275, 1993. P. Rodrigues, M. G. M. Gomes, and C. Rebelo, “ Drug resistance in tuberculosis-a reinfection model ,” Theoretical Population Biology , vol. 71, no. 2, pp. 196–212, 2007. M. A. Espinal, et al., “Standard short-course chemotherapy for drug-resistant tuberculosis: treatment outcomes in 6 countries,” Journal of the American Medical Association , vol. 283, no. 19, pp. 2537–2545, 2000. H. W. Hethcote, “ The mathematics of infectious diseases ,” SIAM Review , vol. 42, no. 4, pp. 599–653, 2000. F. Brauer and C. Castillo-Chávez, Mathematical models in population biology and epidemiology , vol. 40 of Texts in Applied Mathematics , Springer, New York, NY, USA, 2001. B. M. Murphy, B. H. Singer, S. Anderson, and D. Kirschner, “ Comparing epidemic tuberculosis in demographically distinct heterogeneous populations ,” Mathematical Biosciences , vol. 180, pp. 161–185, 2002. C. Castillo-Chavez and B. Song, “Dynamical models of tuberculosis and their applications,” Mathematical Biosciences and Engineering , vol. 1, no. 2, pp. 361–404, 2004. B. M. Murphy, B. H. Singer, and D. Kirschner, “ On treatment of tuberculosis in heterogeneous populations ,” Journal of Theoretical Biology , vol. 223, no. 4, pp. 391–404, 2003. C. Castillo-Chavez and Z. Feng, “ Global stability of an age-structure model for TB and its applications to optimal vaccination strategies ,” Mathematical Biosciences , vol. 151, no. 2, pp. 135–154, 1998. C. C. McCluskey, “Lyapunov functions for tuberculosis models with fast and slow progression,” Mathematical Biosciences and Engineering , vol. 3, no. 4, pp. 603–614, 2006. C. Castillo-Chavez and Z. Feng, “ To treat or not to treat: the case of tuberculosis ,” Journal of Mathematical Biology , vol. 35, no. 6, pp. 629–656, 1997. S. Bowong, “ Optimal control of the transmission dynamics of tuberculosis ,” Nonlinear Dynamics , vol. 61, no. 4, pp. 729–748, 2010. S. Bowong and J. J. Tewa, “ Mathematical analysis of a tuberculosis model with differential infectivity ,” Communications in Nonlinear Science and Numerical Simulation , vol. 14, no. 11, pp. 4010–4021, 2009. H. McCallum, N. Barlow, and J. Hone, “ How should pathogen transmission be modelled? ” Trends in Ecology and Evolution , vol. 16, no. 6, pp. 295–300, 2001. S. Ruan and W. Wang, “ Dynamical behavior of an epidemic model with a nonlinear incidence rate ,” Journal of Differential Equations , vol. 188, no. 1, pp. 135–163, 2003. D. Xiao and S. Ruan, “ Global analysis of an epidemic model with nonmonotone incidence rate ,” Mathematical Biosciences , vol. 208, no. 2, pp. 419–429, 2007. Z. Yuan and L. Wang, “ Global stability of epidemiological models with group mixing and nonlinear incidence rates ,” Nonlinear Analysis: Real World Applications , vol. 11, no. 2, pp. 995–1004, 2010. Y. Tang, D. Huang, S. Ruan, and W. Zhang, “ Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate ,” SIAM Journal on Applied Mathematics , vol. 69, no. 2, pp. 621–639, 2008. M. Y. Li and J. S. Muldowney, “ Global stability for the SEIR model in epidemiology ,” Mathematical Biosciences , vol. 125, no. 2, pp. 155–164, 1995. H. R. Thieme and C. Castillo-Chavez, “On the role of variable infectivity in the dynamics of the human immunodeficiency virus epidemic,” in Mathematical and Statistical Approaches to AIDS Epidemiology , vol. 83 of Lecture Notes in Biomathematics , Springer, Berlin, Germany, 1989. J. Zhang and Z. Ma, “ Global dynamics of an SEIR epidemic model with saturating contact rate ,” Mathematical Biosciences , vol. 185, no. 1, pp. 15–32, 2003. J. A. P. Heesterbeek and J. A. J. Metz, “ The saturating contact rate in marriage- and epidemic models ,” Journal of Mathematical Biology , vol. 31, no. 5, pp. 529–539, 1993. S. Bowong and J. J. Tewa, “ Global analysis of a dynamical model for transmission of tuberculosis with a general contact rate ,” Communications in Nonlinear Science and Numerical Simulation , vol. 15, no. 11, pp. 3621–3631, 2010. P. van den Driessche and J. Watmough, “ Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission ,” Mathematical Biosciences , vol. 180, pp. 29–48, 2002. C. Castillo-Chavez and H. R. Thieme, “Asymptotically autonomous epidemic models,” in Mathematical Population Dynamics: Analysis of Heterogeneity , O. Arino, et al., Ed., vol. 1 of Theory of Epidemics , pp. 33–50, Wuetz, 1995. K. Mischaikow, H. Smith, and H. R. Thieme, “ Asymptotically autonomous semiflows: chain recurrence and Lyapunov functions ,” Transactions of the American Mathematical Society , vol. 347, no. 5, pp. 1669–1685, 1995. H. R. Thieme, “ Persistence under relaxed point-dissipativity (with application to an endemic model) ,” SIAM Journal on Mathematical Analysis , vol. 24, no. 2, pp. 407–435, 1993. J. P. LaSalle, The Stability of Dynamical Systems , SIAM, Philadelphia, Pa, USA, 1976, With an appendix: “Limiting equations and stability of nonautonomous ordinary differential equations” by Z. Artstein, Regional Conference Series in Applied Mathematics. J. P. LaSalle, “ Stability theory for ordinary differential equations ,” Journal of Differential Equations , vol. 4, pp. 57–65, 1968. A. Korobeinikov and P. K. Maini, “A Lyapunov function and global properties for SIR and SEIR epidemiological models with nonlinear incidence,” Mathematical Biosciences and Engineering , vol. 1, no. 1, pp. 57–60, 2004. // http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Abstract and Applied Analysis Hindawi Publishing Corporation

Stability of a Two-Strain Tuberculosis Model with General Contact Rate

Loading next page...
 
/lp/hindawi-publishing-corporation/stability-of-a-two-strain-tuberculosis-model-with-general-contact-rate-E0DaPv7Glt

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Hindawi Publishing Corporation
Copyright
Copyright © 2010 Hai-Feng Huo et al.
ISSN
1085-3375
eISSN
1687-0409
Publisher site
See Article on Publisher Site

Abstract

Stability of a Two-Strain Tuberculosis Model with General Contact Rate //// Hindawi Publishing Corporation Home Journals About Us About this Journal Submit a Manuscript Table of Contents Journal Menu Abstracting and Indexing Aims and Scope Annual Issues Article Processing Charges Articles in Press Author Guidelines Bibliographic Information Contact Information Editorial Board Editorial Workflow Free eTOC Alerts Reviewers Acknowledgment Subscription Information Open Special Issues Published Special Issues Special Issue Guidelines Abstract Full-Text PDF Full-Text HTML Linked References How to Cite this Article Abstract and Applied Analysis Volume 2010 (2010), Article ID 293747, 31 pages doi:10.1155/2010/293747 Research Article <h2>Stability of a Two-Strain Tuberculosis Model with General Contact Rate</h2> Hai-Feng Huo , 1 Shuai-Jun Dang , 1 and Yu-Ning Li 2 1 Institute of Applied Mathematics, Lanzhou University of Technology, Lanzhou, Gansu 730050, China 2 Department of Pediatrics, First Hospital of Lanzhou University, Lanzhou, Gansu 730000, China Received 1 November 2010; Accepted 31 December 2010 Academic Editor: D. Anderson Copyright © 2010 Hai-Feng Huo et al. This is an open access article distributed under the Creative Commons Attribution License , which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract A two-strain tuberculosis model with general contact rate which allows tuberculosis patients with the drug-sensitive Mycobacterium tuberculosis strain to be treated is presented. The model includes both drug-sensitive and drug-resistant strains. A detailed qualitative analysis about positivity, boundedness, existence, uniqueness and global stability of the equilibria of this model is carried out. Analytical results of the model show that the quantities 𝑅 1 and 𝑅 2 , which represent the basic reproduction numbers of the sensitive and resistant strains, respectively, provide the threshold conditions which determine the competitive outcomes of the two strains. Numerical simulations are also conducted to confirm and extend the analytic results. 1. Introduction Tuberculosis (TB), an infectious disease caused by Mycobacterium tuberculosis ( M. tuberculosis ), is one of the world's leading causes of infectious mortality. According to the World Health Organization(WHO), one third of the world's population is infected with M. tuberculosis , leading to between two and three millions death each year. At present, about 95% of the estimated 8 million new cases of TB occurring each year are in developing countries, where 80% occur among people between the ages of 15 and 59 years [ 1 ]. For the time being, TB is becoming a world-wide problem. War, famine, homelessness, and a lack of medical care all contribute to the increasing incidence of tuberculosis among disadvantaged persons. Since TB is easily transmissible between persons, the increase in TB in any segment of the population represents a threat to all segments of the population. Sub-Saharan Africa remains the epicenter of the epidemic, but India, China, Indonesia, Bangladesh and Pakistan together account for more than half of the cases in the world [ 2 ]. TB was assumed to be on its way out in developed countries until the number of TB cases began to increase in the 1980s. With this return, we face the paradox of a well-known bacteria, fully treatable with efficient and affordable drugs according to internationally recommended guidelines, which yet causes increasing human suffering and death. Active TB cases may be pulmonary or extrapulmonary, but pulmonary cases are more infectious and form the bulk of most cases of active TB. The usual symptoms of active TB include fatigue, high fever, night sweats, loss of appetite, and a cough, but confirmation of active TB requires a positive sputum culture. Extrapulmonary accounts for between 5% and 30% of the total cases and may affect any part of the body. M. tuberculosis droplets are released into the air by sneezing and coughing of infectious individuals. Tubercle bacillus carried by such droplets live in the air for a short period of time [ 3 ]; it is believed that occasional contacts with an infectious individual rarely lead to transmission. TB is described as a slow disease because of its long and variable latency period distribution and its short and relatively narrow infectious period distribution. Individuals who are latently infected are neither clinically ill nor capable of transmitting TB [ 4 ]. Most latently infected individuals do not become infectious (active TB). About 5%–10% of the latently infected individuals develop active TB, that is, about 90%–95% remain latently infected. Antituberculosis drugs are a two-edged sword. While they destroy pathogenic M. tuberculosis , they also select for drug-resistant bacteria against with those drugs which are then ineffective. Global surveillance has shown that drug-resistant TB is widespread and is now a threat to TB control programs in many countries. The WHO distinguishes between two types of resistance: acquired resistance—resistance among previously treated patients; and primary resistance—resistance among new cases. In all regions studied, prevalence of acquired resistance is higher than prevalence of primary resistance, but the size of this difference varies between regions [ 5 ]. Treatment of TB consists of a combination of different drugs to avoid acquisition of resistance. Despite these precautions, drug resistance continues to emerge being favoured by the long duration of treatment and improper use of the antibiotics. Drug resistant TB has higher rates of treatment failure and longer periods of infectiousness in part due to the time lapse between TB diagnosis and obtaining drug-sensitivity test results [ 6 ]. In general, tuberculosis can be treated with antibiotics. However, most worrisome is resistance to the two first-line drugs, isoniazid and rifampicin, defined as multidrug resistance(MDR)-TB, which is an emerging problem of great importance to public health, with higher mortality rates than drug-sensitive TB, particularly in immunocompromised patients. MDR-TB patients require treatment with more toxic second line drugs and remain infectious for longer than patients infected with drug-sensitive strains, incurring higher costs due to prolonged hospitalization. Mathematical models can provide useful tools to analyze the spread and control of infectious diseases [ 7 ]. Different mathematical models for TB have been formulated and studied [ 3 , 5 , 8 – 16 ]. In most of the epidemics models for TB discussed in the literature, the question of contact rate has not been a central one. Nevertheless, the mode of transmission is crucially important for two reasons. First, it determines the probable response of the disease to control, second, the objective in many models of disease in animals is to predict what will happen when a pathogen is introduced into a system in which it does not currently exist [ 17 ]. Even more important, several laboratory studies have found that the mass action incidence rate is inadequate for describing pathogen transmission and standard incidence rate is considered, for human disease, more accurate than mass action incidence rate [ 17 ]. Moreover, epidemic models with nonlinear incidence rate have recently attracted a great deal of attention from mathematical models [ 18 – 22 ]. In fact, the incidence of a disease which is the number of new cases per unit time plays an important role in the study of mathematical epidemiology. Thieme and Castillo-Chavez [ 23 ] argued that the general form of a population-size-dependent incidence should be written as 𝛽 𝐶 ( 𝑁 ) 𝑆 𝐼 / 𝑁 , where 𝑆 and 𝐼 are the numbers of susceptible and infective at time 𝑡 , respectively, 𝛽 is the probability of transmitting the infection between two individuals taking part in a contact per unit time, 𝐶 ( 𝑁 ) is the probability for an individual to take part in a contact and 𝑁 is the size of the total population. In [ 24 ], Zhang and Ma think the above incidence frequently takes two forms in most of the literature. When 𝐶 ( 𝑁 ) = 𝑁 , the corresponding incidence is the bilinear form 𝛽 𝑁 𝑆 𝐼 / 𝑁 = 𝛽 𝑆 𝐼 . When 𝐶 ( 𝑁 ) = 𝜆 , the corresponding incidence 𝜆 𝑆 𝐼 / 𝑁 is called standard form. When the total population size 𝑁 is not quite large, since the number of contacts made by an individual per unit time should increase as the total population size 𝑁 increases, the bilinear form would be suitable, but when the total population size is too large, since the number of contacts made by an infective per unit time should be limited, or should grow less rapidly as the total population size 𝑁 increases, the bilinear form is not suitable and the standard form may be more realistic [ 24 ]. Therefore, the two forms of incidence mentioned above are actually two extreme cases for the total population size 𝑁 being very small and very large, respectively. Moreover, Heesterbeek and Metz [ 25 ] derived the expression for the saturating contact rate of individuals contacts in a population that mixes randomly, that is, 𝐶 ( 𝑁 ) = 𝑏 𝑁 √ 1 + 𝑏 𝑁 + , 1 + 2 𝑏 𝑁 ( 1 . 1 ) where 𝐶 ( 𝑁 ) is nondecreasing and 𝐶 ( 𝑁 ) / 𝑁 is nonincreasing. In Bowong and Tewa [ 26 ], a SEI type of tuberculosis model with a general contact rate is proposed to study the global asymptotically stability of disease-free equilibrium and endemic equilibrium. The stability of equilibria is derived through the use of Lyapunov stability theory and LaSalle's invariant set theorem. Nevertheless, the effect of drug-resistant strain of tuberculosis is not taken into account in their article. In this paper, we study a tuberculosis model with general contact rate that includes explicitly both drug-sensitive strain and drug-resistant strain. We adapt the approach of Bowong and Tewa [ 26 ] for modeling the effective chemoprophylaxis (given to latently infected individuals) and therapeutic treatments (given to infectious). We extend their model by including the latently infected with resistant strain class and infectious with resistant strain class. The new model allows us to examine the effects of drug treatment on the prevalence of both drug-sensitive and drug-resistant strains. Mathematical properties of the model system are studied both analytically and numerically. It is shown that the system has three possible equilibrium points including an endemic equilibrium at which both strains are present. A detailed analysis of global asymptotically stability is conducted, which shows that the dynamic behaviors of the system are determined by two quantities, 𝑅 1 and 𝑅 2 , which represent the basic reproduction numbers of the sensitive and resistant strains. Without regard to disease-induced death rate, we prove that the disease-free equilibrium is globally asymptotically stable when 𝑅 1 < 1 , 𝑅 2 < 1 , the unique drug-resistant-TB-strain-only equilibrium exists and is globally asymptotically stable when 𝑅 1 < 1 < 𝑅 2 . As for the unique interior equilibrium, existence, uniqueness and global stability can be proved under certain conditions. The paper is organized as follows. In Section 2 , we formulate a two-strain tuberculosis model with general contact rate. The threshold conditions for the existence and uniqueness equilibria are derived and the global asymptotically stability of equilibria are proved in Sections 3 and 4 is devoted to numerical simulations. In Section 5 , we summarize the findings and conclusions. 2. Model Description According to the transmitted features of TB, we subdivide the population into susceptible individuals ( 𝑆 ) , those exposed to drug-sensitive TB ( 𝐸 1 ) , individuals with symptoms of TB and drug-sensitive ( 𝐼 1 ) , those exposed to drug-resistant TB ( 𝐸 2 ) and those displaying symptoms of TB and drug-resistant ( 𝐼 2 ) , the total number of population at time 𝑡 is given by 𝑁 ( 𝑡 ) = 𝑆 ( 𝑡 ) + 𝐸 1 ( 𝑡 ) + 𝐼 1 ( 𝑡 ) + 𝐸 2 ( 𝑡 ) + 𝐼 2 ( 𝑡 ) . ( 2 . 1 ) The model is represented by the transfer diagram in Figure 1 . Here, Λ is the recruitment rate, 𝜇 is nondisease related death rate. Since exposed individuals are not capable of transmitting the disease, we assume that susceptible may become infected only through contacts with active infectious individual at a rate 𝛽 𝑖 ( 𝑁 ) 𝐼 𝑖 , 𝑖 = 1 , 2 , where 𝑖 = 1 and 𝑖 = 2 represent rates of infection by drug-sensitive strain and drug-resistant strain, and the transmission coefficient 𝛽 𝑖 ( 𝑁 ) is a nonnegative 𝐶 2 function of the total population 𝑁 . Susceptible are infected with drug-sensitive strain and drug-resistant strain entering classes 𝐸 1 and 𝐸 2 , respectively. For simplicity, we only account for treatment of sensitive strain, we assume that chemoprophylaxis of individuals in 𝐸 1 reduces their reactivation at a constant rate 𝑟 1 and that the initiation of therapeutics immediately removes individuals in 𝐼 1 from active status and places them into exposed class at a constant rate 𝑟 2 . The time before individuals in 𝐸 1 who does not received effective chemoprophylaxis become infectious is assumed to satisfy an exponential distribution, with mean waiting time 1 / 𝑘 1 . Thus, individuals leave the class 𝐸 1 to the class 𝐼 1 at a constant rate 𝑘 1 ( 1 − 𝑟 1 ) . The time before individuals in 𝐸 2 become infectious is also assumed to satisfy an exponential distribution, with mean waiting time 1 / 𝑘 2 . Thus, individuals leave the class 𝐸 2 to the class 𝐼 2 at a constant rate 𝑘 2 . Individuals sick with drug-sensitive strain receive treatment at rate 𝑟 2 and a proportion 𝑝 𝑟 2 respond to treatment and move into exposed class 𝐸 1 , in the remaining proportion ( 1 − 𝑝 ) 𝑟 2 , inappropriate treatment results in the development of drug-resistant strain and the individuals move into class 𝐸 2 . Individuals in the infectious class 𝐼 1 and 𝐼 2 suffer additional disease-induced death at rate 𝑑 1 and 𝑑 2 , where 𝑑 2 > 𝑑 1 . Then the following two strains tuberculosis model is formulated 𝑆  = Λ − 𝛽 1 ( 𝑁 ) 𝑆 𝐼 1 − 𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 𝐸 − 𝜇 𝑆 ,  1 = 𝛽 1 ( 𝑁 ) 𝑆 𝐼 1 + 𝑝 𝑟 2 𝐼 1 − 𝑘 1  1 − 𝑟 1  𝐸 1 − 𝜇 𝐸 1 , 𝐼  1 = 𝑘 1  1 − 𝑟 1  𝐸 1 −  𝑟 2 + 𝜇 + 𝑑 1  𝐼 1 , 𝐸  2 = 𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 1 −  𝜇 + 𝑘 2  𝐸 2 , 𝐼  2 = 𝑘 2 𝐸 2 −  𝜇 + 𝑑 2  𝐼 2 , ( 2 . 2 ) where Λ , 𝑘 1 , 𝑘 2 , 𝑟 1 , and 𝑟 2 are assumed to be positive and 𝑝 ∈ [ 0 , 1 ] . It is natural to assume that the transmission coefficient 𝛽 1 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) satisfies the following conditions: 𝛽 1 ( 𝑁 ) > 0 , 𝛽  1  ( 𝑁 ) ≤ 0 , 𝑁 𝛽 1  ( 𝑁 )  𝛽 ≥ 0 , 2 ( 𝑁 ) > 0 , 𝛽  2  ( 𝑁 ) ≤ 0 , 𝑁 𝛽 2  ( 𝑁 )  ≥ 0 . ( 2 . 3 ) Figure 1: Transfer diagram for the dynamics of tuberculosis with general contact rate. Obviously, we can see that 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 corresponds to the standard incidence rate, that 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 corresponds to the mass action incidence rate, and that 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) corresponds to the saturating contact rate, where 𝐶 ( 𝑁 ) = 𝑏 𝑁 √ 1 + 𝑏 𝑁 + . 1 + 2 𝑏 𝑁 ( 2 . 4 ) 2.1. Positivity and Boundedness System ( 2.2 ) describes human population and therefore it is necessary to prove that all the variables 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 are nonnegative for all time. Solutions of system ( 2.2 ) with positive initial data remains positive for all time 𝑡 ≥ 0 and are bounded. Theorem 2.1. If 𝑆 ( 0 ) ≥ 0 , 𝐸 1 ( 0 ) ≥ 0 , 𝐼 1 ( 0 ) ≥ 0 , 𝐸 2 ( 0 ) ≥ 0 , 𝐼 2 ( 0 ) ≥ 0 , the solutions 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2 ( 𝑡 ) of system ( 2.2 ) are positive for 𝑡 ≥ 0 . For system ( 2.2 ), the region Γ is positively invariant and all solutions starting in Γ approach, enter, or stay in Γ , where   Γ = 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2  ( 𝑡 ) ∈ ℝ 5 + Λ ∶ 𝑁 ( 𝑡 ) ≤ 𝜇  . ( 2 . 5 ) Proof. Under the given initial conditions, it is easy to prove that the solutions of the system ( 2.2 ) are positive; if not, we assume a contradiction: that there exists a first time 𝑡 1 such that 𝑆  𝑡 1  = 0 , 𝑆   𝑡 1  < 0 , 𝐸 1 ( 𝑡 ) ≥ 0 , 𝐼 1 ( 𝑡 ) ≥ 0 , 𝐸 2 ( 𝑡 ) ≥ 0 , 𝐼 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 1 , ( 2 . 6 ) there exists a 𝑡 2 𝐸 1  𝑡 2  = 0 , 𝐸  1  𝑡 2  < 0 , 𝑆 ( 𝑡 ) ≥ 0 , 𝐼 1 ( 𝑡 ) ≥ 0 , 𝐸 2 ( 𝑡 ) ≥ 0 , 𝐼 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 2 , ( 2 . 7 ) there exists a 𝑡 3 𝐼 1  𝑡 3  = 0 , 𝐼  1  𝑡 3  < 0 , 𝑆 ( 𝑡 ) ≥ 0 , 𝐸 1 ( 𝑡 ) ≥ 0 , 𝐸 2 ( 𝑡 ) ≥ 0 , 𝐼 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 3 , ( 2 . 8 ) there exists a 𝑡 4 𝐸 2  𝑡 4  = 0 , 𝐸  2  𝑡 4  < 0 , 𝑆 ( 𝑡 ) ≥ 0 , 𝐸 1 ( 𝑡 ) ≥ 0 , 𝐼 1 ( 𝑡 ) ≥ 0 , 𝐼 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 4 , ( 2 . 9 ) or there exists a 𝑡 5 𝐼 2  𝑡 5  = 0 , 𝐼  2  𝑡 5  < 0 , 𝑆 ( 𝑡 ) ≥ 0 , 𝐸 1 ( 𝑡 ) ≥ 0 , 𝐼 1 ( 𝑡 ) ≥ 0 , 𝐸 2 ( 𝑡 ) ≥ 0 , f o r 0 < 𝑡 < 𝑡 5 . ( 2 . 1 0 ) In the first case, we have 𝑆   𝑡 1  = Λ > 0 , ( 2 . 1 1 ) which is a contradiction meaning that 𝑆 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . In the second case, we have 𝐸  1  𝑡 2  = 𝐼 1  𝑡 2 𝛽   1  𝑡 ( 𝑁 ) 𝑆 2  + 𝑝 𝑟 2  ≥ 0 , ( 2 . 1 2 ) which is a contradiction meaning that 𝐸 1 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . In the third case, we have 𝐼  1  𝑡 3  = 𝑘 1  1 − 𝑟 1  𝐸 1  𝑡 3  ≥ 0 , ( 2 . 1 3 ) which is a contradiction meaning that 𝐼 1 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . In the fourth case, we have 𝐸  2  𝑡 4  = 𝛽 2  𝑡 ( 𝑁 ) 𝑆 4  𝐼 2  𝑡 4  + ( 1 − 𝑝 ) 𝑟 2 𝐼 1  𝑡 4  ≥ 0 , ( 2 . 1 4 ) which is a contradiction meaning that 𝐸 2 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . In the fifth case, we have 𝐼  2  𝑡 5  = 𝑘 2 𝐸 2  𝑡 5  ≥ 0 , ( 2 . 1 5 ) which is a contradiction meaning that 𝐼 2 ( 𝑡 ) ≥ 0 , 𝑡 ≥ 0 . Thus, in all cases, 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2 ( 𝑡 ) remain positive for 𝑡 ≥ 0 . Let ( 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2 ( 𝑡 ) ) ∈ ℝ 5 + be any solution with nonnegative initial condition, adding all equations in system ( 2.2 ) gives 𝑁  ( 𝑡 ) = Λ − 𝜇 𝑁 − 𝑑 1 𝐼 1 − 𝑑 2 𝐼 2 , ≤ Λ − 𝜇 𝑁 . ( 2 . 1 6 ) It follows that Λ 0 ≤ 𝑁 ( 𝑡 ) ≤ 𝜇 + 𝑁 ( 0 ) 𝑒 − 𝜇 𝑡 , ( 2 . 1 7 ) where 𝑁 ( 0 ) represents initial values of the total population. Thus 0 ≤ 𝑁 ( 𝑡 ) ≤ Λ / 𝜇 , as 𝑡 → ∞ . Therefore all feasible solutions of system ( 2.2 ) enter the region   Γ = 𝑆 ( 𝑡 ) , 𝐸 1 ( 𝑡 ) , 𝐼 1 ( 𝑡 ) , 𝐸 2 ( 𝑡 ) , 𝐼 2  ( 𝑡 ) ∈ ℝ 5 + Λ ∶ 𝑁 ( 𝑡 ) ≤ 𝜇  . ( 2 . 1 8 ) Hence, Γ is positively invariant and it is sufficient to consider solutions of system ( 2.2 ) in Γ . Existence, uniqueness and continuation results for system ( 2.2 ) hold in this region. It can be shown that 𝑁 ( 𝑡 ) is bounded and all the solutions starting in Γ approach, enter or stay in Γ . 3. Model Analysis There are one disease-free equilibrium 𝑃 0 and two possible endemic equilibria for system ( 2.2 ), the endemic equilibria include boundary equilibrium 𝑃 1 (when only the drug-resistant TB-strain is present) and the interior equilibrium 𝑃 2 (when both strains exist). 3.1. Global Stability of the Disease-Free Equilibrium The disease-free-equilibrium is given as 𝑃 0 =  𝑆 0  , 0 , 0 , 0 , 0 , 𝑆 0 = Λ 𝜇 . ( 3 . 1 ) The basic reproduction number is defined as the number of secondary infections produced by a single infectious individual during the entire infectious period. In our case the reproduction number defines the number of secondary TB infections produced by a single active TB individual during the entire infectious period. Mathematically it is defined as the spectral radius of the next generation matrix [ 27 ]. To determine the stability of the disease-free steady state 𝑃 0 , we use the results in van den Driessche and Watmough [ 27 ]. Reorder the components of 𝑃 0 as 𝐸 1 = 0 , 𝐼 1 = 0 , 𝐸 2 = 0 , 𝐼 2 = 0 , 𝑆 = Λ / 𝜇 . Set ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐹 ℱ = 1 𝐹 2 𝐹 3 𝐹 4 𝐹 5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝛽 1 ( 𝑁 ) 𝑆 𝐼 1 0 𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 𝑘 𝒱 = 1  1 − 𝑟 1  𝐸 1 + 𝜇 𝐸 1 − 𝑝 𝑟 2 𝐼 1 𝑟 2 𝐼 1 +  𝜇 + 𝑑 1  𝐼 1 − 𝑘 1  1 − 𝑟 1  𝐸 1  𝜇 + 𝑘 2  𝐸 2 − ( 1 − 𝑝 ) 𝑟 2 𝐼 1  𝜇 + 𝑑 2  𝐼 2 − 𝑘 2 𝐸 2 𝜇 𝑆 + 𝛽 1 ( 𝑁 ) 𝑆 𝐼 1 + 𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . − Λ ( 3 . 2 ) The infected compartments are 𝐸 1 , 𝐼 1 , 𝐸 2 , and 𝐼 2 . Thus ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 𝐹 = 0 𝛽 1  𝑆 0  𝑆 0 0 0 0 0 0 0 0 0 0 𝛽 1  𝑆 0  𝑆 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 𝑉 = 𝜇 + 𝑘 1  1 − 𝑟 1  − 𝑝 𝑟 2 0 0 − 𝑘 1  1 − 𝑟 1  𝜇 + 𝑑 1 + 𝑟 2 0 0 0 − ( 1 − 𝑝 ) 𝑟 2 𝜇 + 𝑘 2 0 0 0 − 𝑘 2 𝜇 + 𝑑 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ . ( 3 . 3 ) The dominant eigenvalues of 𝐹 𝑉 − 1 are given by 𝑅 1 = 𝑘 1  1 − 𝑟 1  𝛽 1  𝑆 0  𝑆 0  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  , 𝑅 2 = 𝑘 2 𝛽 2  𝑆 0  𝑆 0  𝜇 + 𝑘 2   𝜇 + 𝑑 2  , ( 3 . 4 ) where 𝑅 1 and 𝑅 2 are reproduction numbers for drug-sensitive TB strain only and drug-resistant TB strain only, respectively. It implies that the spectral radius of the matrix 𝐹 𝑉 − 1 is 𝜌  𝐹 𝑉 − 1   𝑅 = m a x 1 , 𝑅 2  . ( 3 . 5 ) If 𝑅 1 < 1 , and 𝑅 2 < 1 , then 𝜌 ( 𝐹 𝑉 − 1 ) < 1 . By Theorem 3.2 in van den Driessche and Watmough [ 27 ], we know that the disease-free steady state 𝑃 0 is locally asymptotically stable, 𝑃 0 is unstable if 𝑅 1 > 1 or 𝑅 2 > 1 . The following theorems provide the global stability of the disease-free equilibrium. To conduct an analytical analysis of global asymptotical behaviors of the disease-free equilibrium point, first of all, we assume that there is no disease-induced death rate, that is, 𝑑 1 = 0 , 𝑑 2 = 0 . Consequently, the total population size 𝑁 ( 𝑡 ) satisfies the equation 𝑑 𝑁 / 𝑑 𝑡 = Λ − 𝜇 𝑁 and 𝑁 ( 𝑡 ) → Λ / 𝜇 as 𝑡 → ∞ . Using results from Castillo-Chavez and Thieme [ 28 ] and Mischaikow et al. [ 29 ], we can obtain analytical results by considering the following limiting system of ( 2.2 ) in which the total population is assumed to be constant 𝑁 = 𝑆 0 = Λ / 𝜇 : 𝐸  1 = 𝛽 1  𝑆 0  𝐼 1  𝑆 0 − 𝐸 1 − 𝐼 1 − 𝐸 2 − 𝐼 2  + 𝑝 𝑟 2 𝐼 1 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝐸   1 , 𝐼  1 = 𝑘 1  1 − 𝑟 1  𝐸 1 −  𝑟 2  𝐼 + 𝜇 1 , 𝐸  2 = 𝛽 2  𝑆 0  𝐼 2  𝑆 0 − 𝐸 1 − 𝐼 1 − 𝐸 2 − 𝐼 2  + ( 1 − 𝑝 ) 𝑟 2 𝐼 1 −  𝜇 + 𝑘 2  𝐸 2 , 𝐼  2 = 𝑘 2 𝐸 2 − 𝜇 𝐼 2 . ( 3 . 6 ) Notice that the 𝑆 equation is eliminated from ( 3.6 ) and the variable 𝑆 is replaced by 𝑆 0 − 𝐸 1 − 𝐼 1 − 𝐸 2 − 𝐼 2 . By introducing the fractions 𝑥 1 = 𝐸 1 𝑁 , 𝑥 2 = 𝐼 1 𝑁 , 𝑦 1 = 𝐸 2 𝑁 , 𝑦 2 = 𝐼 2 𝑁 , ( 3 . 7 ) we obtain the equivalent limiting system from ( 3.6 ) 𝑥  1 = 𝑆 0 𝛽 1  𝑆 0   1 − 𝑥 1 − 𝑥 2 − 𝑦 1 − 𝑦 2  𝑥 2 + 𝑝 𝑟 2 𝑥 2 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   1 , 𝑥  2 = 𝑘 1  1 − 𝑟 1  𝑥 1 −  𝜇 + 𝑟 2  𝑥 2 , 𝑦  1 = 𝑆 0 𝛽 2  𝑆 0   1 − 𝑥 1 − 𝑥 2 − 𝑦 1 − 𝑦 2  𝑦 2 + ( 1 − 𝑝 ) 𝑟 2 𝑥 2 −  𝜇 + 𝑘 2  𝑦 1 , 𝑦  2 = 𝑘 2 𝑦 1 − 𝜇 𝑦 2 . ( 3 . 8 ) Obviously 0 ≤ 𝑥 1 + 𝑥 2 + 𝑦 1 + 𝑦 2 ≤ 1 , ( 3 . 9 ) for all time 𝑡 ≥ 0 . For a bounded real-valued function 𝑓 on [ 0 , ∞ ) , we define 𝑓 ∞ = l i m i n f 𝑡 → ∞ 𝑓 ( 𝑡 ) , 𝑓 ∞ = l i m s u p 𝑡 → ∞ 𝑓 ( 𝑡 ) . ( 3 . 1 0 ) Lemma 3.1 (see [ 30 ]). Let 𝑓 ∶ [ 0 , ∞ ) → ∞ be bounded and continuously differentiable. Then there are sequences 𝑡 𝑛 , 𝑠 𝑛 → ∞ with the following properties: 𝑓  𝑡 𝑛  ⟶ 𝑓 ∞ , 𝑓   𝑡 𝑛  𝑓  𝑠 ⟶ 0 , 𝑛  ⟶ 𝑓 ∞ , 𝑓   𝑠 𝑛  ⟶ 0 , ( 3 . 1 1 ) for 𝑛 → ∞ . Theorem 3.2. In the absence of disease-induced death rate, the disease-free equilibrium 𝑃 0 of system ( 2.2 ) is globally asymptotically stable if 𝑅 1 < 1 , 𝑅 2 < 1 . Proof. By Lemma 3.1 and the 𝑥  2 , 𝑦  2 equations in ( 3.8 ) we have 𝑥 ∞ 2 ≤ 𝑘 1  1 − 𝑟 1  𝜇 + 𝑟 2 𝑥 ∞ 1 , 𝑦 ∞ 2 ≤ 𝑘 2 𝜇 𝑦 ∞ 1 . ( 3 . 1 2 ) Using the 𝑥  1 equation in ( 3.8 ) and choosing 𝑡 𝑛 → ∞ such that 𝑥 1  𝑡 𝑛  ⟶ 𝑥 ∞ 1 , 𝑥  1  𝑡 𝑛  ⟶ 0 , 𝑡 ⟶ ∞ , ( 3 . 1 3 ) we get 0 ≤ 𝑆 0 𝛽 1  𝑆 0   1 − 𝑥 1 − 𝑥 2 − 𝑦 1 − 𝑦 2  ∞ 𝑥 ∞ 2 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1 ≤ 𝑆 0 𝛽 1  𝑆 0  𝑥 ∞ 2 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1 ≤ 𝑆 0 𝛽 1  𝑆 0  𝑘 1  1 − 𝑟 1  𝜇 + 𝑟 2 𝑥 ∞ 1 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1 ≤  𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1  𝑘 1  1 − 𝑟 1  𝑆 0 𝛽 1  𝑆 0   𝜇 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1  ≤    − 1 𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1  𝑘 1  1 − 𝑟 1  𝑆 0 𝛽 1  𝑆 0   𝜇 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1   =  − 1 𝜇 + 𝑘 1  1 − 𝑟 1 𝑥   ∞ 1  𝑅 1  . − 1 ( 3 . 1 4 ) It is shown that 𝑥 ∞ 1 ≤ 0 since 𝑅 1 < 1 , as 𝑥 ∞ 1 ≥ 0 , we have that 𝑥 ∞ 1 = 0 . The inequalities in ( 3.12 ) also imply that 𝑥 ∞ 2 = 0 . Using the 𝑦  1 equation in ( 3.8 ) and choosing 𝑠 𝑛 → ∞ such that 𝑦 1  𝑠 𝑛  ⟶ 𝑦 ∞ 1 , 𝑦  1  𝑠 𝑛  ⟶ 0 , 𝑡 → ∞ , ( 3 . 1 5 ) we get 0 ≤ 𝑆 0 𝛽 2  𝑆 0   1 − 𝑥 1 − 𝑥 2 − 𝑦 1 − 𝑦 2  ∞ 𝑦 ∞ 2 −  𝜇 + 𝑘 2  𝑦 ∞ 1 ≤ 𝑆 0 𝛽 2  𝑆 0  𝑦 ∞ 2 −  𝜇 + 𝑘 2  𝑦 ∞ 1 ≤ 𝑆 0 𝛽 2  𝑆 0  𝑘 2 𝜇 𝑦 ∞ 1 −  𝜇 + 𝑘 2  𝑦 ∞ 1 ≤  𝜇 + 𝑘 2  𝑦 ∞ 1  𝑘 2 𝑆 0 𝛽 2  𝑆 0  𝜇  𝜇 + 𝑘 2   =  − 1 𝜇 + 𝑘 2  𝑦 ∞ 1  𝑅 2  . − 1 ( 3 . 1 6 ) It is shown that 𝑦 ∞ 1 ≤ 0 since 𝑅 2 < 1 , as 𝑦 ∞ 1 ≥ 0 , we have that 𝑦 ∞ 1 = 0 . The inequalities in ( 3.12 ) also imply that 𝑦 ∞ 2 = 0 . Hence l i m 𝑡 → ∞ 𝑥 1 ( 𝑡 ) = l i m 𝑡 → ∞ 𝑥 2 ( 𝑡 ) = l i m 𝑡 → ∞ 𝑦 1 ( 𝑡 ) = l i m 𝑡 → ∞ 𝑦 2 ( 𝑡 ) = 0 . ( 3 . 1 7 ) Therefore, the disease-free equilibrium 𝑃 0 is globally asymptotically stable. Next, we consider that there is disease-induced death rate, which means that 𝑑 1 > 0 , 𝑑 2 > 0 . Then, we can get the following theorem. Theorem 3.3. For 𝑑 1 > 0 , 𝑑 2 > 0 , the disease-free equilibrium 𝑃 0 of system ( 2.2 ) is globally asymptotically stable if 𝑅 1 < 1 − 𝑘 2 𝑟 2 ( 1 − 𝑝 ) / ( ( 𝜇 + 𝑑 1 + 𝑟 2 ) ( 𝜇 + 𝑘 1 ( 1 − 𝑟 1 ) ) − 𝑘 1 𝑝 𝑟 2 ( 1 − 𝑟 1 ) ) , 𝑅 2 < 1 . Proof. Let us consider the following Lyapunov function: 𝑉  𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2  = 𝑘 1  1 − 𝑟 1  𝐸 1 +  𝜇 + 𝑘 1  1 − 𝑟 1 𝐼   1 + 𝑘 2 𝐸 2 +  𝜇 + 𝑘 2  𝐼 2 . ( 3 . 1 8 ) Its time derivative along the solutions of system ( 2.2 ) satisfies ̇ 𝑉 = 𝑘 1  1 − 𝑟 1  ̇ 𝐸 1 +  𝜇 + 𝑘 1  1 − 𝑟 1 ̇ 𝐼   1 + 𝑘 2 ̇ 𝐸 2 +  𝜇 + 𝑘 2  ̇ 𝐼 2 = 𝑘 1  1 − 𝑟 1 𝛽   1 ( 𝑁 ) 𝑆 𝐼 1 + 𝑝 𝑟 2 𝐼 1 −  𝜇 + 𝑘 1  1 − 𝑟 1 𝐸   1  +  𝜇 + 𝑘 1  1 − 𝑟 1 𝑘    1  1 − 𝑟 1  𝐸 1 −  𝜇 + 𝑑 1 + 𝑟 2  𝐼 1  + 𝑘 2  𝛽 2 ( 𝑁 ) 𝑆 𝐼 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 1 −  𝜇 + 𝑘 2  𝐸 2  +  𝜇 + 𝑘 2 𝑘   2 𝐸 2 −  𝜇 + 𝑑 2  𝐼 2  , =  𝑘 1  1 − 𝑟 1  𝛽 1 ( 𝑁 ) 𝑆 + 𝑘 1  1 − 𝑟 1  𝑝 𝑟 2 −  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1 𝐼    1 +  𝑘 2 𝛽 2  ( 𝑁 ) 𝑆 − 𝜇 + 𝑑 2   𝜇 + 𝑘 2 𝐼   2 + 𝑘 2 𝑟 2 ( 1 − 𝑝 ) 𝐼 1 . ( 3 . 1 9 ) Now, using ( 2.3 ), one has 𝛽 1 ( 𝑁 ) 𝑆 ≤ 𝛽 1 ( 𝑆 0 ) 𝑆 0 , 𝛽 2 ( 𝑁 ) 𝑆 ≤ 𝛽 2 ( 𝑆 0 ) 𝑆 0 , then ̇ 𝑉 ≤   𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1  1 − 𝑟 1  𝑝 𝑟 2 𝑅   1  𝐼 − 1 1 + 𝑘 2 𝑟 2 ( 1 − 𝑝 ) 𝐼 1 +  𝜇 + 𝑑 2   𝜇 + 𝑘 2 𝑅   2  𝐼 − 1 2 , =    𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1  1 − 𝑟 1  𝑝 𝑟 2 𝑅   1  − 1 + 𝑘 2 𝑟 2  𝐼 ( 1 − 𝑝 ) 1 +  𝜇 + 𝑑 2   𝜇 + 𝑘 2 𝑅   2  𝐼 − 1 2 . ( 3 . 2 0 ) Since  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1   = 𝜇 𝜇 + 𝑑 1 + 𝑟 2  + 𝑘 1  1 − 𝑟 1   𝜇 + 𝑑 1  + 𝑘 1 𝑟 2  1 − 𝑟 1  − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1   = 𝜇 𝜇 + 𝑑 1 + 𝑟 2  + 𝑘 1  1 − 𝑟 1   𝜇 + 𝑑 1  + 𝑘 1 𝑟 2  1 − 𝑟 1  ( 1 − 𝑝 ) ≥ 0 , ( 3 . 2 1 ) and 𝑅 1 < 1 − 𝑘 2 𝑟 2 ( 1 − 𝑝 ) / ( ( 𝜇 + 𝑑 1 + 𝑟 2 ) ( 𝜇 + 𝑘 1 ( 1 − 𝑟 1 ) ) − 𝑘 1 𝑝 𝑟 2 ( 1 − 𝑟 1 ) ) ≤ 1 , we have   𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1  1 − 𝑟 1  𝑝 𝑟 2 𝑅   1  − 1 + 𝑘 2 𝑟 2 ( 1 − 𝑝 ) ≤ 0 . ( 3 . 2 2 ) Thus, ̇ 𝑉 ≤ 0 if 𝑅 1 < 1 − 𝑘 2 𝑟 2 ( 1 − 𝑝 ) / ( ( 𝜇 + 𝑑 1 + 𝑟 2 ) ( 𝜇 + 𝑘 1 ( 1 − 𝑟 1 ) ) − 𝑘 1 𝑝 𝑟 2 ( 1 − 𝑟 1 ) ) , 𝑅 2 < 1 . It is obvious that ̇ 𝑉 = 0 if and only if 𝐼 1 = 0 , 𝐼 2 = 0 . Then, the largest invariant set of system ( 2.2 ) on the set { ( 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) ∈ ℝ 5 + , ̇ 𝑉 ( 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) = 0 } is the disease-free equilibrium 𝑃 0 . Therefore, it follows from LaSalle's invariance principle [ 31 , 32 ] that 𝑃 0 is globally stable if 𝑅 1 < 1 − 𝑘 2 𝑟 2 ( 1 − 𝑝 ) / ( ( 𝜇 + 𝑑 1 + 𝑟 2 ) ( 𝜇 + 𝑘 1 ( 1 − 𝑟 1 ) ) − 𝑘 1 𝑝 𝑟 2 ( 1 − 𝑟 1 ) ) , 𝑅 2 < 1 . 3.2. Existence and Uniqueness of the Drug-Resistant-TB-Strain-Only Equilibrium Here, we present a result concerning the existence and uniqueness of the drug-resistant-TB-strain-only equilibrium for the system ( 2.2 ). To do this, we shall make use of the basic reproduction ratio 𝑅 1 and 𝑅 2 . Let 𝑃 1   𝐸 = ( 𝑆 , 0 , 0 , 2 ,  𝐼 2 ) be boundary equilibrium of system ( 2.2 ). Then, the boundary equilibrium can be obtained by setting the right-hand side of each differential equations in system ( 2.2 ) equal to zero with the exception of second and third equation, giving Λ − 𝛽 2   𝑁   𝑆  𝐼 2  𝛽 − 𝜇 𝑆 = 0 , 2   𝑁   𝑆  𝐼 2 −  𝜇 + 𝑘 2   𝐸 2 𝑘 = 0 , 2  𝐸 2 −  𝜇 + 𝑑 2   𝐼 2 = 0 , ( 3 . 2 3 ) adding the above three equations, we have  Λ − 𝜇 𝑁 − 𝑑 2  𝐼 2 = 0 . ( 3 . 2 4 ) Using ( 3.24 ), the first and second equation of ( 3.23 ), we can easily express  𝑆 ,  𝐸 2 and  𝐼 2 in terms of  𝑁 in the form:  𝑆 = Λ 𝑑 2 𝜇 𝑑 2 + 𝛽 2   𝑁  𝑁  ,  𝐸   Λ − 𝜇 2 = Λ 𝛽 2   𝑁  𝑁    Λ − 𝜇  𝜇 + 𝑘 2   𝜇 𝑑 2 + 𝛽 2   𝑁  𝑁 ,  𝐼   Λ − 𝜇   2 =  𝑁 Λ − 𝜇 𝑑 2 . ( 3 . 2 5 ) Substituting  𝑆 ,  𝐸 2 ,  𝐼 2 in the third equations of ( 3.23 ) yields   𝑁  𝐹   𝑁  Λ − 𝜇 = 0 , ( 3 . 2 6 ) where 𝐹   𝑁  = Λ 𝑘 2 𝑑 2 𝛽 2   𝑁  −  𝜇 + 𝑘 2   𝜇 + 𝑑 2   𝜇 𝑑 2 + 𝛽 2   𝑁  𝑁   Λ − 𝜇   . ( 3 . 2 7 ) Clearly,  Λ − 𝜇 𝑁 = 0 is a fixed point of ( 3.23 ), which corresponds to the disease-free equilibrium 𝑃 0 . Since  𝑁 ∈ [ 0 , 𝑆 0 ] , we get 𝐹 ( 0 ) = Λ 𝑘 2 𝑑 2 𝛽 2  ( 0 ) − 𝜇 + 𝑘 2   𝜇 + 𝑑 2   𝜇 𝑑 2 + Λ 𝛽 2  ( 0 ) = Λ 𝑘 2 𝑑 2 𝛽 2  𝜇 ( 0 ) − 2  𝑑 + 𝜇 2 + 𝑘 2  + 𝑑 2 𝑘 2   𝜇 𝑑 2 + Λ 𝛽 2  ( 0 ) = Λ 𝑘 2 𝑑 2 𝛽 2  𝜇 ( 0 ) − 2  𝑑 + 𝜇 2 + 𝑘 2    𝜇 𝑑 2 + Λ 𝛽 2  ( 0 ) − Λ 𝑑 2 𝑘 2 𝛽 2 ( 0 ) − 𝜇 𝑘 2 𝑑 2 2  𝜇 = − 2  𝑘 + 𝜇 2 + 𝑑 2    𝜇 𝑑 2 + Λ 𝛽 2  ( 0 ) − 𝜇 𝑘 2 𝑑 2 2 , 𝐹  𝑆 0  = Λ 𝑘 2 𝑑 2 𝛽 2  𝑆 0  − 𝜇 𝑑 2  𝜇 + 𝑘 2   𝜇 + 𝑑 2   𝑘 = 𝜇 2 𝑑 2 𝛽 2  𝑆 0  𝑆 0 − 𝑑 2  𝜇 + 𝑘 2   𝜇 + 𝑑 2   = 𝜇 𝑑 2  𝜇 + 𝑘 2   𝜇 + 𝑑 2   𝑘 2 𝛽 2  𝑆 0  𝑆 0  𝜇 + 𝑘 2   𝜇 + 𝑑 2   − 1 = 𝜇 𝑑 2  𝜇 + 𝑘 2   𝜇 + 𝑑 2 𝑅   2  . − 1 ( 3 . 2 8 ) It appears that 𝐹 ( 0 ) < 0 and 𝐹 ( 𝑆 0 ) > 0 when 𝑅 2 > 1 . The existence of fixed point follows from the intermediate value-theorem. Now,  𝐹 ( 𝑁 ) is monotone increasing, so that  𝐹 ( 𝑁 ) = 0 has only one positive root in the interval [ 0 , 𝑆 0 ] . Thus, we have established the following result Lemma 3.4. When 𝑅 2 > 1 , the system ( 2.2 ) has a unique drug-resistant-TB-strain-only equilibrium 𝑃 1   𝐸 = ( 𝑆 , 0 , 0 , 2 ,  𝐼 2 ) with  𝑆 ,  𝐸 2 , and  𝐼 2 all nonnegative. 3.3. Global Stability of the Drug-Resistant-TB-Strain-Only Equilibrium Theorem 3.5. If 𝑅 1 < 1 < 𝑅 2 , the drug-resistant-TB-strain-only equilibrium 𝑃 1 exists and is locally asymptotically stable. Proof. Reorder the equilibrium 𝑃 1 as  𝐸 ( 0 , 0 , 2 ,  𝐼 2 ,  𝑆 ) . Similarly as in the proof of locally stability of the disease-free steady state 𝑃 0 , we have 𝐹 1 =  0 𝛽 1   𝑁   𝑆  , 𝑉 0 0 1 =  𝜇 + 𝑘 1  1 − 𝑟 1  − 𝑝 𝑟 2 − 𝑘 1  1 − 𝑟 1  𝜇 + 𝑑 1 + 𝑟 2  , 𝐹 1 𝑉 1 − 1  𝑘 = 𝐴 1  1 − 𝑟 1  𝛽 1   𝑁   𝑆  𝜇 + 𝑘 1  1 − 𝑟 1 𝛽   1   𝑁   𝑆  , 0 0 ( 3 . 2 9 ) where 1 𝐴 =  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  . ( 3 . 3 0 ) The spectral radius of 𝐹 1 𝑉 1 − 1 is given by 𝜌  𝐹 1 𝑉 1 − 1  = 𝑘 1  1 − 𝑟 1  𝛽 1   𝑁   𝑆  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  , ≤ 𝑘 1  1 − 𝑟 1  𝛽 1   𝑁   𝑁  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  . ( 3 . 3 1 ) Since ( 𝛽 1 ( 𝑁 ) 𝑁 )  ≥ 0 , and  𝑁 ≤ 𝑆 0 , then 𝛽 1   𝑁   𝑁 ≤ 𝛽 1  𝑆 0  𝑆 0 , ( 3 . 3 2 ) Thus 𝜌  𝐹 1 𝑉 1 − 1  ≤ 𝑘 1  1 − 𝑟 1  𝛽 1  𝑆 0  𝑆 0  𝜇 + 𝑑 1 + 𝑟 2   𝜇 + 𝑘 1  1 − 𝑟 1   − 𝑘 1 𝑝 𝑟 2  1 − 𝑟 1  = 𝑅 1 . ( 3 . 3 3 ) As 𝑅 1 < 1 , we have 𝜌 ( 𝐹 1 𝑉 1 − 1 ) < 1 , which implies that the drug-resistant-TB-strain-only equilibrium 𝑃 1 is locally asymptotically stable by Theorem 3.2 in van den Driessche and Watmough [ 27 ]. The following theorem provides the global stability of drug-resistant-TB-strain-only equilibrium. Theorem 3.6. If there is no disease-induced death rate, the drug-resistant-TB-strain-only equilibrium 𝑃 1 is global asymptotically stable when 𝑅 1 < 1 < 𝑅 2 . Proof. If 𝑅 1 < 1 < 𝑅 2 , from Theorem 3.5 , 𝑃 1 is locally asymptotically stable. In the following, we only need to prove that the 𝑃 1 is a global attractor. There is no disease-induced death rate, which means 𝑑 1 = 𝑑 2 = 0 , using a similar argument as in the proof of Theorem 3.2 , we can show that if 𝑅 1 < 1 , l i m 𝑡 → ∞ 𝑥 1 ( 𝑡 ) = l i m 𝑡 → ∞ 𝑥 2 ( 𝑡 ) = 0 . ( 3 . 3 4 ) Then the equivalent limiting system of ( 3.6 ) is 𝐸  2 = 𝛽 2  𝑆 0  𝐼 2  𝑆 0 − 𝐸 2 − 𝐼 2  −  𝜇 + 𝑘 2  𝐸 2 , 𝐼  2 = 𝑘 2 𝐸 2 − 𝜇 𝐼 2 . ( 3 . 3 5 ) Let ( 𝑃 , 𝑄 ) be the vector field defined by system ( 3.35 ). Then for the Dulac function 𝐷 ( 𝐸 2 , 𝐼 2 ) = 1 / 𝐸 2 𝐼 2 , there holds 𝜕 𝐸 2 𝜕 ( 𝐷 𝑃 ) + 𝐼 2 𝛽 ( 𝐷 𝑄 ) = − 2  𝑆 0 𝑆   0 − 𝐼 2  𝐸 2 2 − 𝑘 2 𝐼 2 2 < 0 . ( 3 . 3 6 ) From Dulac's criterion, we can see equivalent limiting system ( 3.35 ) does not have a limit cycle. Therefore, the local stability of 𝑃 1 implies the global stability in Int Γ . This completes the proof of Theorem 3.6 . 3.4. Existence and Uniqueness of the Interior Equilibrium First of all, we assume that (W 1 ) l i m 𝑁 → 0 ( 1 / 𝛽 1 ( 𝑁 ) ) = + ∞ , (W 2 ) 𝛽 1 ( 𝑁 ) / 𝛽 2 ( 𝑁 ) = 𝛽 1 / 𝛽 2 , then, existence, uniqueness and global stability of interior equilibrium can be obtained. Let 𝑃 2 = ( 𝑆 ∗ , 𝐸 ∗ 1 , 𝐼 ∗ 1 , 𝐸 ∗ 2 , 𝐼 ∗ 2 ) be interior equilibrium of system ( 2.2 ). Then, the interior equilibrium (steady state with 𝐼 1 > 0 , 𝐼 2 > 0 ) can be obtained by setting the right-hand side of each of the five differential equations in system ( 2.2 ) equal to zero, giving Λ − 𝛽 1  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1 − 𝛽 2  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 2 − 𝜇 𝑆 ∗ 𝛽 = 0 , 1  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1 + 𝑝 𝑟 2 𝐼 ∗ 1 − 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 − 𝜇 𝐸 ∗ 1 𝑘 = 0 , 1  1 − 𝑟 1  𝐸 ∗ 1 −  𝑟 2 + 𝜇 + 𝑑 1  𝐼 ∗ 1 𝛽 = 0 , 2  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 −  𝜇 + 𝑘 2  𝐸 ∗ 2 𝑘 = 0 , 2 𝐸 ∗ 2 −  𝜇 + 𝑑 2  𝐼 ∗ 2 = 0 , ( 3 . 3 7 ) adding the above five equations, we have Λ − 𝜇 𝑁 ∗ − 𝑑 1 𝐼 ∗ 1 − 𝑑 2 𝐼 ∗ 2 = 0 . ( 3 . 3 8 ) Using ( 3.37 ) and ( 3.38 ), we get the following equivalent equations set Λ − 𝛽 1  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1 − 𝛽 2  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 2 − 𝜇 𝑆 ∗ 𝛽 = 0 , 1  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1 + 𝑝 𝑟 2 𝐼 ∗ 1 − 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 − 𝜇 𝐸 ∗ 1 𝑘 = 0 , 1  1 − 𝑟 1  𝐸 ∗ 1 −  𝑟 2 + 𝜇 + 𝑑 1  𝐼 ∗ 1 𝑘 = 0 , 2 𝐸 ∗ 2 −  𝜇 + 𝑑 2  𝐼 ∗ 2 = 0 , Λ − 𝜇 𝑁 ∗ − 𝑑 1 𝐼 ∗ 1 − 𝑑 2 𝐼 ∗ 2 = 0 . ( 3 . 3 9 ) From ( 3.39 ), we have 𝑆 ∗ = 𝛽  𝑆 0  𝑆 0 𝑅 1 𝛽 1 ( 𝑁 ∗ ) , 𝐼 ∗ 1 = 𝑑 2  Λ − 𝜇 𝑆 ∗  − 𝑆 ∗ 𝛽 2  𝑁 ∗   Λ − 𝜇 𝑁 ∗  𝑑 2 𝑆 ∗ 𝛽 1 ( 𝑁 ∗ ) − 𝑑 1 𝑆 ∗ 𝛽 2 ( 𝑁 ∗ ) , 𝐸 ∗ 1 = 𝜇 + 𝑑 1 + 𝑟 2 𝑘 1  1 − 𝑟 1  𝐼 ∗ 1 , 𝐼 ∗ 2 =  Λ − 𝜇 𝑁 ∗  𝑑 2 − 𝑑 1 𝑑 2 𝐼 ∗ 1 , 𝐸 ∗ 2 =  𝜇 + 𝑑 2   Λ − 𝜇 𝑁 ∗  𝑘 2 𝑑 2 − 𝑑 1  𝜇 + 𝑑 2  𝑘 2 𝑑 2 𝐼 ∗ 1 . ( 3 . 4 0 ) Substituting ( 3.40 ) into 𝛽 2 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 − ( 𝜇 + 𝑘 2 ) 𝐸 ∗ 2 = 0 , 𝑁 ∗ satisfies the following equation: 𝐺  𝑁 ∗  = 0 , ( 3 . 4 1 ) where 𝐺  𝑁 ∗  = 𝛽 1  𝑆 0  𝑆 0  Λ − 𝜇 𝑁 ∗  𝑅 1  𝑘 2 𝑑 2 𝑆 0 𝛽 1  𝑆 0  𝛽 2  𝑁 ∗  𝑅 1 𝛽 1 ( 𝑁 ∗ ) − 𝑑 2 𝜇 1 − ( 1 − 𝑝 ) 𝑘 2 𝑑 2 𝑟 2 𝛽 2  𝑁 ∗  𝛽 1 ( 𝑁 ∗ )  + 𝑑 2 Λ  𝛽 1 − 1  𝑆 0  𝛽 1 ( 𝑁 ∗ ) 𝑅 1 𝑑   1 𝜇 1 + ( 1 − 𝑝 ) 𝑘 2 𝑑 2 𝑟 2 − 𝑑 1 𝑘 2 𝑆 0 𝛽 1  𝑆 0  𝛽 2  𝑁 ∗  𝑅 1 𝛽 1 ( 𝑁 ∗ )  , 𝜇 1 =  𝜇 + 𝑘 2   𝜇 + 𝑑 2  . ( 3 . 4 2 ) Since 𝛽 1 ( 𝑁 ) / 𝛽 2 ( 𝑁 ) = 𝛽 1 / 𝛽 2 , l i m 𝑁 → 0 ( 1 / 𝛽 1 ( 𝑁 ) ) = + ∞ , if 𝑅 1 > 𝑑 1 𝑘 2 𝛽 2 𝑆 0 𝛽 1 ( 𝑆 0 ) / 𝛽 1 [ 𝑑 1 𝜇 1 + ( 1 − 𝑝 ) 𝑘 2 𝑑 2 𝑟 2 ] , then l i m 𝑁 ∗ → 0 𝐺  𝑁 ∗  = − ∞ . ( 3 . 4 3 ) It means that 𝐺 ( 𝑁 ∗ ) = 0 has only one positive root in the interval [ 0 , 𝑆 0 ] if 𝐺 ( 𝑁 ∗ ) is monotone increasing and 𝐺 ( 𝑆 0 ) ≥ 0 , which is equivalent to 𝑅 1 > 𝐻 1 + 𝐻 2 𝐻 3 + 𝐻 4 , 𝑅 2 <  1 + ( 1 − 𝑝 ) 𝑑 2 𝑘 2 𝑟 2 𝑑 1 𝜇 1  𝑅 1 , ( 3 . 4 4 ) where 𝐻 1 = 𝑘 2 𝑑 2 𝛽 2 𝑆 0 𝛽 1  𝑆 0  𝛽 2 1  𝑁 ∗  , 𝐻 2 = 𝑑 1 𝑑 2 𝑘 2 𝛽 2 𝑆 0 𝛽 1  𝑆 0  𝛽  1  𝑁 ∗  , 𝐻 3 = 𝛽 2 1  𝑁 ∗ 𝑑   2 𝛽 1 𝜇 1 + 𝑑 2 𝑘 2 𝛽 2 𝑟 2  , 𝐻 ( 1 − 𝑝 ) 4 = 𝑑 2 𝛽  1  𝑁 ∗ 𝑑   1 𝛽 1 𝜇 1 + 𝑑 2 𝑘 2 𝛽 1 𝑟 2  . ( 1 − 𝑝 ) ( 3 . 4 5 ) Let  𝐻 𝐻 = m a x 1 , 1 + 𝐻 2 𝐻 3 + 𝐻 4 , 𝑑 1 𝑘 2 𝛽 2 𝑆 0 𝛽 1  𝑆 0  𝛽 1  𝑑 1 𝜇 1 + 𝑘 2 𝑑 2 𝑟 2   , 𝑓  𝑅 ( 1 − 𝑝 ) 1  =  1 + ( 1 − 𝑝 ) 𝑑 2 𝑘 2 𝑟 2 𝑑 1 𝜇 1  𝑅 1 . ( 3 . 4 6 ) When 𝑅 1 > 𝐻 , 𝑅 2  𝑅 < 𝑓 1  , ( 3 . 4 7 ) it shows that 𝐺 ( 𝑆 0 ) ≥ 0 , and 𝐺 ( 𝑁 ∗ ) is monotone increasing, by the intermediate value-theorem, 𝐺 ( 𝑁 ∗ ) = 0 has only one positive root in the interval [ 0 , 𝑆 0 ] , which signify that 𝑃 2 is unique. Obviously, 𝑆 ∗ > 0 . In order to make certain that 𝐸 ∗ 1 > 0 , 𝐼 ∗ 1 > 0 , 𝐸 ∗ 2 > 0 , 𝐼 ∗ 2 > 0 , we obtain that 𝐸 ∗ 1 > 0 , 𝐼 ∗ 1 > 0 , 𝐸 ∗ 2 > 0 , 𝐼 ∗ 2 > 0 if and only if 0 < 𝐼 ∗ 1 < Λ − 𝜇 𝑁 ∗ 𝑑 1 . ( 3 . 4 8 ) Due to 𝑑 2 > 𝑑 1 > 0 , introducing ( 3.40 ) into ( 3.48 ), we acquire 𝐵 1 < 𝑅 1 < 𝐵 2 , ( 3 . 4 9 ) where 𝐵 1 = 𝑆 0 𝛽 1  𝑆 0 𝛽   2  𝑁 ∗   Λ − 𝜇 𝑁 ∗  + 𝜇 𝑑 2  Λ 𝑑 2 𝛽 1 ( 𝑁 ∗ ) , 𝐵 2 = 𝑆 0 𝛽 1  𝑆 0 𝛽   1  𝑁 ∗   Λ − 𝜇 𝑁 ∗  + 𝜇 𝑑 1  Λ 𝑑 1 𝛽 1 ( 𝑁 ∗ ) , ( 3 . 5 0 ) and 𝐵 1 < 𝐵 2 . Now, we can get the following conclusion. Lemma 3.7. When conditions (W 1 ), (W 2 ), ( 3.47 ) and ( 3.49 ) are satisfied, the system ( 2.2 ) has a unique interior equilibrium 𝑃 2 = ( 𝑆 ∗ , 𝐸 ∗ 1 , 𝐼 ∗ 1 , 𝐸 ∗ 2 , 𝐼 ∗ 2 ) with 𝑆 ∗ , 𝐸 ∗ 1 , 𝐼 ∗ 1 , 𝐸 ∗ 2 , and 𝐼 ∗ 2 all nonnegative. Remark 3.8. For 𝛽 1 ( 𝑁 ) = 𝛽 1 and 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , condition (W 1 ) is not satisfied. However, at least when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , condition (W 1 ) is satisfied, we also can find our conclusions are still correct from Figures 7 and 10 . 3.5. Global Stability of the Interior Equilibrium Herein, we study the global stability of the interior equilibrium 𝑃 2 of system ( 2.2 ). For simplicity, we only consider the special case when the transmission coefficients of drug-sensitive TB-strain and drug-resistant TB-strain are equal, that is, 𝛽 1 ( 𝑁 ) = 𝛽 2 ( 𝑁 ) = 𝛽 ( 𝑁 ) , then, we get the following result. Theorem 3.9. If conditions (W 1 ), (W 2 ), ( 3.47 ) and ( 3.49 ) are satisfied, interior equilibrium 𝑃 2 of system ( 2.2 ) is is globally asymptotically stable when 𝑆 𝑆 ∗ = 𝐸 1 𝐸 ∗ 1 = 𝐸 2 𝐸 ∗ 2 𝐼 ≥ 1 , 2 𝐼 ∗ 2 ≥ 𝐼 1 𝐼 ∗ 1 𝑆 o r 𝑆 ∗ = 𝐸 1 𝐸 ∗ 1 = 𝐸 2 𝐸 ∗ 2 𝐼 ≤ 1 , 2 𝐼 ∗ 2 ≤ 𝐼 1 𝐼 ∗ 1 . ( 3 . 5 1 ) Proof. Following [ 33 ], we consider the Lyapunov function: 𝑈  𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2  =  𝑆 − 𝑆 ∗  +  𝐸 l n 𝑆 1 − 𝐸 ∗ 1 l n 𝐸 ∗ 1   𝐼 + 𝐶 1 − 𝐼 ∗ 1 l n 𝐼 ∗ 1  +  𝐸 2 − 𝐸 ∗ 2 l n 𝐸 ∗ 2   𝐼 + 𝐷 2 − 𝐼 ∗ 2 l n 𝐼 ∗ 2  , ( 3 . 5 2 ) where 𝐶 and 𝐷 are positive constants to be determined later. Differentiating this function with respect to time along the solutions of system ( 2.2 ) yields ̇  𝑆 𝑈 = 1 − ∗ 𝑆  ̇  𝐸 𝑆 + 1 − ∗ 1 𝐸 1  ̇ 𝐸 1  𝐼 + 𝐶 1 − ∗ 1 𝐼 1  ̇ 𝐼 1 +  𝐸 1 − ∗ 2 𝐸 2  ̇ 𝐸 2  𝐼 + 𝐷 1 − ∗ 2 𝐼 2  ̇ 𝐼 2 =  𝑆 1 − ∗ 𝑆   Λ − 𝛽 ( 𝑁 ) 𝑆 𝐼 1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 2  +  𝐸 − 𝜇 𝑆 1 − ∗ 1 𝐸 1   𝛽 ( 𝑁 ) 𝑆 𝐼 1 + 𝑝 𝑟 2 𝐼 1 − 𝑘 1  1 − 𝑟 1  𝐸 1 − 𝜇 𝐸 1   𝐼 + 𝐶 1 − ∗ 1 𝐼 1   𝑘 1  1 − 𝑟 1  𝐸 1 −  𝑟 2 + 𝜇 + 𝑑 1  𝐼 1  +  𝐸 1 − ∗ 2 𝐸 2   𝛽 ( 𝑁 ) 𝑆 𝐼 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 1 −  𝜇 + 𝑘 2  𝐸 2   𝐼 + 𝐷 1 − ∗ 2 𝐼 2   𝑘 2 𝐸 2 −  𝜇 + 𝑑 2  𝐼 2  . ( 3 . 5 3 ) Considering ( 3.37 ), we can deduce that  𝑁 Λ = 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 + 𝜇 𝑆 ∗ , 𝜇 + 𝑘 1  1 − 𝑟 1   𝑁 = 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1 𝐸 ∗ 1 + 𝑝 𝑟 2 𝐼 ∗ 1 𝐸 ∗ 1 , 𝜇 + 𝑑 1 + 𝑟 2 = 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 𝐼 ∗ 1 , 𝜇 + 𝑘 2  𝑁 = 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 𝐸 ∗ 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 𝐸 ∗ 2 , 𝜇 + 𝑑 2 = 𝑘 2 𝐸 ∗ 2 𝐼 ∗ 2 . ( 3 . 5 4 ) Then, ( 3.53 ) becomes ̇  𝑆 𝑈 = 1 − ∗ 𝑆   𝛽  𝑁 ∗  𝑆 ∗ 𝐼 ∗ 1  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 + 𝜇 𝑆 ∗ − 𝛽 ( 𝑁 ) 𝑆 𝐼 1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 2  +  𝐸 − 𝜇 𝑆 1 − ∗ 1 𝐸 1   𝛽 ( 𝑁 ) 𝑆 𝐼 1 + 𝑝 𝑟 2 𝐼 1  𝑁 − 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1 𝐸 1 𝐸 ∗ 1 − 𝑝 𝑟 2 𝐼 ∗ 1 𝐸 1 𝐸 ∗ 1   𝐼 + 𝐶 1 − ∗ 1 𝐼 1 𝑘   1  1 − 𝑟 1  𝐸 1 − 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 𝐼 1 𝐼 ∗ 1  +  𝐸 1 − ∗ 2 𝐸 2   𝛽 ( 𝑁 ) 𝑆 𝐼 2 + ( 1 − 𝑝 ) 𝑟 2 𝐼 1  𝑁 − 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 𝐸 2 𝐸 ∗ 2 − ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 𝐸 2 𝐸 ∗ 2   𝐼 + 𝐷 1 − ∗ 2 𝐼 2 𝑘   2 𝐸 2 − 𝑘 2 𝐸 ∗ 2 𝐼 2 𝐼 ∗ 2  𝜇  = − 𝑆 − 𝑆 ∗  2 𝑆  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1  1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 1 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 1 𝑆   1 − ∗ 𝑆   𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2  1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 2 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 𝑆   1 − ∗ 𝑆  +  𝐸 1 − ∗ 1 𝐸 1 𝛽  𝑁   ∗  𝑆 ∗ 𝐼 ∗ 1  𝛽 ( 𝑁 ) 𝑆 𝐼 1 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 1 − 𝐸 1 𝐸 ∗ 1  + 𝑝 𝑟 2 𝐼 ∗ 1  𝐼 1 𝐼 ∗ 1 − 𝐸 1 𝐸 ∗ 1  𝐼   + 𝐶 1 − ∗ 1 𝐼 1 𝑘   1  1 − 𝑟 1  𝐸 ∗ 1  𝐸 1 𝐸 ∗ 1 − 𝐼 1 𝐼 ∗ 1 +  𝐸   1 − ∗ 2 𝐸 2 𝛽  𝑁   ∗  𝑆 ∗ 𝐼 ∗ 2  𝛽 ( 𝑁 ) 𝑆 𝐼 2 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 − 𝐸 2 𝐸 ∗ 2  + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  𝐼 1 𝐼 ∗ 1 − 𝐸 2 𝐸 ∗ 2  𝐼   + 𝐷 1 − ∗ 2 𝐼 2 𝑘   2 𝐸 ∗ 2  𝐸 2 𝐸 ∗ 2 − 𝐼 2 𝐼 ∗ 2 .   ( 3 . 5 5 ) Now, using ( 3.37 ), we have 𝑘 1  1 − 𝑟 1  𝐸 ∗ 1 =  𝑟 2 + 𝜇 + 𝑑 1  𝐼 ∗ 1 , 𝑘 2 𝐸 ∗ 2 =  𝜇 + 𝑑 2  𝐼 ∗ 2 . ( 3 . 5 6 ) Then, ( 3.55 ) can be rewritten as follows: ̇ 𝜇  𝑈 = − 𝑆 − 𝑆 ∗  2 𝑆  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1   1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 1 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 1 𝑆   1 − ∗ 𝑆  +  𝐸 1 − ∗ 1 𝐸 1   𝛽 ( 𝑁 ) 𝑆 𝐼 1 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 1 − 𝐸 1 𝐸 ∗ 1   + 𝑝 𝑟 2 𝐼 ∗ 1   𝐸 1 − ∗ 1 𝐸 1 𝐼   1 𝐼 ∗ 1 − 𝐸 1 𝐸 ∗ 1  + 𝐶  𝜇 + 𝑟 2 + 𝑑 1  𝑝 𝑟 2  𝐼 1 − ∗ 1 𝐼 1 𝐸   1 𝐸 ∗ 1 − 𝐼 1 𝐼 ∗ 1    𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2   1 − 𝛽 ( 𝑁 ) 𝑆 𝐼 2 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 𝑆   1 − ∗ 𝑆  +  𝛽 ( 𝑁 ) 𝑆 𝐼 2 𝛽 ( 𝑁 ∗ ) 𝑆 ∗ 𝐼 ∗ 2 − 𝐸 2 𝐸 ∗ 2 𝐸   1 − ∗ 2 𝐸 2   + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1   𝐼 1 𝐼 ∗ 1 − 𝐸 2 𝐸 ∗ 2 𝐸   1 − ∗ 2 𝐸 2  + 𝐷  𝜇 + 𝑑 2  𝐼 ∗ 2 ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  𝐼 1 − ∗ 2 𝐼 2 𝐸   2 𝐸 ∗ 2 − 𝐼 2 𝐼 ∗ 2   . ( 3 . 5 7 ) Now, let ( 𝑥 , 𝑤 1 , 𝑧 1 , 𝑤 2 , 𝑧 2 ) = ( 𝑆 / 𝑆 ∗ , 𝐸 1 / 𝐸 ∗ 1 , 𝐼 1 / 𝐼 ∗ 1 , 𝐸 2 / 𝐸 ∗ 2 , 𝐼 2 / 𝐼 ∗ 2 ) and 𝑔 ( 𝑁 ) = 𝛽 ( 𝑁 ) / 𝛽 ( 𝑁 ∗ ) , then we get ̇ 𝜇  𝑈 = − 𝑆 − 𝑆 ∗  2 𝑆  𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1   1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 1   1 1 − 𝑥  +  1 1 − 𝑤 1   𝑔 ( 𝑁 ) 𝑥 𝑧 1 − 𝑤 1   + 𝑝 𝑟 2 𝐼 ∗ 1   1 1 − 𝑤 1   𝑧 1 − 𝑤 1  + 𝐶  𝜇 + 𝑟 2 + 𝑑 1  𝑝 𝑟 2  1 1 − 𝑧 1   𝑤 1 − 𝑧 1    𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2   1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 2   1 1 − 𝑥  +  𝑔 ( 𝑁 ) 𝑥 𝑧 2 − 𝑤 2   1 1 − 𝑤 2   + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1   𝑧 1 − 𝑤 2   1 1 − 𝑤 2  + 𝐷  𝜇 + 𝑑 2  𝐼 ∗ 2 ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  1 1 − 𝑧 2   𝑤 2 − 𝑧 2   , 𝜇  = − 𝑆 − 𝑆 ∗  2 𝑆  + 𝑓 𝑥 , 𝑤 1 , 𝑧 1 , 𝑤 2 , 𝑧 2  , ( 3 . 5 8 ) where 𝑓  𝑥 , 𝑤 1 , 𝑧 1 , 𝑤 2 , 𝑧 2   𝑁 = 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 1 𝑓 1  𝑥 , 𝑤 1 , 𝑧 1  + 𝑝 𝑟 2 𝐼 ∗ 1 𝑓 2  𝑤 1 , 𝑧 1   𝑁 + 𝛽 ∗  𝑆 ∗ 𝐼 ∗ 2 𝑓 3  𝑥 , 𝑤 2 , 𝑧 2  + ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1 𝑓 4  𝑧 1 , 𝑤 2 , 𝑧 2  , 𝑓 1  𝑥 , 𝑤 1 , 𝑧 1  =  1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 1   1 1 − 𝑥  +  1 1 − 𝑤 1   𝑔 ( 𝑁 ) 𝑥 𝑧 1 − 𝑤 1  , 𝑓 2  𝑤 1 , 𝑧 1  =  1 1 − 𝑤 1   𝑧 1 − 𝑤 1  + 𝐶  𝜇 + 𝑟 2 + 𝑑 1  𝑝 𝑟 2  1 1 − 𝑧 1   𝑤 1 − 𝑧 1  , 𝑓 3  𝑥 , 𝑤 2 , 𝑧 2  =  1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 2   1 1 − 𝑥  +  𝑔 ( 𝑁 ) 𝑥 𝑧 2 − 𝑤 2   1 1 − 𝑤 2  , 𝑓 4  𝑧 1 , 𝑤 2 , 𝑧 2  =  𝑧 1 − 𝑤 2   1 1 − 𝑤 2  + 𝐷  𝜇 + 𝑑 2  𝐼 ∗ 2 ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  1 1 − 𝑧 2   𝑤 2 − 𝑧 2  . ( 3 . 5 9 ) We can choose 𝐶 = 𝑝 𝑟 2  𝜇 + 𝑟 2 + 𝑑 1  , 𝐷 = ( 1 − 𝑝 ) 𝑟 2 𝐼 ∗ 1  𝜇 + 𝑑 2  𝐼 ∗ 2 . ( 3 . 6 0 ) For 𝑥 = 𝑤 1 = 𝑤 2 , we have 𝑓 1  𝑥 , 𝑤 1 , 𝑧 1  =  1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 1   1 1 − 𝑥  +  1 1 − 𝑤 1   𝑔 ( 𝑁 ) 𝑥 𝑧 1 − 𝑤 1  = 𝑔 ( 𝑁 ) 𝑥 𝑧 1  1 𝑥 1 − 1 + 1 − 𝑤 1  1 + 1 − 𝑥 − 𝑤 1  1 1 − 𝑤 1  1 = 1 − 𝑥 − 𝑤 1 1 + 1 = 2 − 𝑥 − 𝑥 , ( 3 . 6 1 ) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑥 = 1 𝑓 2  𝑤 1 , 𝑧 1  =  1 1 − 𝑤 1   𝑧 1 − 𝑤 1  +  1 1 − 𝑧 1   𝑤 1 − 𝑧 1  =  𝑧 1 − 𝑤 1   1 1 − 𝑤 1 1 − 1 + 𝑧 1  =  𝑧 1 − 𝑤 1   1 𝑧 1 − 1 𝑤 1  𝑤 = 2 − 1 𝑧 1 − 𝑧 1 𝑤 1 , ( 3 . 6 2 ) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑤 1 = 𝑧 1 𝑓 3  𝑥 , 𝑤 2 , 𝑧 2  =  1 − 𝑔 ( 𝑁 ) 𝑥 𝑧 2   1 1 − 𝑥  +  𝑔 ( 𝑁 ) 𝑥 𝑧 2 − 𝑤 2   1 1 − 𝑤 2  = 𝑔 ( 𝑁 ) 𝑥 𝑧 2  1 1 − 𝑤 2 + 1 𝑥  1 − 1 + 1 − 𝑥 − 𝑤 2 1 + 1 = 1 − 𝑥 − 𝑤 2 1 + 1 = 2 − 𝑥 − 𝑥 , ( 3 . 6 3 ) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑥 = 1 . For 𝑤 2 ≥ 1 , 𝑧 2 ≥ 𝑧 1 or 𝑤 2 ≤ 1 , 𝑧 2 ≤ 𝑧 1 , we have 𝑓 4  𝑧 1 , 𝑤 2 , 𝑧 2  =  𝑧 1 − 𝑤 2   1 1 − 𝑤 2  +  1 1 − 𝑧 2   𝑤 2 − 𝑧 2  ≤  𝑧 2 − 𝑤 2   1 1 − 𝑤 2 1 − 1 + 𝑧 2  =  𝑧 2 − 𝑤 2   1 𝑧 2 − 1 𝑤 2  𝑤 = 2 − 2 𝑧 2 − 𝑧 2 𝑤 2 , ( 3 . 6 4 ) which is less than or equal to zero by the arithmetic-geometric mean inequality, with equality if and only if 𝑤 2 = 𝑧 2 . Thus, ̇ 𝑈 ( 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) is less or equal to zero with equality only if 𝑥 = 1 , 𝑤 1 = 𝑧 1 , 𝑤 2 = 𝑧 2 . Since 𝑥 = 𝑤 1 = 𝑤 2 , the largest invariant set of system ( 2.2 ) on the set { ( 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) ∈ ℝ 5 + , ̇ 𝑈 ( 𝑆 , 𝐸 1 , 𝐼 1 , 𝐸 2 , 𝐼 2 ) = 0 } is the endemic equilibrium point 𝑃 2 . By LaSalle's invariance principle [ 31 , 32 ], we can conclude that the interior equilibrium 𝑃 2 is globally asymptotically stable when condition ( 3.51 ) is satisfied. Remark 3.10. It is possible for condition ( 3.51 ) to fail, in which case the global stability of the interior equilibrium of system ( 2.2 ) has not been established. Figures 7 , 10 and 13 , however, seem to support the idea that the interior equilibrium of system ( 2.2 ) is still global asymptotically stable even in this case. 4. Numerical Simulations To illustrate the theoretical results contained in this paper, we give some simulations using the parameter value in Table 1 . Numerical results are displayed in the following figures. Table 1: Description and estimation of parameters. Taking no account of disease-induced death rate, we assume that 𝑑 1 = 𝑑 2 = 0 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 0 . 0 0 0 1 (so that 𝑅 1 < 1 , 𝑅 2 < 1 ). For 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 and 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑏 𝑁 / ( 1 + 𝑏 𝑁 + 1 + 2 𝑏 𝑁 ) , 𝑏 = 1 , Figures 2 , 3 , and 4 show the trajectories plots and their plane figures, respectively, which are in agreement with Theorem 3.2 . Figure 2: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 0 . 0 0 0 1 , 𝑑 1 = 𝑑 2 = 0 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 3: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 0 . 0 0 0 1 , 𝑑 1 = 𝑑 2 = 0 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 4: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑁 / ( 1 + 𝑁 + 1 + 2 𝑁 ) , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 0 . 0 0 0 1 , 𝑑 1 = 𝑑 2 = 0 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. When there is disease-induced death rate, which means that 𝑑 1 > 0 , 𝑑 2 > 0 , we take 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 . Firstly, we consider system ( 2.2 ) with 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 and choose 𝛽 1 = 𝛽 2 = 0 . 0 0 1 , so that 𝑅 1 < 1 , 𝑅 2 < 1 , Figure 5 presents the trajectories plot and plane figure. From this figure, we can see that the trajectories of system ( 2.2 ) converge to the disease-free equilibrium. We also choose 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 , so that 𝑅 1 < 1 , 𝑅 2 > 1 , Figure 6 presents the trajectories plot and plane figure. It shows that the drug-resistant-TB-strain-only equilibrium is global asymptotically stable. Then, we choose 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 , so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) , the trajectories plot and its plane figure are depicted in Figure 7 , we can observe that the trajectories converge to the unique interior equilibrium. This implies that both the sensitive and resistant strains coexist and prevail in the host population as shown in Theorem 3.9 . Figure 5: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 = 𝛽 2 = 0 . 0 0 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 6: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 > 1 . It shows that the boundary equilibrium 𝑃 1 is stable. Figure 7: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 , 𝛽 2 ( 𝑁 ) = 𝛽 2 , 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) . It shows that the interior equilibrium 𝑃 2 is stable. Next, we need to think the case when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 . Figures 8 , 9 and 10 show the trajectories plot and its plane figure of system ( 2.2 ) when 𝛽 1 = 𝛽 2 = 0 . 0 0 1 (so that 𝑅 1 < 1 , 𝑅 2 < 1 ), 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 (so that 𝑅 1 < 1 , 𝑅 2 > 1 ), and 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 (so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) ), respectively. As the previous cases, we can find that the trajectories of system ( 2.2 ) converge to the disease-free equilibrium when 𝑅 1 < 1 , 𝑅 2 < 1 (see Figure 8 ), the drug-resistant-TB-strain-only equilibrium is global asymptotically stable when 𝑅 1 < 1 , 𝑅 2 > 1 (see Figure 9 ), while the trajectories of system ( 2.2 ) converge to the unique interior equilibrium when 𝑅 1 > 1 , 𝑅 2 < 𝑓 ( 𝑅 1 ) (see Figure 10 ). Figure 8: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 , 𝛽 1 = 𝛽 2 = 0 . 0 0 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 9: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 > 1 . It shows that the boundary equilibrium 𝑃 1 is stable. Figure 10: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 / 𝑁 , 𝛽 2 ( 𝑁 ) = 𝛽 2 / 𝑁 , 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) . It shows that the interior equilibrium 𝑃 2 is stable. Eventually, we assume that 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑏 𝑁 / ( 1 + 𝑏 𝑁 + 1 + 2 𝑏 𝑁 ) in system ( 2.2 ), of which 𝑏 = 1 . Numerical results are depicted in Figures 11 , 12 , and 13 . Figure 11 presents the trajectories and its plane figure when 𝛽 1 = 𝛽 2 = 0 . 0 0 1 (so that 𝑅 1 < 1 , 𝑅 2 < 1 ). From this figure, it clearly appears that the disease disappears in the host population. Figure 12 presents the trajectories and its plane figure when 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 (so that 𝑅 1 < 1 , 𝑅 2 > 1 ), we can see that the drug-resistant-TB-strain-only equilibrium is global asymptotically stable. Figure 13 gives a description of the trajectories and its plane figure when 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 (so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) ). By observation, we conclude that both the sensitive and resistant strains will coexist and prevail in the host population. Figure 11: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑁 / ( 1 + 𝑁 + 1 + 2 𝑁 ) , 𝛽 1 = 𝛽 2 = 0 . 0 0 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 < 1 . It shows that the disease-free equilibrium 𝑃 0 is stable. Figure 12: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑁 / ( 1 + 𝑁 + 1 + 2 𝑁 ) , 𝛽 1 = 0 . 0 0 1 , 𝛽 2 = 1 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 < 1 , 𝑅 2 > 1 . It shows that the boundary equilibrium 𝑃 1 is stable. Figure 13: Trajectories of system ( 2.2 ) when 𝛽 1 ( 𝑁 ) = 𝛽 1 𝐶 ( 𝑁 ) , 𝛽 2 ( 𝑁 ) = 𝛽 2 𝐶 ( 𝑁 ) , √ 𝐶 ( 𝑁 ) = 𝑁 / ( 1 + 𝑁 + 1 + 2 𝑁 ) , 𝛽 1 = 1 , 𝛽 2 = 0 . 2 5 , 𝑑 1 = 0 . 0 2 , 𝑑 2 = 0 . 0 7 , so that 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) . It shows that the interior equilibrium 𝑃 2 is stable. 5. Discussion In this paper, we studied a two-strain tuberculosis model that includes both drug-sensitive and drug-resistant strains. The novelty of our model is that we incorporate general contact rate 𝛽 𝑖 ( 𝑁 ) , 𝑖 = 1 , 2 into a two-strain tuberculosis model. By using the Lyapunov stability theory, LaSalle's invariant set theorem and Dulac's criterion, the global stability of equilibria of the proposed model is proved. The basic reproduction number 𝑅 1 , 𝑅 2 provide the threshold conditions which determine the competitive outcomes of the two strains. When 𝑅 1 < 1 , 𝑅 2 < 1 , the disease-free equilibrium 𝑃 0 is globally stable and the disease will die out. When 𝑅 1 < 1 , 𝑅 2 > 1 , a unique boundary equilibrium 𝑃 1 exists and is globally stable. While when 𝑅 1 > 𝐻 , 𝑅 2 < 𝑓 ( 𝑅 1 ) , there exists a unique interior equilibrium 𝑃 2 which is globally asymptotically stable under certain conditions. A fairly good agreement is obtained between the analytical and numerical results. With the appearance of numerous second-line antituberculous drugs, such as capreomycin, amikacin, and kanamycin, the treatment of drug-resistant-TB has become possible. We will leave this for future study. Acknowledgments This work was supported by the NNSF of China (10961018), the Key Project of Chinese Ministry of Education (209131), the Project Sponsored by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry, the NSF of Gansu Province of China (3ZS042-B25-013), the NSF of Bureau of Education of Gansu Province of China (0803-01) and the Development Program for Outstanding Young Teachers in Lanzhou University of Technology (Q200703 ) and the Doctor's Foundation of Lanzhou University of Technology. <h4>References</h4> C. Dye, S. Scheele, P. Dolin, V. Pathania, and M. C. Raviglione, “ Global burden of tuberculosis: estimated incidence, prevalence, and mortality by country ,” Journal of the American Medical Association , vol. 282, no. 7, pp. 677–686, 1999. T. R. Frieden, T. R. Sterling, S. S. Munsiff, C. J. Watt, and C. Dye, “ Tuberculosis ,” Lancet , vol. 362, no. 9387, pp. 887–899, 2003. B. Song, C. Castillo-Chavez, and J. P. Aparicio, “ Tuberculosis models with fast and slow dynamics: the role of close and casual contacts ,” Mathematical Biosciences , vol. 180, pp. 187–205, 2002. B. Miller, “Preventive therapy for tuberculosis,” Medical Clinics of North America , vol. 77, no. 6, pp. 1263–1275, 1993. P. Rodrigues, M. G. M. Gomes, and C. Rebelo, “ Drug resistance in tuberculosis-a reinfection model ,” Theoretical Population Biology , vol. 71, no. 2, pp. 196–212, 2007. M. A. Espinal, et al., “Standard short-course chemotherapy for drug-resistant tuberculosis: treatment outcomes in 6 countries,” Journal of the American Medical Association , vol. 283, no. 19, pp. 2537–2545, 2000. H. W. Hethcote, “ The mathematics of infectious diseases ,” SIAM Review , vol. 42, no. 4, pp. 599–653, 2000. F. Brauer and C. Castillo-Chávez, Mathematical models in population biology and epidemiology , vol. 40 of Texts in Applied Mathematics , Springer, New York, NY, USA, 2001. B. M. Murphy, B. H. Singer, S. Anderson, and D. Kirschner, “ Comparing epidemic tuberculosis in demographically distinct heterogeneous populations ,” Mathematical Biosciences , vol. 180, pp. 161–185, 2002. C. Castillo-Chavez and B. Song, “Dynamical models of tuberculosis and their applications,” Mathematical Biosciences and Engineering , vol. 1, no. 2, pp. 361–404, 2004. B. M. Murphy, B. H. Singer, and D. Kirschner, “ On treatment of tuberculosis in heterogeneous populations ,” Journal of Theoretical Biology , vol. 223, no. 4, pp. 391–404, 2003. C. Castillo-Chavez and Z. Feng, “ Global stability of an age-structure model for TB and its applications to optimal vaccination strategies ,” Mathematical Biosciences , vol. 151, no. 2, pp. 135–154, 1998. C. C. McCluskey, “Lyapunov functions for tuberculosis models with fast and slow progression,” Mathematical Biosciences and Engineering , vol. 3, no. 4, pp. 603–614, 2006. C. Castillo-Chavez and Z. Feng, “ To treat or not to treat: the case of tuberculosis ,” Journal of Mathematical Biology , vol. 35, no. 6, pp. 629–656, 1997. S. Bowong, “ Optimal control of the transmission dynamics of tuberculosis ,” Nonlinear Dynamics , vol. 61, no. 4, pp. 729–748, 2010. S. Bowong and J. J. Tewa, “ Mathematical analysis of a tuberculosis model with differential infectivity ,” Communications in Nonlinear Science and Numerical Simulation , vol. 14, no. 11, pp. 4010–4021, 2009. H. McCallum, N. Barlow, and J. Hone, “ How should pathogen transmission be modelled? ” Trends in Ecology and Evolution , vol. 16, no. 6, pp. 295–300, 2001. S. Ruan and W. Wang, “ Dynamical behavior of an epidemic model with a nonlinear incidence rate ,” Journal of Differential Equations , vol. 188, no. 1, pp. 135–163, 2003. D. Xiao and S. Ruan, “ Global analysis of an epidemic model with nonmonotone incidence rate ,” Mathematical Biosciences , vol. 208, no. 2, pp. 419–429, 2007. Z. Yuan and L. Wang, “ Global stability of epidemiological models with group mixing and nonlinear incidence rates ,” Nonlinear Analysis: Real World Applications , vol. 11, no. 2, pp. 995–1004, 2010. Y. Tang, D. Huang, S. Ruan, and W. Zhang, “ Coexistence of limit cycles and homoclinic loops in a SIRS model with a nonlinear incidence rate ,” SIAM Journal on Applied Mathematics , vol. 69, no. 2, pp. 621–639, 2008. M. Y. Li and J. S. Muldowney, “ Global stability for the SEIR model in epidemiology ,” Mathematical Biosciences , vol. 125, no. 2, pp. 155–164, 1995. H. R. Thieme and C. Castillo-Chavez, “On the role of variable infectivity in the dynamics of the human immunodeficiency virus epidemic,” in Mathematical and Statistical Approaches to AIDS Epidemiology , vol. 83 of Lecture Notes in Biomathematics , Springer, Berlin, Germany, 1989. J. Zhang and Z. Ma, “ Global dynamics of an SEIR epidemic model with saturating contact rate ,” Mathematical Biosciences , vol. 185, no. 1, pp. 15–32, 2003. J. A. P. Heesterbeek and J. A. J. Metz, “ The saturating contact rate in marriage- and epidemic models ,” Journal of Mathematical Biology , vol. 31, no. 5, pp. 529–539, 1993. S. Bowong and J. J. Tewa, “ Global analysis of a dynamical model for transmission of tuberculosis with a general contact rate ,” Communications in Nonlinear Science and Numerical Simulation , vol. 15, no. 11, pp. 3621–3631, 2010. P. van den Driessche and J. Watmough, “ Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission ,” Mathematical Biosciences , vol. 180, pp. 29–48, 2002. C. Castillo-Chavez and H. R. Thieme, “Asymptotically autonomous epidemic models,” in Mathematical Population Dynamics: Analysis of Heterogeneity , O. Arino, et al., Ed., vol. 1 of Theory of Epidemics , pp. 33–50, Wuetz, 1995. K. Mischaikow, H. Smith, and H. R. Thieme, “ Asymptotically autonomous semiflows: chain recurrence and Lyapunov functions ,” Transactions of the American Mathematical Society , vol. 347, no. 5, pp. 1669–1685, 1995. H. R. Thieme, “ Persistence under relaxed point-dissipativity (with application to an endemic model) ,” SIAM Journal on Mathematical Analysis , vol. 24, no. 2, pp. 407–435, 1993. J. P. LaSalle, The Stability of Dynamical Systems , SIAM, Philadelphia, Pa, USA, 1976, With an appendix: “Limiting equations and stability of nonautonomous ordinary differential equations” by Z. Artstein, Regional Conference Series in Applied Mathematics. J. P. LaSalle, “ Stability theory for ordinary differential equations ,” Journal of Differential Equations , vol. 4, pp. 57–65, 1968. A. Korobeinikov and P. K. Maini, “A Lyapunov function and global properties for SIR and SEIR epidemiological models with nonlinear incidence,” Mathematical Biosciences and Engineering , vol. 1, no. 1, pp. 57–60, 2004. //

Journal

Abstract and Applied AnalysisHindawi Publishing Corporation

Published: Feb 21, 2011

There are no references for this article.