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

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 𝑅₁ and 𝑅₂, 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, 𝑅₁ and 𝑅₂, 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, the unique drug-resistant-TB-strain-only equilibrium exists and is globally asymptotically stable when 𝑅₁<1<𝑅₂. 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 (𝐸₁), individuals with symptoms of TB and drug-sensitive (𝐼₁), those exposed to drug-resistant TB (𝐸₂) and those displaying symptoms of TB and drug-resistant (𝐼₂), the total number of population at time 𝑡 is given by𝑁(𝑡)=𝑆(𝑡)+𝐸₁(𝑡)+𝐼₁(𝑡)+𝐸₂(𝑡)+𝐼₂(𝑡). (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 𝐶² function of the total population 𝑁. Susceptible are infected with drug-sensitive strain and drug-resistant strain entering classes 𝐸₁ and 𝐸₂, respectively. For simplicity, we only account for treatment of sensitive strain, we assume that chemoprophylaxis of individuals in 𝐸₁ reduces their reactivation at a constant rate 𝑟₁ and that the initiation of therapeutics immediately removes individuals in 𝐼₁ from active status and places them into exposed class at a constant rate 𝑟₂. The time before individuals in 𝐸₁ who does not received effective chemoprophylaxis become infectious is assumed to satisfy an exponential distribution, with mean waiting time 1/𝑘₁. Thus, individuals leave the class 𝐸₁ to the class 𝐼₁ at a constant rate 𝑘₁(1−𝑟₁). The time before individuals in 𝐸₂ become infectious is also assumed to satisfy an exponential distribution, with mean waiting time 1/𝑘₂. Thus, individuals leave the class 𝐸₂ to the class 𝐼₂ at a constant rate 𝑘₂. Individuals sick with drug-sensitive strain receive treatment at rate 𝑟₂ and a proportion 𝑝𝑟₂ respond to treatment and move into exposed class 𝐸₁, in the remaining proportion (1−𝑝)𝑟₂, inappropriate treatment results in the development of drug-resistant strain and the individuals move into class 𝐸₂. Individuals in the infectious class 𝐼₁ and 𝐼₂ suffer additional disease-induced death at rate 𝑑₁ and 𝑑₂, where 𝑑₂>𝑑₁. Then the following two strains tuberculosis model is formulated𝑆^′=Λ−𝛽₁(𝑁)𝑆𝐼₁−𝛽₂(𝑁)𝑆𝐼₂𝐸−𝜇𝑆,_′1=𝛽₁(𝑁)𝑆𝐼₁+𝑝𝑟₂𝐼₁−𝑘₁1−𝑟₁𝐸₁−𝜇𝐸₁,𝐼^′₁=𝑘₁1−𝑟₁𝐸₁−𝑟₂+𝜇+𝑑₁𝐼₁,𝐸^′₂=𝛽₂(𝑁)𝑆𝐼₂+1−𝑝𝑟₂𝐼₁−𝜇+𝑘₂𝐸₂,𝐼^′₂=𝑘₂𝐸₂−𝜇+𝑑₂𝐼₂, (2.2) where Λ, 𝑘₁, 𝑘₂, 𝑟₁, and 𝑟₂ are assumed to be positive and 𝑝∈[0,1]. It is natural to assume that the transmission coefficient 𝛽₁(𝑁),𝛽₂(𝑁) satisfies the following conditions:𝛽₁(𝑁)>0,𝛽^′₁(𝑁)≤0,𝑁𝛽₁(𝑁)^′𝛽≥0,₂(𝑁)>0,𝛽^′₂(𝑁)≤0,𝑁𝛽₂(𝑁)^′≥0. (2.3)

Figure 1 

Transfer diagram for the dynamics of tuberculosis with general contact rate.

Obviously, we can see that 𝛽₁(𝑁)=𝛽₁/𝑁, 𝛽₂(𝑁)=𝛽₂/𝑁 corresponds to the standard incidence rate, that 𝛽₁(𝑁)=𝛽₁, 𝛽₂(𝑁)=𝛽₂ corresponds to the mass action incidence rate, and that 𝛽₁(𝑁)=𝛽₁𝐶(𝑁), 𝛽₂(𝑁)=𝛽₂𝐶(𝑁) 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 𝑆, 𝐸₁, 𝐼₁, 𝐸₂, 𝐼₂ 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, 𝐸₁(0)≥0, 𝐼₁(0)≥0, 𝐸₂(0)≥0, 𝐼₂(0)≥0, the solutions 𝑆(𝑡), 𝐸₁(𝑡), 𝐼₁(𝑡), 𝐸₂(𝑡), 𝐼₂(𝑡) 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 Γ=𝑆(𝑡),𝐸₁(𝑡),𝐼₁(𝑡),𝐸₂(𝑡),𝐼₂(𝑡)∈ℝ⁵₊(𝑡)≤Λ∶𝑁𝜇. (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 𝑡₁ such that 𝑆𝑡₁=0,𝑆^′𝑡₁<0,𝐸₁(𝑡)≥0,𝐼₁(𝑡)≥0,𝐸₂(𝑡)≥0,𝐼₂(𝑡)≥0,for0<𝑡<𝑡₁, (2.6) there exists a 𝑡₂𝐸₁𝑡₂=0,𝐸^′₁𝑡₂<0,𝑆(𝑡)≥0,𝐼₁(𝑡)≥0,𝐸₂(𝑡)≥0,𝐼₂(𝑡)≥0,for0<𝑡<𝑡₂, (2.7) there exists a 𝑡₃𝐼₁𝑡₃=0,𝐼^′₁𝑡₃<0,𝑆(𝑡)≥0,𝐸₁(𝑡)≥0,𝐸₂(𝑡)≥0,𝐼₂(𝑡)≥0,for0<𝑡<𝑡₃, (2.8) there exists a 𝑡₄𝐸₂𝑡₄=0,𝐸^′₂𝑡₄<0,𝑆(𝑡)≥0,𝐸₁(𝑡)≥0,𝐼₁(𝑡)≥0,𝐼₂(𝑡)≥0,for0<𝑡<𝑡₄, (2.9) or there exists a 𝑡₅𝐼₂𝑡₅=0,𝐼^′₂𝑡₅<0,𝑆(𝑡)≥0,𝐸₁(𝑡)≥0,𝐼₁(𝑡)≥0,𝐸₂(𝑡)≥0,for0<𝑡<𝑡₅. (2.10) In the first case, we have 𝑆^′𝑡₁=Λ>0, (2.11) which is a contradiction meaning that 𝑆(𝑡)≥0, 𝑡≥0. In the second case, we have 𝐸^′₁𝑡₂=𝐼₁𝑡₂𝛽₁(𝑁)𝑆𝑡₂+𝑝𝑟₂≥0, (2.12) which is a contradiction meaning that 𝐸₁(𝑡)≥0, 𝑡≥0. In the third case, we have 𝐼^′₁𝑡₃=𝑘₁1−𝑟₁𝐸₁𝑡₃≥0, (2.13) which is a contradiction meaning that 𝐼₁(𝑡)≥0, 𝑡≥0. In the fourth case, we have 𝐸^′₂𝑡₄=𝛽₂(𝑁)𝑆𝑡₄𝐼₂𝑡₄+1−𝑝𝑟₂𝐼₁𝑡₄≥0, (2.14) which is a contradiction meaning that 𝐸₂(𝑡)≥0, 𝑡≥0. In the fifth case, we have 𝐼^′₂𝑡₅=𝑘₂𝐸₂𝑡₅≥0, (2.15) which is a contradiction meaning that 𝐼₂(𝑡)≥0, 𝑡≥0.
Thus, in all cases, 𝑆(𝑡), 𝐸₁(𝑡), 𝐼₁(𝑡), 𝐸₂(𝑡), 𝐼₂(𝑡) remain positive for 𝑡≥0.
Let (𝑆(𝑡),𝐸₁(𝑡),𝐼₁(𝑡),𝐸₂(𝑡),𝐼₂(𝑡))∈ℝ⁵₊ be any solution with nonnegative initial condition, adding all equations in system (2.2) gives 𝑁^′(𝑡)=Λ−𝜇𝑁−𝑑₁𝐼₁−𝑑₂𝐼₂,≤Λ−𝜇𝑁. (2.16) It follows that (𝑡)≤Λ0≤𝑁𝜇+𝑁(0)𝑒^−𝜇𝑡, (2.17) where 𝑁(0) represents initial values of the total population. Thus 0≤𝑁(𝑡)≤Λ/𝜇, as 𝑡→∞. Therefore all feasible solutions of system (2.2) enter the region Γ=𝑆(𝑡),𝐸₁(𝑡),𝐼₁(𝑡),𝐸₂(𝑡),𝐼₂(𝑡)∈ℝ⁵₊(𝑡)≤Λ∶𝑁𝜇. (2.18) 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 𝑃₀ and two possible endemic equilibria for system (2.2), the endemic equilibria include boundary equilibrium 𝑃₁ (when only the drug-resistant TB-strain is present) and the interior equilibrium 𝑃₂ (when both strains exist).

3.1. Global Stability of the Disease-Free Equilibrium

The disease-free-equilibrium is given as𝑃₀=𝑆₀,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 𝑃₀, we use the results in van den Driessche and Watmough [27]. Reorder the components of 𝑃₀ as 𝐸₁=0, 𝐼₁=0, 𝐸₂=0, 𝐼₂=0, 𝑆=Λ/𝜇. Set𝐹ℱ=¹𝐹₂𝐹₃𝐹₄𝐹₅=𝛽₁(𝑁)𝑆𝐼₁0𝛽₂(𝑁)𝑆𝐼₂00,𝑘𝒱=¹1−𝑟₁𝐸₁+𝜇𝐸₁−𝑝𝑟₂𝐼₁𝑟₂𝐼₁+𝜇+𝑑₁𝐼₁−𝑘₁1−𝑟₁𝐸₁𝜇+𝑘₂𝐸₂−1−𝑝𝑟₂𝐼₁𝜇+𝑑₂𝐼₂−𝑘₂𝐸₂𝜇𝑆+𝛽₁(𝑁)𝑆𝐼₁+𝛽₂(𝑁)𝑆𝐼₂.−Λ (3.2) The infected compartments are 𝐸₁, 𝐼₁, 𝐸₂, and 𝐼₂. Thus𝐹=0𝛽₁𝑆₀𝑆₀000000000𝛽₁𝑆₀𝑆₀,0000𝑉=𝜇+𝑘₁1−𝑟₁−𝑝𝑟₂00−𝑘₁1−𝑟₁𝜇+𝑑₁+𝑟₂𝑟000−1−𝑝₂𝜇+𝑘₂000−𝑘₂𝜇+𝑑₂. (3.3) The dominant eigenvalues of 𝐹𝑉⁻¹ are given by 𝑅₁=𝑘₁1−𝑟₁𝛽₁𝑆₀𝑆₀𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁𝑝𝑟₂1−𝑟₁,𝑅₂=𝑘₂𝛽₂𝑆₀𝑆₀𝜇+𝑘₂𝜇+𝑑₂, (3.4) where 𝑅₁ and 𝑅₂ 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 𝐹𝑉⁻¹ is𝜌𝐹𝑉⁻¹=max𝑅₁,𝑅₂. (3.5) If 𝑅₁<1, and 𝑅₂<1, then 𝜌(𝐹𝑉⁻¹)<1. By Theorem 3.2 in van den Driessche and Watmough [27], we know that the disease-free steady state 𝑃₀ is locally asymptotically stable, 𝑃₀ is unstable if 𝑅₁>1 or 𝑅₂>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, 𝑑₁=0, 𝑑₂=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 𝑁=𝑆₀=Λ/𝜇: 𝐸^′₁=𝛽₁𝑆₀𝐼₁𝑆₀−𝐸₁−𝐼₁−𝐸₂−𝐼₂+𝑝𝑟₂𝐼₁−𝜇+𝑘₁1−𝑟₁𝐸₁,𝐼^′₁=𝑘₁1−𝑟₁𝐸₁−𝑟₂+𝜇𝐼₁,𝐸^′₂=𝛽₂𝑆₀𝐼₂𝑆₀−𝐸₁−𝐼₁−𝐸₂−𝐼₂+1−𝑝𝑟₂𝐼₁−𝜇+𝑘₂𝐸_2,𝐼^′₂=𝑘₂𝐸₂−𝜇𝐼₂. (3.6) Notice that the 𝑆 equation is eliminated from (3.6) and the variable 𝑆 is replaced by 𝑆₀−𝐸₁−𝐼₁−𝐸₂−𝐼₂. By introducing the fractions𝑥₁=𝐸₁𝑁,𝑥₂=𝐼₁𝑁,𝑦₁=𝐸₂𝑁,𝑦₂=𝐼₂𝑁, (3.7) we obtain the equivalent limiting system from (3.6) 𝑥^′₁=𝑆₀𝛽₁𝑆₀1−𝑥₁−𝑥₂−𝑦₁−𝑦₂𝑥₂+𝑝𝑟₂𝑥₂−𝜇+𝑘₁1−𝑟₁𝑥₁,𝑥^′₂=𝑘₁1−𝑟₁𝑥₁−𝜇+𝑟₂𝑥₂,𝑦^′₁=𝑆₀𝛽₂𝑆₀1−𝑥₁−𝑥₂−𝑦₁−𝑦₂𝑦₂+1−𝑝𝑟₂𝑥₂−𝜇+𝑘₂𝑦₁,𝑦^′₂=𝑘₂𝑦₁−𝜇𝑦₂. (3.8) Obviously0≤𝑥₁+𝑥₂+𝑦₁+𝑦₂≤1, (3.9) for all time 𝑡≥0.

For a bounded real-valued function 𝑓 on [0,∞), we define𝑓^∞=liminf_𝑡→∞𝑓(𝑡),𝑓_∞=limsup_𝑡→∞𝑓(𝑡). (3.10)

Lemma 3.1 (see [30]). Let 𝑓∶[0,∞)→∞ be bounded and continuously differentiable. Then there are sequences 𝑡_𝑛,𝑠_𝑛→∞ with the following properties: 𝑓𝑡_𝑛⟶𝑓_∞,𝑓^′𝑡_𝑛𝑓𝑠⟶0,_𝑛⟶𝑓^∞,𝑓^′𝑠_𝑛⟶0, (3.11) for 𝑛→∞.

Theorem 3.2. In the absence of disease-induced death rate, the disease-free equilibrium 𝑃₀ of system (2.2) is globally asymptotically stable if 𝑅₁<1, 𝑅₂<1.

Proof. By Lemma 3.1 and the 𝑥^′₂, 𝑦^′₂ equations in (3.8) we have 𝑥^∞₂≤𝑘₁1−𝑟₁𝜇+𝑟₂𝑥^∞₁,𝑦^∞₂≤𝑘₂𝜇𝑦^∞₁. (3.12)
Using the 𝑥^′₁ equation in (3.8) and choosing 𝑡_𝑛→∞ such that𝑥₁𝑡_𝑛⟶𝑥^∞₁,𝑥^′₁𝑡_𝑛⟶0,𝑡⟶∞, (3.13) we get 0≤𝑆₀𝛽₁𝑆₀1−𝑥₁−𝑥₂−𝑦₁−𝑦₂^∞𝑥^∞₂−𝜇+𝑘₁1−𝑟₁𝑥^∞₁≤𝑆₀𝛽₁𝑆₀𝑥^∞₂−𝜇+𝑘₁1−𝑟₁𝑥^∞₁≤𝑆₀𝛽₁𝑆₀𝑘₁1−𝑟₁𝜇+𝑟₂𝑥^∞₁−𝜇+𝑘₁1−𝑟₁𝑥^∞₁≤𝜇+𝑘₁1−𝑟₁𝑥^∞₁𝑘₁1−𝑟₁𝑆₀𝛽₁𝑆₀𝜇+𝑟₂𝜇+𝑘₁1−𝑟₁≤−1𝜇+𝑘₁1−𝑟₁𝑥^∞₁𝑘₁1−𝑟₁𝑆₀𝛽₁𝑆₀𝜇+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁𝑝𝑟₂1−𝑟₁=−1𝜇+𝑘₁1−𝑟₁𝑥^∞₁𝑅₁−1. (3.14) It is shown that 𝑥^∞₁≤0 since 𝑅₁<1, as 𝑥^∞₁≥0, we have that 𝑥^∞₁=0. The inequalities in (3.12) also imply that 𝑥^∞₂=0.
Using the 𝑦^′₁ equation in (3.8) and choosing 𝑠_𝑛→∞ such that𝑦₁𝑠_𝑛⟶𝑦^∞₁,𝑦^′₁𝑠_𝑛⟶0,𝑡→∞, (3.15) we get 0≤𝑆₀𝛽₂𝑆₀1−𝑥₁−𝑥₂−𝑦₁−𝑦₂^∞𝑦^∞₂−𝜇+𝑘₂𝑦^∞₁≤𝑆₀𝛽₂𝑆₀𝑦^∞₂−𝜇+𝑘₂𝑦^∞₁≤𝑆₀𝛽₂𝑆₀𝑘₂𝜇𝑦^∞₁−𝜇+𝑘₂𝑦^∞₁≤𝜇+𝑘₂𝑦^∞₁𝑘₂𝑆₀𝛽₂𝑆₀𝜇𝜇+𝑘₂=−1𝜇+𝑘₂𝑦^∞₁𝑅₂−1. (3.16) It is shown that 𝑦^∞₁≤0 since 𝑅₂<1, as 𝑦^∞₁≥0, we have that 𝑦^∞₁=0. The inequalities in (3.12) also imply that 𝑦^∞₂=0.
Hencelim_𝑡→∞𝑥₁(𝑡)=lim_𝑡→∞𝑥₂(𝑡)=lim_𝑡→∞𝑦₁(𝑡)=lim_𝑡→∞𝑦₂(𝑡)=0. (3.17) Therefore, the disease-free equilibrium 𝑃₀ is globally asymptotically stable.

Next, we consider that there is disease-induced death rate, which means that 𝑑₁>0, 𝑑₂>0. Then, we can get the following theorem.

Theorem 3.3. For 𝑑₁>0, 𝑑₂>0, the disease-free equilibrium 𝑃₀ of system (2.2) is globally asymptotically stable if 𝑅₁<1−𝑘₂𝑟₂(1−𝑝)/((𝜇+𝑑₁+𝑟₂)(𝜇+𝑘₁(1−𝑟₁))−𝑘₁𝑝𝑟₂(1−𝑟₁)), 𝑅₂<1.

Proof. Let us consider the following Lyapunov function: 𝑉𝐸₁,𝐼₁,𝐸₂,𝐼₂=𝑘₁1−𝑟₁𝐸₁+𝜇+𝑘₁1−𝑟₁𝐼₁+𝑘₂𝐸₂+𝜇+𝑘₂𝐼₂. (3.18) Its time derivative along the solutions of system (2.2) satisfies 𝑉̇=𝑘₁1−𝑟₁𝐸̇₁+𝜇+𝑘₁1−𝑟₁𝐼̇₁+𝑘₂𝐸̇₂+𝜇+𝑘₂𝐼̇₂=𝑘₁1−𝑟₁𝛽₁(𝑁)𝑆𝐼₁+𝑝𝑟₂𝐼₁−𝜇+𝑘₁1−𝑟₁𝐸₁+𝜇+𝑘₁1−𝑟₁𝑘₁1−𝑟₁𝐸₁−𝜇+𝑑₁+𝑟₂𝐼₁+𝑘₂𝛽₂(𝑁)𝑆𝐼₂+1−𝑝𝑟₂𝐼₁−𝜇+𝑘₂𝐸₂+𝜇+𝑘₂𝑘₂𝐸₂−𝜇+𝑑₂𝐼₂,=𝑘₁1−𝑟₁𝛽₁(𝑁)𝑆+𝑘₁1−𝑟₁𝑝𝑟₂−𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁𝐼₁+𝑘₂𝛽₂(𝑁)𝑆−𝜇+𝑑₂𝜇+𝑘₂𝐼₂+𝑘₂𝑟₂1−𝑝𝐼₁. (3.19) Now, using (2.3), one has 𝛽₁(𝑁)𝑆≤𝛽₁(𝑆₀)𝑆₀, 𝛽₂(𝑁)𝑆≤𝛽₂(𝑆₀)𝑆₀, then 𝑉̇≤𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁1−𝑟₁𝑝𝑟₂𝑅₁−1𝐼₁+𝑘₂𝑟₂1−𝑝𝐼₁+𝜇+𝑑₂𝜇+𝑘₂𝑅₂−1𝐼₂,=𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁1−𝑟₁𝑝𝑟₂𝑅₁−1+𝑘₂𝑟₂1−𝑝𝐼₁+𝜇+𝑑₂𝜇+𝑘₂𝑅₂−1𝐼₂. (3.20) Since 𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁𝑝𝑟₂1−𝑟₁=𝜇𝜇+𝑑₁+𝑟₂+𝑘₁1−𝑟₁𝜇+𝑑₁+𝑘₁𝑟₂1−𝑟₁−𝑘₁𝑝𝑟₂1−𝑟₁=𝜇𝜇+𝑑₁+𝑟₂+𝑘₁1−𝑟₁𝜇+𝑑₁+𝑘₁𝑟₂1−𝑟₁1−𝑝≥0, (3.21) and 𝑅₁<1−𝑘₂𝑟₂(1−𝑝)/((𝜇+𝑑₁+𝑟₂)(𝜇+𝑘₁(1−𝑟₁))−𝑘₁𝑝𝑟₂(1−𝑟₁))≤1, we have 𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁1−𝑟₁𝑝𝑟₂𝑅₁−1+𝑘₂𝑟₂1−𝑝≤0. (3.22) Thus, 𝑉̇≤0 if 𝑅₁<1−𝑘₂𝑟₂(1−𝑝)/((𝜇+𝑑₁+𝑟₂)(𝜇+𝑘₁(1−𝑟₁))−𝑘₁𝑝𝑟₂(1−𝑟₁)), 𝑅₂<1. It is obvious that 𝑉̇=0 if and only if 𝐼₁=0, 𝐼₂=0. Then, the largest invariant set of system (2.2) on the set {(𝑆,𝐸₁,𝐼₁,𝐸₂,𝐼₂)∈ℝ⁵₊, 𝑉̇(𝐸₁,𝐼₁,𝐸₂,𝐼₂)=0} is the disease-free equilibrium 𝑃₀. Therefore, it follows from LaSalle's invariance principle [31, 32] that 𝑃₀ is globally stable if 𝑅₁<1−𝑘₂𝑟₂(1−𝑝)/((𝜇+𝑑₁+𝑟₂)(𝜇+𝑘₁(1−𝑟₁))−𝑘₁𝑝𝑟₂(1−𝑟₁)), 𝑅₂<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 𝑅₁ and 𝑅₂.

Let 𝑃₁=(𝑆,0,0,𝐸₂,𝐼₂̂) 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 Λ−𝛽₂𝑁𝑆𝐼₂̂𝛽−𝜇𝑆=0,₂𝑁𝑆𝐼₂̂−𝜇+𝑘₂𝐸₂𝑘=0,₂𝐸₂−𝜇+𝑑₂𝐼₂̂=0, (3.23) adding the above three equations, we haveΛ−𝜇𝑁−𝑑₂𝐼₂̂=0. (3.24)

Using (3.24), the first and second equation of (3.23), we can easily express 𝑆, 𝐸₂ and 𝐼₂̂ in terms of 𝑁 in the form: 𝑆=Λ𝑑₂𝜇𝑑₂+𝛽₂𝑁𝑁,𝐸Λ−𝜇₂=Λ𝛽₂𝑁Λ−𝜇𝑁𝜇+𝑘₂𝜇𝑑₂+𝛽₂𝑁𝑁,𝐼̂Λ−𝜇₂=Λ−𝜇𝑁𝑑₂. (3.25)

Substituting 𝑆, 𝐸₂, 𝐼₂̂ in the third equations of (3.23) yieldsΛ−𝜇𝑁𝐹𝑁=0, (3.26) where𝐹𝑁=Λ𝑘₂𝑑₂𝛽₂𝑁−𝜇+𝑘₂𝜇+𝑑₂𝜇𝑑₂+𝛽₂𝑁Λ−𝜇𝑁. (3.27)

Clearly, Λ−𝜇𝑁=0 is a fixed point of (3.23), which corresponds to the disease-free equilibrium 𝑃₀. Since 𝑁∈[0,𝑆₀], we get𝐹(0)=Λ𝑘₂𝑑₂𝛽₂(0)−𝜇+𝑘₂𝜇+𝑑₂𝜇𝑑₂+Λ𝛽₂(0)=Λ𝑘₂𝑑₂𝛽₂(0)−𝜇²+𝜇𝑑₂₊𝑘₂+𝑑₂𝑘₂𝜇𝑑₂+Λ𝛽₂(0)=Λ𝑘₂𝑑₂𝛽₂(0)−𝜇²+𝜇𝑑₂₊𝑘₂𝜇𝑑₂+Λ𝛽₂(0)−Λ𝑑₂𝑘₂𝛽₂(0)−𝜇𝑘₂𝑑²₂=−𝜇²+𝜇𝑘₂+𝑑₂𝜇𝑑₂+Λ𝛽₂(0)−𝜇𝑘₂𝑑²₂,𝐹𝑆₀=Λ𝑘₂𝑑₂𝛽₂𝑆₀−𝜇𝑑₂𝜇+𝑘₂𝜇+𝑑₂=𝜇𝑘₂𝑑₂𝛽₂𝑆₀𝑆₀−𝑑₂𝜇+𝑘₂𝜇+𝑑₂=𝜇𝑑₂𝜇+𝑘₂𝜇+𝑑₂𝑘₂𝛽₂𝑆₀𝑆₀𝜇+𝑘₂𝜇+𝑑₂−1=𝜇𝑑₂𝜇+𝑘₂𝜇+𝑑₂𝑅₂−1. (3.28)

It appears that 𝐹(0)<0 and 𝐹(𝑆₀)>0 when 𝑅₂>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,𝑆₀].

Thus, we have established the following result

Lemma 3.4. When𝑅₂>1, the system (2.2) has a unique drug-resistant-TB-strain-only equilibrium 𝑃₁=(𝑆,0,0,𝐸₂,𝐼₂̂) with 𝑆, 𝐸₂, and 𝐼₂̂ all nonnegative.

3.3. Global Stability of the Drug-Resistant-TB-Strain-Only Equilibrium

Theorem 3.5. If 𝑅₁<1<𝑅₂, the drug-resistant-TB-strain-only equilibrium 𝑃₁ exists and is locally asymptotically stable.

Proof. Reorder the equilibrium 𝑃₁ as (0,0,𝐸₂,𝐼₂̂,𝑆). Similarly as in the proof of locally stability of the disease-free steady state 𝑃₀, we have 𝐹₁=0𝛽₁𝑁𝑆,𝑉00₁=𝜇+𝑘₁1−𝑟₁−𝑝𝑟₂−𝑘₁1−𝑟₁𝜇+𝑑₁+𝑟₂,𝐹₁𝑉₁⁻¹𝑘=𝐴¹1−𝑟₁𝛽₁𝑁𝑆𝜇+𝑘₁1−𝑟₁𝛽₁𝑁𝑆,00 (3.29) where 1𝐴=𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁𝑝𝑟₂1−𝑟₁. (3.30) The spectral radius of 𝐹₁𝑉₁⁻¹ is given by 𝜌𝐹₁𝑉₁⁻¹=𝑘₁1−𝑟₁𝛽₁𝑁𝑆𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁𝑝𝑟₂1−𝑟₁,≤𝑘₁1−𝑟₁𝛽₁𝑁𝑁𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁𝑝𝑟₂1−𝑟₁. (3.31) Since (𝛽₁(𝑁)𝑁)^′≥0, and 𝑁≤𝑆₀, then 𝛽₁𝑁𝑁≤𝛽₁𝑆₀𝑆₀, (3.32) Thus 𝜌𝐹₁𝑉₁⁻¹≤𝑘₁1−𝑟₁𝛽₁𝑆₀𝑆₀𝜇+𝑑₁+𝑟₂𝜇+𝑘₁1−𝑟₁−𝑘₁𝑝𝑟₂1−𝑟₁=𝑅₁. (3.33) As 𝑅₁<1, we have 𝜌(𝐹₁𝑉₁⁻¹)<1, which implies that the drug-resistant-TB-strain-only equilibrium 𝑃₁ 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 𝑃₁ is global asymptotically stable when 𝑅₁<1<𝑅₂.

Proof. If 𝑅₁<1<𝑅₂, from Theorem 3.5, 𝑃₁ is locally asymptotically stable. In the following, we only need to prove that the 𝑃₁ is a global attractor.
There is no disease-induced death rate, which means 𝑑₁=𝑑₂=0, using a similar argument as in the proof of Theorem 3.2, we can show that if 𝑅₁<1,lim_𝑡→∞𝑥₁(𝑡)=lim_𝑡→∞𝑥₂(𝑡)=0. (3.34) Then the equivalent limiting system of (3.6) is 𝐸^′₂=𝛽₂𝑆₀𝐼₂𝑆₀−𝐸₂−𝐼₂−𝜇+𝑘₂𝐸₂,𝐼^′₂=𝑘₂𝐸₂−𝜇𝐼₂. (3.35)
Let (𝑃,𝑄) be the vector field defined by system (3.35). Then for the Dulac function 𝐷(𝐸₂,𝐼₂)=1/𝐸₂𝐼₂, there holds𝜕𝐸₂()+𝜕𝐷𝑃𝐼₂()𝛽𝐷𝑄=−²𝑆₀𝑆₀−𝐼₂𝐸²₂−𝑘₂𝐼²₂<0. (3.36)
From Dulac's criterion, we can see equivalent limiting system (3.35) does not have a limit cycle. Therefore, the local stability of 𝑃₁ 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₁)lim_𝑁→0(1/𝛽₁(𝑁))=+∞, (W₂)𝛽₁(𝑁)/𝛽₂(𝑁)=𝛽₁/𝛽₂,

then, existence, uniqueness and global stability of interior equilibrium can be obtained. Let 𝑃₂=(𝑆^∗,𝐸^∗₁,𝐼^∗₁,𝐸^∗₂,𝐼^∗₂) be interior equilibrium of system (2.2). Then, the interior equilibrium (steady state with 𝐼₁>0, 𝐼₂>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 Λ−𝛽₁𝑁^∗𝑆^∗𝐼^∗₁−𝛽₂𝑁^∗𝑆^∗𝐼^∗₂−𝜇𝑆^∗𝛽=0,₁𝑁^∗𝑆^∗𝐼^∗₁+𝑝𝑟₂𝐼^∗₁−𝑘₁1−𝑟₁𝐸^∗₁−𝜇𝐸^∗₁𝑘=0,₁1−𝑟₁𝐸^∗₁−𝑟₂+𝜇+𝑑₁𝐼^∗₁𝛽=0,₂𝑁^∗𝑆^∗𝐼^∗₂+1−𝑝𝑟₂𝐼^∗₁−𝜇+𝑘₂𝐸^∗₂𝑘=0,₂𝐸^∗₂−𝜇+𝑑₂𝐼^∗₂=0, (3.37) adding the above five equations, we haveΛ−𝜇𝑁^∗−𝑑₁𝐼^∗₁−𝑑₂𝐼^∗₂=0. (3.38)

Using (3.37) and (3.38), we get the following equivalent equations set Λ−𝛽₁𝑁^∗𝑆^∗𝐼^∗₁−𝛽₂𝑁^∗𝑆^∗𝐼^∗₂−𝜇𝑆^∗𝛽=0,₁𝑁^∗𝑆^∗𝐼^∗₁+𝑝𝑟₂𝐼^∗₁−𝑘₁1−𝑟₁𝐸^∗₁−𝜇𝐸^∗₁𝑘=0,₁1−𝑟₁𝐸^∗₁−𝑟₂+𝜇+𝑑₁𝐼^∗₁𝑘=0,₂𝐸^∗₂−𝜇+𝑑₂𝐼^∗₂=0,Λ−𝜇𝑁^∗−𝑑₁𝐼^∗₁−𝑑₂𝐼^∗₂=0. (3.39) From (3.39), we have 𝑆^∗=𝛽𝑆₀𝑆₀𝑅₁𝛽₁(𝑁^∗),𝐼^∗₁=𝑑₂Λ−𝜇𝑆^∗−𝑆^∗𝛽₂𝑁^∗Λ−𝜇𝑁^∗𝑑₂𝑆^∗𝛽₁(𝑁^∗)−𝑑₁𝑆^∗𝛽₂(𝑁^∗),𝐸^∗₁=𝜇+𝑑₁+𝑟₂𝑘₁1−𝑟₁𝐼^∗₁,𝐼^∗₂=Λ−𝜇𝑁^∗𝑑₂−𝑑₁𝑑₂𝐼^∗₁,𝐸^∗₂=𝜇+𝑑₂Λ−𝜇𝑁^∗𝑘₂𝑑₂−𝑑₁𝜇+𝑑₂𝑘₂𝑑₂𝐼^∗₁. (3.40)

Substituting (3.40) into 𝛽₂(𝑁^∗)𝑆^∗𝐼^∗₂+(1−𝑝)𝑟₂𝐼^∗₁−(𝜇+𝑘₂)𝐸^∗₂=0, 𝑁^∗satisfies the following equation:𝐺𝑁^∗=0, (3.41) where 𝐺𝑁^∗=𝛽₁𝑆₀𝑆₀Λ−𝜇𝑁^∗𝑅₁𝑘₂𝑑₂𝑆₀𝛽₁𝑆₀𝛽₂𝑁^∗𝑅₁𝛽₁(𝑁^∗)−𝑑₂𝜇₁−1−𝑝𝑘₂𝑑₂𝑟₂𝛽₂𝑁^∗𝛽₁(𝑁^∗)+𝑑₂Λ𝛽1−¹𝑆₀𝛽₁(𝑁^∗)𝑅₁𝑑₁𝜇₁+1−𝑝𝑘₂𝑑₂𝑟₂−𝑑₁𝑘₂𝑆₀𝛽₁𝑆₀𝛽₂𝑁^∗𝑅₁𝛽₁(𝑁^∗),𝜇₁=𝜇+𝑘₂𝜇+𝑑₂. (3.42) Since 𝛽₁(𝑁)/𝛽₂(𝑁)=𝛽₁/𝛽₂, lim_𝑁→0(1/𝛽₁(𝑁))=+∞, if 𝑅₁>𝑑₁𝑘₂𝛽₂𝑆₀𝛽₁(𝑆₀)/𝛽₁[𝑑₁𝜇₁+(1−𝑝)𝑘₂𝑑₂𝑟₂], thenlim_𝑁∗→0𝐺𝑁^∗=−∞. (3.43) It means that 𝐺(𝑁^∗)=0 has only one positive root in the interval [0,𝑆₀] if 𝐺(𝑁^∗) is monotone increasing and 𝐺(𝑆₀)≥0, which is equivalent to 𝑅₁>𝐻₁+𝐻₂𝐻₃+𝐻₄,𝑅₂<𝑑1+1−𝑝₂𝑘₂𝑟₂𝑑₁𝜇₁𝑅₁, (3.44) where 𝐻₁=𝑘₂𝑑₂𝛽₂𝑆₀𝛽₁𝑆₀𝛽²₁𝑁^∗,𝐻₂=𝑑₁𝑑₂𝑘₂𝛽₂𝑆₀𝛽₁𝑆₀𝛽^′₁𝑁^∗,𝐻₃=𝛽²₁𝑁^∗𝑑₂𝛽₁𝜇₁+𝑑₂𝑘₂𝛽₂𝑟₂,𝐻1−𝑝₄=𝑑₂𝛽^′₁𝑁^∗𝑑₁𝛽₁𝜇₁+𝑑₂𝑘₂𝛽₁𝑟₂1−𝑝. (3.45) Let 𝐻𝐻=max1,¹+𝐻₂𝐻₃+𝐻₄,𝑑₁𝑘₂𝛽₂𝑆₀𝛽₁𝑆₀𝛽₁𝑑₁𝜇₁+𝑘₂𝑑₂𝑟₂,𝑓𝑅1−𝑝₁=𝑑1+1−𝑝₂𝑘₂𝑟₂𝑑₁𝜇₁𝑅₁. (3.46) When𝑅₁>𝐻,𝑅₂<𝑓𝑅₁, (3.47) it shows that 𝐺(𝑆₀)≥0, and 𝐺(𝑁^∗) is monotone increasing, by the intermediate value-theorem, 𝐺(𝑁^∗)=0 has only one positive root in the interval [0,𝑆₀], which signify that 𝑃₂ is unique. Obviously, 𝑆^∗>0. In order to make certain that 𝐸^∗₁>0, 𝐼^∗₁>0, 𝐸^∗₂>0, 𝐼^∗₂>0, we obtain that 𝐸^∗₁>0, 𝐼^∗₁>0, 𝐸^∗₂>0, 𝐼^∗₂>0 if and only if0<𝐼^∗₁<Λ−𝜇𝑁^∗𝑑₁. (3.48) Due to 𝑑₂>𝑑₁>0, introducing (3.40) into (3.48), we acquire𝐵₁<𝑅₁<𝐵₂, (3.49) where 𝐵₁=𝑆₀𝛽₁𝑆₀𝛽₂𝑁^∗Λ−𝜇𝑁^∗+𝜇𝑑₂Λ𝑑₂𝛽₁(𝑁^∗),𝐵₂=𝑆₀𝛽₁𝑆₀𝛽₁𝑁^∗Λ−𝜇𝑁^∗+𝜇𝑑₁Λ𝑑₁𝛽₁(𝑁^∗), (3.50) and 𝐵₁<𝐵₂.

Now, we can get the following conclusion.

Lemma 3.7. When conditions (W₁), (W₂), (3.47) and (3.49) are satisfied, the system (2.2) has a unique interior equilibrium 𝑃₂=(𝑆^∗,𝐸^∗₁,𝐼^∗₁,𝐸^∗₂,𝐼^∗₂) with 𝑆^∗, 𝐸^∗₁, 𝐼^∗₁, 𝐸^∗₂, and 𝐼^∗₂ all nonnegative.

Remark 3.8. For 𝛽₁(𝑁)=𝛽₁ and 𝛽₁(𝑁)=𝛽₁/𝑁, condition (W₁) is not satisfied. However, at least when 𝛽₁(𝑁)=𝛽₁𝐶(𝑁), condition (W₁) 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 𝑃₂ 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, 𝛽₁(𝑁)=𝛽₂(𝑁)=𝛽(𝑁), then, we get the following result.

Theorem 3.9. If conditions (W₁), (W₂), (3.47) and (3.49) are satisfied, interior equilibrium 𝑃₂ of system (2.2) is is globally asymptotically stable when 𝑆𝑆^∗=𝐸₁𝐸^∗₁=𝐸₂𝐸^∗₂𝐼≥1,²𝐼^∗₂≥𝐼₁𝐼^∗₁𝑆or𝑆^∗=𝐸₁𝐸^∗₁=𝐸₂𝐸^∗₂𝐼≤1,²𝐼^∗₂≤𝐼₁𝐼^∗₁. (3.51)

Proof. Following [33], we consider the Lyapunov function: 𝑈𝑆,𝐸₁,𝐼₁,𝐸₂,𝐼₂=𝑆−𝑆^∗ln𝑆+𝐸₁−𝐸^∗₁ln𝐸^∗₁+𝐶𝐼₁−𝐼^∗₁ln𝐼^∗₁+𝐸₂−𝐸^∗₂ln𝐸^∗₂+𝐷𝐼₂−𝐼^∗₂ln𝐼^∗₂, (3.52) 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−^∗2𝐸₂𝐸̇₂+𝐷1−𝐼^∗₂𝐼₂𝐼̇₂=𝑆1−^∗𝑆Λ−𝛽(𝑁)𝑆𝐼₁−𝛽(𝑁)𝑆𝐼₂+𝐸−𝜇𝑆1−^∗1𝐸₁𝛽(𝑁)𝑆𝐼₁+𝑝𝑟₂𝐼₁−𝑘₁1−𝑟₁𝐸₁−𝜇𝐸₁𝐼+𝐶1−^∗1𝐼₁𝑘₁1−𝑟₁𝐸₁−𝑟₂+𝜇+𝑑₁𝐼₁+𝐸1−^∗2𝐸₂𝛽(𝑁)𝑆𝐼₂+1−𝑝𝑟₂𝐼₁−𝜇+𝑘₂𝐸₂𝐼+𝐷1−^∗2𝐼₂𝑘₂𝐸₂−𝜇+𝑑₂𝐼₂. (3.53) Considering (3.37), we can deduce that Λ=𝛽𝑁^∗𝑆^∗𝐼^∗₁+𝛽𝑁^∗𝑆^∗𝐼^∗₂+𝜇𝑆^∗,𝜇+𝑘₁1−𝑟₁=𝛽𝑁^∗𝑆^∗𝐼^∗₁𝐸^∗₁+𝑝𝑟₂𝐼^∗₁𝐸^∗₁,𝜇+𝑑₁+𝑟₂=𝑘₁1−𝑟₁𝐸^∗₁𝐼^∗₁,𝜇+𝑘₂=𝛽𝑁^∗𝑆^∗𝐼^∗₂𝐸^∗₂+1−𝑝𝑟₂𝐼^∗₁𝐸^∗₂,𝜇+𝑑₂=𝑘₂𝐸^∗₂𝐼^∗₂. (3.54) Then, (3.53) becomes 𝑈̇=1−𝑆^∗𝑆𝛽𝑁^∗𝑆^∗𝐼^∗₁+𝛽𝑁^∗𝑆^∗𝐼^∗₂+𝜇𝑆^∗−𝛽(𝑁)𝑆𝐼₁−𝛽(𝑁)𝑆𝐼₂+𝐸−𝜇𝑆1−^∗1𝐸₁𝛽(𝑁)𝑆𝐼₁+𝑝𝑟₂𝐼₁−𝛽𝑁^∗𝑆^∗𝐼^∗₁𝐸₁𝐸^∗₁−𝑝𝑟₂𝐼^∗₁𝐸₁𝐸^∗₁𝐼+𝐶1−^∗1𝐼₁𝑘₁1−𝑟₁𝐸₁−𝑘₁1−𝑟₁𝐸^∗₁𝐼₁𝐼^∗₁+𝐸1−^∗2𝐸₂𝛽(𝑁)𝑆𝐼₂+1−𝑝𝑟₂𝐼₁−𝛽𝑁^∗𝑆^∗𝐼^∗₂𝐸₂𝐸^∗₂−1−𝑝𝑟₂𝐼^∗₁𝐸₂𝐸^∗₂𝐼+𝐷1−^∗2𝐼₂𝑘₂𝐸₂−𝑘₂𝐸^∗₂𝐼₂𝐼^∗₂𝜇=−𝑆−𝑆^∗2𝑆+𝛽𝑁^∗𝑆^∗𝐼^∗₁𝛽(𝑁)1−𝑆𝐼₁𝛽(𝑁^∗)𝑆^∗𝐼^∗₁𝑆1−^∗𝑆+𝛽𝑁^∗𝑆^∗𝐼^∗₂𝛽(𝑁)1−𝑆𝐼₂𝛽(𝑁^∗)𝑆^∗𝐼^∗₂𝑆1−^∗𝑆+𝐸1−^∗1𝐸₁𝛽𝑁^∗𝑆^∗𝐼^∗₁𝛽(𝑁)𝑆𝐼₁𝛽(𝑁^∗)𝑆^∗𝐼^∗₁−𝐸₁𝐸^∗₁+𝑝𝑟₂𝐼^∗₁𝐼₁𝐼^∗₁−𝐸₁𝐸^∗₁𝐼+𝐶1−^∗1𝐼₁𝑘₁1−𝑟₁𝐸^∗₁𝐸₁𝐸^∗₁−𝐼₁𝐼^∗₁+𝐸1−^∗2𝐸₂𝛽𝑁^∗𝑆^∗𝐼^∗₂𝛽(𝑁)𝑆𝐼₂𝛽(𝑁^∗)𝑆^∗𝐼^∗₂−𝐸₂𝐸^∗₂+1−𝑝𝑟₂𝐼^∗₁𝐼₁𝐼^∗₁−𝐸₂𝐸^∗₂𝐼+𝐷1−^∗2𝐼₂𝑘₂𝐸^∗₂𝐸₂𝐸^∗₂−𝐼₂𝐼^∗₂. (3.55) Now, using (3.37), we have 𝑘₁1−𝑟₁𝐸^∗₁=𝑟₂+𝜇+𝑑₁𝐼^∗₁,𝑘₂𝐸^∗₂=𝜇+𝑑₂𝐼^∗₂. (3.56) Then, (3.55) can be rewritten as follows: 𝑈̇=−𝜇𝑆−𝑆^∗2𝑆+𝛽𝑁^∗𝑆^∗𝐼^∗₁𝛽(𝑁)1−𝑆𝐼₁𝛽(𝑁^∗)𝑆^∗𝐼^∗₁𝑆1−^∗𝑆+𝐸1−^∗1𝐸₁𝛽(𝑁)𝑆𝐼₁𝛽(𝑁^∗)𝑆^∗𝐼^∗₁−𝐸₁𝐸^∗₁+𝑝𝑟₂𝐼^∗₁𝐸1−^∗1𝐸₁𝐼₁𝐼^∗₁−𝐸₁𝐸^∗₁+𝐶𝜇+𝑟₂+𝑑₁𝑝𝑟₂𝐼1−^∗1𝐼₁𝐸₁𝐸^∗₁−𝐼₁𝐼^∗₁+𝛽𝑁^∗𝑆^∗𝐼^∗₂𝛽(𝑁)1−𝑆𝐼₂𝛽(𝑁^∗)𝑆^∗𝐼^∗₂𝑆1−^∗𝑆+𝛽(𝑁)𝑆𝐼₂𝛽(𝑁^∗)𝑆^∗𝐼^∗₂−𝐸₂𝐸^∗₂𝐸1−^∗2𝐸₂+1−𝑝𝑟₂𝐼^∗₁𝐼₁𝐼^∗₁−𝐸₂𝐸^∗₂𝐸1−^∗2𝐸₂+𝐷𝜇+𝑑₂𝐼^∗₂1−𝑝𝑟₂𝐼^∗₁𝐼1−^∗2𝐼₂𝐸₂𝐸^∗₂−𝐼₂𝐼^∗₂. (3.57) Now, let (𝑥,𝑤₁,𝑧₁,𝑤₂,𝑧₂)=(𝑆/𝑆^∗,𝐸₁/𝐸^∗₁,𝐼₁/𝐼^∗₁,𝐸₂/𝐸^∗₂,𝐼₂/𝐼^∗₂) and 𝑔(𝑁)=𝛽(𝑁)/𝛽(𝑁^∗), then we get 𝑈̇=−𝜇𝑆−𝑆^∗2𝑆+𝛽𝑁^∗𝑆^∗𝐼^∗₁1−𝑔(𝑁)𝑥𝑧₁11−𝑥+11−𝑤₁𝑔(𝑁)𝑥𝑧₁−𝑤₁+𝑝𝑟₂𝐼^∗₁11−𝑤₁𝑧₁−𝑤₁+𝐶𝜇+𝑟₂+𝑑₁𝑝𝑟₂11−𝑧₁𝑤₁−𝑧₁+𝛽𝑁^∗𝑆^∗𝐼^∗₂1−𝑔(𝑁)𝑥𝑧₂11−𝑥+𝑔(𝑁)𝑥𝑧₂−𝑤₂11−𝑤₂+1−𝑝𝑟₂𝐼^∗₁𝑧₁−𝑤₂11−𝑤₂+𝐷𝜇+𝑑₂𝐼^∗₂1−𝑝𝑟₂𝐼^∗₁11−𝑧₂𝑤₂−𝑧₂,𝜇=−𝑆−𝑆^∗2𝑆+𝑓𝑥,𝑤₁,𝑧₁,𝑤₂,𝑧₂, (3.58) where 𝑓𝑥,𝑤₁,𝑧₁,𝑤₂,𝑧₂=𝛽𝑁^∗𝑆^∗𝐼^∗₁𝑓₁𝑥,𝑤₁,𝑧₁+𝑝𝑟₂𝐼^∗₁𝑓₂𝑤₁,𝑧₁+𝛽𝑁^∗𝑆^∗𝐼^∗₂𝑓₃𝑥,𝑤₂,𝑧₂+1−𝑝𝑟₂𝐼^∗₁𝑓₄𝑧₁,𝑤₂,𝑧₂,𝑓₁𝑥,𝑤₁,𝑧₁=1−𝑔(𝑁)𝑥𝑧₁11−𝑥+11−𝑤₁𝑔(𝑁)𝑥𝑧₁−𝑤₁,𝑓₂𝑤₁,𝑧₁=11−𝑤₁𝑧₁−𝑤₁+𝐶𝜇+𝑟₂+𝑑₁𝑝𝑟₂11−𝑧₁𝑤₁−𝑧₁,𝑓₃𝑥,𝑤₂,𝑧₂=1−𝑔(𝑁)𝑥𝑧₂11−𝑥+𝑔(𝑁)𝑥𝑧₂−𝑤₂11−𝑤₂,𝑓₄𝑧₁,𝑤₂,𝑧₂=𝑧₁−𝑤₂11−𝑤₂+𝐷𝜇+𝑑₂𝐼^∗₂1−𝑝𝑟₂𝐼^∗₁11−𝑧₂𝑤₂−𝑧₂. (3.59) We can choose 𝐶=𝑝𝑟₂𝜇+𝑟₂+𝑑₁𝑟,𝐷=1−𝑝₂𝐼^∗₁𝜇+𝑑₂𝐼^∗₂. (3.60) For 𝑥=𝑤₁=𝑤₂, we have 𝑓₁𝑥,𝑤₁,𝑧₁=1−𝑔(𝑁)𝑥𝑧₁11−𝑥+11−𝑤₁𝑔(𝑁)𝑥𝑧₁−𝑤₁=𝑔(𝑁)𝑥𝑧₁1𝑥1−1+1−𝑤₁1+1−𝑥−𝑤₁11−𝑤₁1=1−𝑥−𝑤₁1+1=2−𝑥−𝑥, (3.61) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑥=1𝑓₂𝑤₁,𝑧₁=11−𝑤₁𝑧₁−𝑤₁+11−𝑧₁𝑤₁−𝑧₁=𝑧₁−𝑤₁11−𝑤₁1−1+𝑧₁=𝑧₁−𝑤₁1𝑧₁−1𝑤₁𝑤=2−¹𝑧₁−𝑧₁𝑤₁, (3.62) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑤₁=𝑧₁𝑓₃𝑥,𝑤₂,𝑧₂=1−𝑔(𝑁)𝑥𝑧₂11−𝑥+𝑔(𝑁)𝑥𝑧₂−𝑤₂11−𝑤₂=𝑔(𝑁)𝑥𝑧₂11−𝑤₂+1𝑥1−1+1−𝑥−𝑤₂1+1=1−𝑥−𝑤₂1+1=2−𝑥−𝑥, (3.63) which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if 𝑥=1.
For 𝑤₂≥1, 𝑧₂≥𝑧₁ or 𝑤₂≤1, 𝑧₂≤𝑧₁, we have𝑓₄𝑧₁,𝑤₂,𝑧₂=𝑧₁−𝑤₂11−𝑤₂+11−𝑧₂𝑤₂−𝑧₂≤𝑧₂−𝑤₂11−𝑤₂1−1+𝑧₂=𝑧₂−𝑤₂1𝑧₂−1𝑤₂𝑤=2−²𝑧₂−𝑧₂𝑤₂, (3.64) which is less than or equal to zero by the arithmetic-geometric mean inequality, with equality if and only if 𝑤₂=𝑧₂.
Thus, 𝑈̇(𝑆,𝐸₁,𝐼₁,𝐸₂,𝐼₂) is less or equal to zero with equality only if 𝑥=1,𝑤₁=𝑧₁,𝑤₂=𝑧₂. Since 𝑥=𝑤₁=𝑤₂, the largest invariant set of system (2.2) on the set {(𝑆,𝐸₁,𝐼₁,𝐸₂,𝐼₂)∈ℝ⁵₊, 𝑈̇(𝑆,𝐸₁,𝐼₁,𝐸₂,𝐼₂)=0} is the endemic equilibrium point 𝑃₂. By LaSalle's invariance principle [31, 32], we can conclude that the interior equilibrium 𝑃₂ 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 𝑑₁=𝑑₂=0, 𝛽₁=0.001, 𝛽₂=0.0001 (so that 𝑅₁<1, 𝑅₂<1). For 𝛽₁(𝑁)=𝛽₁, 𝛽₂(𝑁)=𝛽₂, 𝛽₁(𝑁)=𝛽₁/𝑁, 𝛽₂(𝑁)=𝛽₂/𝑁 and 𝛽₁(𝑁)=𝛽₁𝐶(𝑁), 𝛽₂(𝑁)=𝛽₂𝐶(𝑁), 𝐶(𝑁)=𝑏𝑁/(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.001, 𝛽2=0.0001, 𝑑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.001, 𝛽2=0.0001, 𝑑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.001, 𝛽2=0.0001, 𝑑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 𝑑₁>0, 𝑑₂>0, we take 𝑑₁=0.02, 𝑑₂=0.07. Firstly, we consider system (2.2) with 𝛽₁(𝑁)=𝛽₁, 𝛽₂(𝑁)=𝛽₂ and choose 𝛽₁=𝛽₂=0.001, so that 𝑅₁<1, 𝑅₂<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 𝛽₁=0.001, 𝛽₂=1, so that 𝑅₁<1, 𝑅₂>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, 𝛽₂=0.25, so that 𝑅₁>𝐻, 𝑅₂<𝑓(𝑅₁), 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.001, 𝑑1=0.02, 𝑑2=0.07, 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.001, 𝛽2=1, 𝑑1=0.02, 𝑑2=0.07, 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.25, 𝑑1=0.02, 𝑑2=0.07, so that 𝑅1>𝐻, 𝑅2<𝑓(𝑅1). It shows that the interior equilibrium 𝑃2 is stable.

Next, we need to think the case when 𝛽₁(𝑁)=𝛽₁/𝑁, 𝛽₂(𝑁)=𝛽₂/𝑁. Figures 8, 9 and 10 show the trajectories plot and its plane figure of system (2.2) when 𝛽₁=𝛽₂=0.001 (so that 𝑅₁<1, 𝑅₂<1), 𝛽₁=0.001, 𝛽₂=1 (so that 𝑅₁<1, 𝑅₂>1), and 𝛽₁=1, 𝛽₂=0.25 (so that 𝑅₁>𝐻, 𝑅₂<𝑓(𝑅₁)), respectively. As the previous cases, we can find that the trajectories of system (2.2) converge to the disease-free equilibrium when 𝑅₁<1, 𝑅₂<1 (see Figure 8), the drug-resistant-TB-strain-only equilibrium is global asymptotically stable when 𝑅₁<1, 𝑅₂>1 (see Figure 9), while the trajectories of system (2.2) converge to the unique interior equilibrium when 𝑅₁>1, 𝑅₂<𝑓(𝑅₁) (see Figure 10).

Figure 8 

Trajectories of system (2.2) when 𝛽1(𝑁)=𝛽1/𝑁, 𝛽2(𝑁)=𝛽2/𝑁, 𝛽1=𝛽2=0.001, 𝑑1=0.02, 𝑑2=0.07, 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.001, 𝛽2=1, 𝑑1=0.02, 𝑑2=0.07, 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.25, 𝑑1=0.02, 𝑑2=0.07, so that 𝑅1>𝐻, 𝑅2<𝑓(𝑅1). It shows that the interior equilibrium 𝑃2 is stable.

Eventually, we assume that 𝛽₁(𝑁)=𝛽₁𝐶(𝑁), 𝛽₂(𝑁)=𝛽₂𝐶(𝑁), 𝐶(𝑁)=𝑏𝑁/(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 𝛽₁=𝛽₂=0.001 (so that 𝑅₁<1, 𝑅₂<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 𝛽₁=0.001, 𝛽₂=1 (so that 𝑅₁<1, 𝑅₂>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, 𝛽₂=0.25 (so that 𝑅₁>𝐻, 𝑅₂<𝑓(𝑅₁)). 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.001, 𝑑1=0.02, 𝑑2=0.07, 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.001, 𝛽2=1, 𝑑1=0.02, 𝑑2=0.07, 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.25, 𝑑1=0.02, 𝑑2=0.07, 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 𝑅₁, 𝑅₂ provide the threshold conditions which determine the competitive outcomes of the two strains. When 𝑅₁<1, 𝑅₂<1, the disease-free equilibrium 𝑃₀ is globally stable and the disease will die out. When 𝑅₁<1, 𝑅₂>1, a unique boundary equilibrium 𝑃₁ exists and is globally stable. While when 𝑅₁>𝐻, 𝑅₂<𝑓(𝑅₁), there exists a unique interior equilibrium 𝑃₂ 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.