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

Authors: Hai-Feng Huo; Shuai-Jun Dang; Yu-Ning Li
Journal: Abstract and Applied Analysis
 (2010)

---

## Abstract

A two-strain tuberculosis model with general contact rate which allows tuberculosis patients with the drug-sensitiveMycobacterium 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
R1 and R2, 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.


---

## Body

## 1. Introduction

Tuberculosis (TB), an infectious disease caused byMycobacterium 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 pathogenicM. 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 βC(N)SI/N, where S and I are the numbers of susceptible and infective at time t, respectively, β is the probability of transmitting the infection between two individuals taking part in a contact per unit time, C(N) is the probability for an individual to take part in a contact and N 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 C(N)=N, the corresponding incidence is the bilinear form βNSI/N=βSI. When C(N)=λ, the corresponding incidence λSI/N is called standard form. When the total population size N is not quite large, since the number of contacts made by an individual per unit time should increase as the total population size N 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 N 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 N 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)C(N)=bN1+bN+1+2bN,
where C(N) is nondecreasing and C(N)/N 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, R1 and R2, 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 R1<1, R2<1, the unique drug-resistant-TB-strain-only equilibrium exists and is globally asymptotically stable when R1<1<R2. As for the unique interior equilibrium, existence, uniqueness and global stability can be proved under certain conditions.The paper is organized as follows. In Section2, 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(S), those exposed to drug-sensitive TB (E1), individuals with symptoms of TB and drug-sensitive (I1), those exposed to drug-resistant TB (E2) and those displaying symptoms of TB and drug-resistant (I2), the total number of population at time t is given by(2.1)N(t)=S(t)+E1(t)+I1(t)+E2(t)+I2(t).
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 βi(N)Ii, i=1,2, where i=1 and i=2 represent rates of infection by drug-sensitive strain and drug-resistant strain, and the transmission coefficient βi(N) is a nonnegative C2 function of the total population N. Susceptible are infected with drug-sensitive strain and drug-resistant strain entering classes E1 and E2, respectively. For simplicity, we only account for treatment of sensitive strain, we assume that chemoprophylaxis of individuals in E1 reduces their reactivation at a constant rate r1 and that the initiation of therapeutics immediately removes individuals in I1 from active status and places them into exposed class at a constant rate r2. The time before individuals in E1 who does not received effective chemoprophylaxis become infectious is assumed to satisfy an exponential distribution, with mean waiting time 1/k1. Thus, individuals leave the class E1 to the class I1 at a constant rate k1(1-r1). The time before individuals in E2 become infectious is also assumed to satisfy an exponential distribution, with mean waiting time 1/k2. Thus, individuals leave the class E2 to the class I2 at a constant rate k2. Individuals sick with drug-sensitive strain receive treatment at rate r2 and a proportion pr2 respond to treatment and move into exposed class E1, in the remaining proportion (1-p)r2, inappropriate treatment results in the development of drug-resistant strain and the individuals move into class E2. Individuals in the infectious class I1 and I2 suffer additional disease-induced death at rate d1 and d2, where d2>d1. Then the following two strains tuberculosis model is formulated(2.2)S′=Λ-β1(N)SI1-β2(N)SI2-μS,E1′=β1(N)SI1+pr2I1-k1(1-r1)E1-μE1,I1′=k1(1-r1)E1-(r2+μ+d1)I1,E2′=β2(N)SI2+(1-p)r2I1-(μ+k2)E2,I2′=k2E2-(μ+d2)I2,
where Λ, k1, k2, r1, and r2 are assumed to be positive and p∈[0,1]. It is natural to assume that the transmission coefficient β1(N),β2(N) satisfies the following conditions:(2.3)β1(N)>0,β1′(N)≤0,(Nβ1(N))′≥0,β2(N)>0,β2′(N)≤0,(Nβ2(N))′≥0.Figure 1
Transfer diagram for the dynamics of tuberculosis with general contact rate.Obviously, we can see thatβ1(N)=β1/N, β2(N)=β2/N corresponds to the standard incidence rate, that β1(N)=β1, β2(N)=β2 corresponds to the mass action incidence rate, and that β1(N)=β1C(N), β2(N)=β2C(N) corresponds to the saturating contact rate, where(2.4)C(N)=bN1+bN+1+2bN.


### 2.1. Positivity and Boundedness

System (2.2) describes human population and therefore it is necessary to prove that all the variables S, E1, I1, E2, I2 are nonnegative for all time. Solutions of system (2.2) with positive initial data remains positive for all time t≥0 and are bounded.Theorem 2.1.
IfS(0)≥0, E1(0)≥0, I1(0)≥0, E2(0)≥0, I2(0)≥0, the solutions S(t), E1(t), I1(t), E2(t), I2(t) of system (2.2) are positive for t≥0. For system (2.2), the region Γ is positively invariant and all solutions starting in Γ approach, enter, or stay in Γ, where
(2.5)Γ={(S(t),E1(t),I1(t),E2(t),I2(t))∈R+5:N(t)≤Λμ}.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 t1 such that
(2.6)S(t1)=0,S′(t1)<0,E1(t)≥0,I1(t)≥0,E2(t)≥0,I2(t)≥0,for0<t<t1,
there exists a t2(2.7)E1(t2)=0,E1′(t2)<0,S(t)≥0,I1(t)≥0,E2(t)≥0,I2(t)≥0,for0<t<t2,
there exists a t3(2.8)I1(t3)=0,I1′(t3)<0,S(t)≥0,E1(t)≥0,E2(t)≥0,I2(t)≥0,for0<t<t3,
there exists a t4(2.9)E2(t4)=0,E2′(t4)<0,S(t)≥0,E1(t)≥0,I1(t)≥0,I2(t)≥0,for0<t<t4,
or there exists a t5(2.10)I2(t5)=0,I2′(t5)<0,S(t)≥0,E1(t)≥0,I1(t)≥0,E2(t)≥0,for0<t<t5.
In the first case, we have
(2.11)S′(t1)=Λ>0,
which is a contradiction meaning that S(t)≥0, t≥0. In the second case, we have
(2.12)E1′(t2)=I1(t2)[β1(N)S(t2)+pr2]≥0,
which is a contradiction meaning that E1(t)≥0, t≥0. In the third case, we have
(2.13)I1′(t3)=k1(1-r1)E1(t3)≥0,
which is a contradiction meaning that I1(t)≥0, t≥0. In the fourth case, we have
(2.14)E2′(t4)=β2(N)S(t4)I2(t4)+(1-p)r2I1(t4)≥0,
which is a contradiction meaning that E2(t)≥0, t≥0. In the fifth case, we have
(2.15)I2′(t5)=k2E2(t5)≥0,
which is a contradiction meaning that I2(t)≥0, t≥0.
Thus, in all cases,S(t), E1(t), I1(t), E2(t), I2(t) remain positive for t≥0.
Let(S(t),E1(t),I1(t),E2(t),I2(t))∈ℝ+5 be any solution with nonnegative initial condition, adding all equations in system (2.2) gives (2.16)N′(t)=Λ-μN-d1I1-d2I2,≤Λ-μN.
It follows that
(2.17)0≤N(t)≤Λμ+N(0)e-μt,
where N(0) represents initial values of the total population. Thus 0≤N(t)≤Λ/μ, as t→∞. Therefore all feasible solutions of system (2.2) enter the region
(2.18)Γ={(S(t),E1(t),I1(t),E2(t),I2(t))∈R+5:N(t)≤Λμ}.
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 N(t) is bounded and all the solutions starting in Γ approach, enter or stay in Γ.


## 2.1. Positivity and Boundedness

System (2.2) describes human population and therefore it is necessary to prove that all the variables S, E1, I1, E2, I2 are nonnegative for all time. Solutions of system (2.2) with positive initial data remains positive for all time t≥0 and are bounded.Theorem 2.1.
IfS(0)≥0, E1(0)≥0, I1(0)≥0, E2(0)≥0, I2(0)≥0, the solutions S(t), E1(t), I1(t), E2(t), I2(t) of system (2.2) are positive for t≥0. For system (2.2), the region Γ is positively invariant and all solutions starting in Γ approach, enter, or stay in Γ, where
(2.5)Γ={(S(t),E1(t),I1(t),E2(t),I2(t))∈R+5:N(t)≤Λμ}.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 t1 such that
(2.6)S(t1)=0,S′(t1)<0,E1(t)≥0,I1(t)≥0,E2(t)≥0,I2(t)≥0,for0<t<t1,
there exists a t2(2.7)E1(t2)=0,E1′(t2)<0,S(t)≥0,I1(t)≥0,E2(t)≥0,I2(t)≥0,for0<t<t2,
there exists a t3(2.8)I1(t3)=0,I1′(t3)<0,S(t)≥0,E1(t)≥0,E2(t)≥0,I2(t)≥0,for0<t<t3,
there exists a t4(2.9)E2(t4)=0,E2′(t4)<0,S(t)≥0,E1(t)≥0,I1(t)≥0,I2(t)≥0,for0<t<t4,
or there exists a t5(2.10)I2(t5)=0,I2′(t5)<0,S(t)≥0,E1(t)≥0,I1(t)≥0,E2(t)≥0,for0<t<t5.
In the first case, we have
(2.11)S′(t1)=Λ>0,
which is a contradiction meaning that S(t)≥0, t≥0. In the second case, we have
(2.12)E1′(t2)=I1(t2)[β1(N)S(t2)+pr2]≥0,
which is a contradiction meaning that E1(t)≥0, t≥0. In the third case, we have
(2.13)I1′(t3)=k1(1-r1)E1(t3)≥0,
which is a contradiction meaning that I1(t)≥0, t≥0. In the fourth case, we have
(2.14)E2′(t4)=β2(N)S(t4)I2(t4)+(1-p)r2I1(t4)≥0,
which is a contradiction meaning that E2(t)≥0, t≥0. In the fifth case, we have
(2.15)I2′(t5)=k2E2(t5)≥0,
which is a contradiction meaning that I2(t)≥0, t≥0.
Thus, in all cases,S(t), E1(t), I1(t), E2(t), I2(t) remain positive for t≥0.
Let(S(t),E1(t),I1(t),E2(t),I2(t))∈ℝ+5 be any solution with nonnegative initial condition, adding all equations in system (2.2) gives (2.16)N′(t)=Λ-μN-d1I1-d2I2,≤Λ-μN.
It follows that
(2.17)0≤N(t)≤Λμ+N(0)e-μt,
where N(0) represents initial values of the total population. Thus 0≤N(t)≤Λ/μ, as t→∞. Therefore all feasible solutions of system (2.2) enter the region
(2.18)Γ={(S(t),E1(t),I1(t),E2(t),I2(t))∈R+5:N(t)≤Λμ}.
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 N(t) is bounded and all the solutions starting in Γ approach, enter or stay in Γ.


## 3. Model Analysis

There are one disease-free equilibriumP0 and two possible endemic equilibria for system (2.2), the endemic equilibria include boundary equilibrium P1 (when only the drug-resistant TB-strain is present) and the interior equilibrium P2 (when both strains exist).


### 3.1. Global Stability of the Disease-Free Equilibrium

The disease-free-equilibrium is given as(3.1)P0=(S0,0,0,0,0),S0=Λμ.
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 P0, we use the results in van den Driessche and Watmough [27]. Reorder the components of P0 as E1=0, I1=0, E2=0, I2=0, S=Λ/μ. Set(3.2)F=[F1F2F3F4F5]=[β1(N)SI10β2(N)SI200],V=[k1(1-r1)E1+μE1-pr2I1r2I1+(μ+d1)I1-k1(1-r1)E1(μ+k2)E2-(1-p)r2I1(μ+d2)I2-k2E2μS+β1(N)SI1+β2(N)SI2-Λ].
The infected compartments are E1, I1, E2, and I2. Thus(3.3)F=[0β1(S0)S0000000000β1(S0)S00000],V=[μ+k1(1-r1)-pr200-k1(1-r1)μ+d1+r2000-(1-p)r2μ+k2000-k2μ+d2].
The dominant eigenvalues of FV-1 are given by (3.4)R1=k1(1-r1)β1(S0)S0(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1),R2=k2β2(S0)S0(μ+k2)(μ+d2),
where R1 and R2 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 FV-1 is(3.5)ρ(FV-1)=max⁡{R1,R2}.
If R1<1, and R2<1, then ρ(FV-1)<1. By Theorem 3.2 in van den Driessche and Watmough [27], we know that the disease-free steady state P0 is locally asymptotically stable, P0 is unstable if R1>1 or R2>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,d1=0, d2=0. Consequently, the total population size N(t) satisfies the equation dN/dt=Λ-μN and N(t)→Λ/μ as t→∞. 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 N=S0=Λ/μ: (3.6)E1′=β1(S0)I1(S0-E1-I1-E2-I2)+pr2I1-(μ+k1(1-r1))E1,I1′=k1(1-r1)E1-(r2+μ)I1,E2′=β2(S0)I2(S0-E1-I1-E2-I2)+(1-p)r2I1-(μ+k2)E2,I2′=k2E2-μI2.
Notice that the S equation is eliminated from (3.6) and the variable S is replaced by S0-E1-I1-E2-I2. By introducing the fractions(3.7)x1=E1N,x2=I1N,y1=E2N,y2=I2N,
we obtain the equivalent limiting system from (3.6) (3.8)x1′=S0β1(S0)(1-x1-x2-y1-y2)x2+pr2x2-(μ+k1(1-r1))x1,x2′=k1(1-r1)x1-(μ+r2)x2,y1′=S0β2(S0)(1-x1-x2-y1-y2)y2+(1-p)r2x2-(μ+k2)y1,y2′=k2y1-μy2.
Obviously(3.9)0≤x1+x2+y1+y2≤1,
for all time t≥0.For a bounded real-valued functionf on [0,∞), we define(3.10)f∞=lim⁡inf⁡t→∞f(t),f∞=limsup⁡t→∞f(t).Lemma 3.1 (see [30]).
Letf:[0,∞)→∞ be bounded and continuously differentiable. Then there are sequences tn,sn→∞ with the following properties:
(3.11)f(tn)⟶f∞,f′(tn)⟶0,f(sn)⟶f∞,f′(sn)⟶0,
for n→∞.Theorem 3.2.
In the absence of disease-induced death rate, the disease-free equilibriumP0 of system (2.2) is globally asymptotically stable if R1<1, R2<1.Proof.
By Lemma3.1 and the x2′, y2′ equations in (3.8) we have
(3.12)x2∞≤k1(1-r1)μ+r2x1∞,y2∞≤k2μy1∞.
Using thex1′ equation in (3.8) and choosing tn→∞ such that(3.13)x1(tn)⟶x1∞,x1′(tn)⟶0,t⟶∞,
we get
(3.14)0≤S0β1(S0)(1-x1-x2-y1-y2)∞x2∞-(μ+k1(1-r1))x1∞≤S0β1(S0)x2∞-(μ+k1(1-r1))x1∞≤S0β1(S0)k1(1-r1)μ+r2x1∞-(μ+k1(1-r1))x1∞≤(μ+k1(1-r1))x1∞[k1(1-r1)S0β1(S0)(μ+r2)(μ+k1(1-r1))-1]≤(μ+k1(1-r1))x1∞[k1(1-r1)S0β1(S0)(μ+r2)(μ+k1(1-r1))-k1pr2(1-r1)-1]=(μ+k1(1-r1))x1∞(R1-1).
It is shown that x1∞≤0 since R1<1, as x1∞≥0, we have that x1∞=0. The inequalities in (3.12) also imply that x2∞=0.
Using they1′ equation in (3.8) and choosing sn→∞ such that(3.15)y1(sn)⟶y1∞,y1′(sn)⟶0,t→∞,
we get
(3.16)0≤S0β2(S0)(1-x1-x2-y1-y2)∞y2∞-(μ+k2)y1∞≤S0β2(S0)y2∞-(μ+k2)y1∞≤S0β2(S0)k2μy1∞-(μ+k2)y1∞≤(μ+k2)y1∞[k2S0β2(S0)μ(μ+k2)-1]=(μ+k2)y1∞(R2-1).
It is shown that y1∞≤0 since R2<1, as y1∞≥0, we have that y1∞=0. The inequalities in (3.12) also imply that y2∞=0.
Hence(3.17)lim⁡t→∞x1(t)=lim⁡t→∞x2(t)=lim⁡t→∞y1(t)=lim⁡t→∞y2(t)=0.
Therefore, the disease-free equilibrium P0 is globally asymptotically stable.Next, we consider that there is disease-induced death rate, which means thatd1>0, d2>0. Then, we can get the following theorem.Theorem 3.3.
Ford1>0, d2>0, the disease-free equilibrium P0 of system (2.2) is globally asymptotically stable if R1<1-k2r2(1-p)/((μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)), R2<1.Proof.
Let us consider the following Lyapunov function:(3.18)V(E1,I1,E2,I2)=k1(1-r1)E1+(μ+k1(1-r1))I1+k2E2+(μ+k2)I2.
Its time derivative along the solutions of system (2.2) satisfies
(3.19)V̇=k1(1-r1)Ė1+(μ+k1(1-r1))İ1+k2Ė2+(μ+k2)İ2=k1(1-r1)[β1(N)SI1+pr2I1-(μ+k1(1-r1))E1]+(μ+k1(1-r1))[k1(1-r1)E1-(μ+d1+r2)I1]+k2[β2(N)SI2+(1-p)r2I1-(μ+k2)E2]+(μ+k2)[k2E2-(μ+d2)I2],=[k1(1-r1)β1(N)S+k1(1-r1)pr2-(μ+d1+r2)(μ+k1(1-r1))]I1+[k2β2(N)S-(μ+d2)(μ+k2)]I2+k2r2(1-p)I1.
Now, using (2.3), one has β1(N)S≤β1(S0)S0, β2(N)S≤β2(S0)S0, then
(3.20)V̇≤[(μ+d1+r2)(μ+k1(1-r1))-k1(1-r1)pr2](R1-1)I1+k2r2(1-p)I1+(μ+d2)(μ+k2)(R2-1)I2,=[[(μ+d1+r2)(μ+k1(1-r1))-k1(1-r1)pr2](R1-1)+k2r2(1-p)]I1+(μ+d2)(μ+k2)(R2-1)I2.
Since
(3.21)(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)=μ(μ+d1+r2)+k1(1-r1)(μ+d1)+k1r2(1-r1)-k1pr2(1-r1)=μ(μ+d1+r2)+k1(1-r1)(μ+d1)+k1r2(1-r1)(1-p)≥0,
and R1<1-k2r2(1-p)/((μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1))≤1, we have
(3.22)[(μ+d1+r2)(μ+k1(1-r1))-k1(1-r1)pr2](R1-1)+k2r2(1-p)≤0.
Thus, V̇≤0 if R1<1-k2r2(1-p)/((μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)), R2<1. It is obvious that V̇=0 if and only if I1=0, I2=0. Then, the largest invariant set of system (2.2) on the set {(S,E1,I1,E2,I2)∈ℝ+5, V̇(E1,I1,E2,I2)=0} is the disease-free equilibrium P0. Therefore, it follows from LaSalle's invariance principle [31, 32] that P0 is globally stable if R1<1-k2r2(1-p)/((μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)), R2<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 R1 and R2.LetP1=(Ŝ,0,0,Ê2,Î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 (3.23)Λ-β2(N̂)ŜÎ2-μŜ=0,β2(N̂)ŜÎ2-(μ+k2)Ê2=0,k2Ê2-(μ+d2)Î2=0,
adding the above three equations, we have(3.24)Λ-μN̂-d2Î2=0.Using (3.24), the first and second equation of (3.23), we can easily express Ŝ, Ê2 and Î2 in terms of N̂ in the form:
(3.25)Ŝ=Λd2μd2+β2(N̂)(Λ-μN̂),Ê2=Λβ2(N̂)(Λ-μN̂)(μ+k2)[μd2+β2(N̂)(Λ-μN̂)],Î2=Λ-μN̂d2.SubstitutingŜ, Ê2, Î2 in the third equations of (3.23) yields(3.26)(Λ-μN̂)F(N̂)=0,
where(3.27)F(N̂)=Λk2d2β2(N̂)-(μ+k2)(μ+d2)[μd2+β2(N̂)(Λ-μN̂)].Clearly,Λ-μN̂=0 is a fixed point of (3.23), which corresponds to the disease-free equilibrium P0. Since N̂∈[0,S0], we get(3.28)F(0)=Λk2d2β2(0)-(μ+k2)(μ+d2)[μd2+Λβ2(0)]=Λk2d2β2(0)-[μ2+μ(d2+k2)+d2k2][μd2+Λβ2(0)]=Λk2d2β2(0)-[μ2+μ(d2+k2)][μd2+Λβ2(0)]-Λd2k2β2(0)-μk2d22=-[μ2+μ(k2+d2)][μd2+Λβ2(0)]-μk2d22,F(S0)=Λk2d2β2(S0)-μd2(μ+k2)(μ+d2)=μ[k2d2β2(S0)S0-d2(μ+k2)(μ+d2)]=μd2(μ+k2)(μ+d2)[k2β2(S0)S0(μ+k2)(μ+d2)-1]=μd2(μ+k2)(μ+d2)(R2-1).It appears thatF(0)<0 and F(S0)>0 when R2>1. The existence of fixed point follows from the intermediate value-theorem. Now, F(N̂) is monotone increasing, so that F(N̂)=0 has only one positive root in the interval [0,S0].Thus, we have established the following resultLemma 3.4.
WhenR2>1, the system (2.2) has a unique drug-resistant-TB-strain-only equilibrium P1=(Ŝ,0,0,Ê2,Î2) with Ŝ, Ê2, and Î2 all nonnegative.


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

Theorem 3.5.
IfR1<1<R2, the drug-resistant-TB-strain-only equilibrium P1 exists and is locally asymptotically stable.Proof.
Reorder the equilibriumP1 as (0,0,Ê2,Î2,Ŝ). Similarly as in the proof of locally stability of the disease-free steady state P0, we have
(3.29)F1=[0β1(N̂)Ŝ00],V1=[μ+k1(1-r1)-pr2-k1(1-r1)μ+d1+r2],F1V1-1=A[k1(1-r1)β1(N̂)Ŝ(μ+k1(1-r1))β1(N̂)Ŝ00],
where
(3.30)A=1(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1).
The spectral radius of F1V1-1 is given by
(3.31)ρ(F1V1-1)=k1(1-r1)β1(N̂)Ŝ(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1),≤k1(1-r1)β1(N̂)N̂(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1).
Since (β1(N)N)′≥0, and N̂≤S0, then
(3.32)β1(N̂)N̂≤β1(S0)S0,
Thus
(3.33)ρ(F1V1-1)≤k1(1-r1)β1(S0)S0(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)=R1.
As R1<1, we have ρ(F1V1-1)<1, which implies that the drug-resistant-TB-strain-only equilibrium P1 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 equilibriumP1 is global asymptotically stable when R1<1<R2.Proof.
IfR1<1<R2, from Theorem 3.5, P1 is locally asymptotically stable. In the following, we only need to prove that the P1 is a global attractor.
There is no disease-induced death rate, which meansd1=d2=0, using a similar argument as in the proof of Theorem 3.2, we can show that if R1<1,(3.34)lim⁡t→∞x1(t)=lim⁡t→∞x2(t)=0.
Then the equivalent limiting system of (3.6) is
(3.35)E2′=β2(S0)I2(S0-E2-I2)-(μ+k2)E2,I2′=k2E2-μI2.
Let(P,Q) be the vector field defined by system (3.35). Then for the Dulac function D(E2,I2)=1/E2I2, there holds(3.36)∂E2(DP)+∂I2(DQ)=-β2(S0)(S0-I2)E22-k2I22<0.
From Dulac's criterion, we can see equivalent limiting system (3.35) does not have a limit cycle. Therefore, the local stability of P1 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(W1)
limN→0(1/β1(N))=+∞,(W2)
β1(N)/β2(N)=β1/β2,then, existence, uniqueness and global stability of interior equilibrium can be obtained. LetP2=(S*,E1*,I1*,E2*,I2*) be interior equilibrium of system (2.2). Then, the interior equilibrium (steady state with I1>0, I2>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 (3.37)Λ-β1(N*)S*I1*-β2(N*)S*I2*-μS*=0,β1(N*)S*I1*+pr2I1*-k1(1-r1)E1*-μE1*=0,k1(1-r1)E1*-(r2+μ+d1)I1*=0,β2(N*)S*I2*+(1-p)r2I1*-(μ+k2)E2*=0,k2E2*-(μ+d2)I2*=0,
adding the above five equations, we have(3.38)Λ-μN*-d1I1*-d2I2*=0.Using (3.37) and (3.38), we get the following equivalent equations set (3.39)Λ-β1(N*)S*I1*-β2(N*)S*I2*-μS*=0,β1(N*)S*I1*+pr2I1*-k1(1-r1)E1*-μE1*=0,k1(1-r1)E1*-(r2+μ+d1)I1*=0,k2E2*-(μ+d2)I2*=0,Λ-μN*-d1I1*-d2I2*=0.
From (3.39), we have (3.40)S*=β(S0)S0R1β1(N*),I1*=d2(Λ-μS*)-S*β2(N*)(Λ-μN*)d2S*β1(N*)-d1S*β2(N*),E1*=μ+d1+r2k1(1-r1)I1*,I2*=(Λ-μN*)d2-d1d2I1*,E2*=(μ+d2)(Λ-μN*)k2d2-d1(μ+d2)k2d2I1*.Substituting (3.40) into β2(N*)S*I2*+(1-p)r2I1*-(μ+k2)E2*=0, N*satisfies the following equation:(3.41)G(N*)=0,
where (3.42)G(N*)=β1(S0)S0(Λ-μN*)R1[k2d2S0β1(S0)β2(N*)R1β1(N*)-d2μ1-(1-p)k2d2r2β2(N*)β1(N*)]+d2Λ(1-β1(S0)β1(N*)R1)[d1μ1+(1-p)k2d2r2-d1k2S0β1(S0)β2(N*)R1β1(N*)],μ1=(μ+k2)(μ+d2).
Since β1(N)/β2(N)=β1/β2, limN→0(1/β1(N))=+∞, if R1>d1k2β2S0β1(S0)/β1[d1μ1+(1-p)k2d2r2], then(3.43)lim⁡N*→0G(N*)=-∞.
It means that G(N*)=0 has only one positive root in the interval [0,S0] if G(N*) is monotone increasing and G(S0)≥0, which is equivalent to (3.44)R1>H1+H2H3+H4,R2<[1+(1-p)d2k2r2d1μ1]R1,
where (3.45)H1=k2d2β2S0β1(S0)β12(N*),H2=d1d2k2β2S0β1(S0)β1′(N*),H3=β12(N*)[d2β1μ1+d2k2β2r2(1-p)],H4=d2β1′(N*)[d1β1μ1+d2k2β1r2(1-p)].
Let (3.46)H=max⁡{1,H1+H2H3+H4,d1k2β2S0β1(S0)β1[d1μ1+k2d2r2(1-p)]},f(R1)=[1+(1-p)d2k2r2d1μ1]R1.
When(3.47)R1>H,R2<f(R1),
it shows that G(S0)≥0, and G(N*) is monotone increasing, by the intermediate value-theorem, G(N*)=0 has only one positive root in the interval [0,S0], which signify that P2 is unique. Obviously, S*>0. In order to make certain that E1*>0, I1*>0, E2*>0, I2*>0, we obtain that E1*>0, I1*>0, E2*>0, I2*>0 if and only if(3.48)0<I1*<Λ-μN*d1.
Due to d2>d1>0, introducing (3.40) into (3.48), we acquire(3.49)B1<R1<B2,
where (3.50)B1=S0β1(S0)[β2(N*)(Λ-μN*)+μd2]Λd2β1(N*),B2=S0β1(S0)[β1(N*)(Λ-μN*)+μd1]Λd1β1(N*),
and B1<B2.Now, we can get the following conclusion.Lemma 3.7.
When conditions (W1), (W2), (3.47) and (3.49) are satisfied, the system (2.2) has a unique interior equilibrium P2=(S*,E1*,I1*,E2*,I2*) with S*, E1*, I1*, E2*, and I2* all nonnegative.Remark 3.8.
Forβ1(N)=β1 and β1(N)=β1/N, condition (W1) is not satisfied. However, at least when β1(N)=β1C(N), condition (W1) 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 equilibriumP2 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(N)=β2(N)=β(N), then, we get the following result.Theorem 3.9.
If conditions (W1), (W2), (3.47) and (3.49) are satisfied, interior equilibrium P2 of system (2.2) is is globally asymptotically stable when
(3.51)SS*=E1E1*=E2E2*≥1,I2I2*≥I1I1*orSS*=E1E1*=E2E2*≤1,I2I2*≤I1I1*.Proof.
Following [33], we consider the Lyapunov function:
(3.52)U(S,E1,I1,E2,I2)=(S-S*ln⁡S)+(E1-E1*ln⁡E1*)+C(I1-I1*ln⁡I1*)+(E2-E2*ln⁡E2*)+D(I2-I2*ln⁡I2*),
where C and D are positive constants to be determined later. Differentiating this function with respect to time along the solutions of system (2.2) yields
(3.53)U̇=(1-S*S)Ṡ+(1-E1*E1)Ė1+C(1-I1*I1)İ1+(1-E2*E2)Ė2+D(1-I2*I2)İ2=(1-S*S)[Λ-β(N)SI1-β(N)SI2-μS]+(1-E1*E1)[β(N)SI1+pr2I1-k1(1-r1)E1-μE1]+C(1-I1*I1)[k1(1-r1)E1-(r2+μ+d1)I1]+(1-E2*E2)[β(N)SI2+(1-p)r2I1-(μ+k2)E2]+D(1-I2*I2)[k2E2-(μ+d2)I2].
Considering (3.37), we can deduce that
(3.54)Λ=β(N*)S*I1*+β(N*)S*I2*+μS*,μ+k1(1-r1)=β(N*)S*I1*E1*+pr2I1*E1*,μ+d1+r2=k1(1-r1)E1*I1*,μ+k2=β(N*)S*I2*E2*+(1-p)r2I1*E2*,μ+d2=k2E2*I2*.
Then, (3.53) becomes
(3.55)U̇=(1-S*S)[β(N*)S*I1*+β(N*)S*I2*+μS*-β(N)SI1-β(N)SI2-μS]+(1-E1*E1)[β(N)SI1+pr2I1-β(N*)S*I1*E1E1*-pr2I1*E1E1*]+C(1-I1*I1)[k1(1-r1)E1-k1(1-r1)E1*I1I1*]+(1-E2*E2)[β(N)SI2+(1-p)r2I1-β(N*)S*I2*E2E2*-(1-p)r2I1*E2E2*]+D(1-I2*I2)[k2E2-k2E2*I2I2*]=-μ(S-S*)2S+β(N*)S*I1*(1-β(N)SI1β(N*)S*I1*)(1-S*S)+β(N*)S*I2*(1-β(N)SI2β(N*)S*I2*)(1-S*S)+(1-E1*E1)[β(N*)S*I1*(β(N)SI1β(N*)S*I1*-E1E1*)+pr2I1*(I1I1*-E1E1*)]+C(1-I1*I1)[k1(1-r1)E1*(E1E1*-I1I1*)]+(1-E2*E2)[β(N*)S*I2*(β(N)SI2β(N*)S*I2*-E2E2*)+(1-p)r2I1*(I1I1*-E2E2*)]+D(1-I2*I2)[k2E2*(E2E2*-I2I2*)].
Now, using (3.37), we have
(3.56)k1(1-r1)E1*=(r2+μ+d1)I1*,k2E2*=(μ+d2)I2*.
Then, (3.55) can be rewritten as follows:
(3.57)U̇=-μ(S-S*)2S+β(N*)S*I1*[(1-β(N)SI1β(N*)S*I1*)(1-S*S)+(1-E1*E1)(β(N)SI1β(N*)S*I1*-E1E1*)]+pr2I1*[(1-E1*E1)(I1I1*-E1E1*)+C(μ+r2+d1)pr2(1-I1*I1)(E1E1*-I1I1*)]+β(N*)S*I2*[(1-β(N)SI2β(N*)S*I2*)(1-S*S)+(β(N)SI2β(N*)S*I2*-E2E2*)(1-E2*E2)]+(1-p)r2I1*[(I1I1*-E2E2*)(1-E2*E2)+D(μ+d2)I2*(1-p)r2I1*(1-I2*I2)(E2E2*-I2I2*)].
Now, let (x,w1,z1,w2,z2)=(S/S*,E1/E1*,I1/I1*,E2/E2*,I2/I2*) and g(N)=β(N)/β(N*), then we get
(3.58)U̇=-μ(S-S*)2S+β(N*)S*I1*[(1-g(N)xz1)(1-1x)+(1-1w1)(g(N)xz1-w1)]+pr2I1*[(1-1w1)(z1-w1)+C(μ+r2+d1)pr2(1-1z1)(w1-z1)]+β(N*)S*I2*[(1-g(N)xz2)(1-1x)+(g(N)xz2-w2)(1-1w2)]+(1-p)r2I1*[(z1-w2)(1-1w2)+D(μ+d2)I2*(1-p)r2I1*(1-1z2)(w2-z2)],=-μ(S-S*)2S+f(x,w1,z1,w2,z2),
where
(3.59)f(x,w1,z1,w2,z2)=β(N*)S*I1*f1(x,w1,z1)+pr2I1*f2(w1,z1)+β(N*)S*I2*f3(x,w2,z2)+(1-p)r2I1*f4(z1,w2,z2),f1(x,w1,z1)=(1-g(N)xz1)(1-1x)+(1-1w1)(g(N)xz1-w1),f2(w1,z1)=(1-1w1)(z1-w1)+C(μ+r2+d1)pr2(1-1z1)(w1-z1),f3(x,w2,z2)=(1-g(N)xz2)(1-1x)+(g(N)xz2-w2)(1-1w2),f4(z1,w2,z2)=(z1-w2)(1-1w2)+D(μ+d2)I2*(1-p)r2I1*(1-1z2)(w2-z2).
We can choose
(3.60)C=pr2(μ+r2+d1),D=(1-p)r2I1*(μ+d2)I2*.
For x=w1=w2, we have
(3.61)f1(x,w1,z1)=(1-g(N)xz1)(1-1x)+(1-1w1)(g(N)xz1-w1)=g(N)xz1(1x-1+1-1w1)+1-1x-w1(1-1w1)=1-1x-w1+1=2-1x-x,
which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if x=1(3.62)f2(w1,z1)=(1-1w1)(z1-w1)+(1-1z1)(w1-z1)=(z1-w1)(1-1w1-1+1z1)=(z1-w1)(1z1-1w1)=2-w1z1-z1w1,
which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if w1=z1(3.63)f3(x,w2,z2)=(1-g(N)xz2)(1-1x)+(g(N)xz2-w2)(1-1w2)=g(N)xz2(1-1w2+1x-1)+1-1x-w2+1=1-1x-w2+1=2-1x-x,
which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if x=1.
Forw2≥1, z2≥z1 or w2≤1, z2≤z1, we have(3.64)f4(z1,w2,z2)=(z1-w2)(1-1w2)+(1-1z2)(w2-z2)≤(z2-w2)(1-1w2-1+1z2)=(z2-w2)(1z2-1w2)=2-w2z2-z2w2,
which is less than or equal to zero by the arithmetic-geometric mean inequality, with equality if and only if w2=z2.
Thus,U̇(S,E1,I1,E2,I2) is less or equal to zero with equality only if x=1,w1=z1,w2=z2. Since x=w1=w2, the largest invariant set of system (2.2) on the set {(S,E1,I1,E2,I2)∈ℝ+5, U̇(S,E1,I1,E2,I2)=0} is the endemic equilibrium point P2. By LaSalle's invariance principle [31, 32], we can conclude that the interior equilibrium P2 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.


## 3.1. Global Stability of the Disease-Free Equilibrium

The disease-free-equilibrium is given as(3.1)P0=(S0,0,0,0,0),S0=Λμ.
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 P0, we use the results in van den Driessche and Watmough [27]. Reorder the components of P0 as E1=0, I1=0, E2=0, I2=0, S=Λ/μ. Set(3.2)F=[F1F2F3F4F5]=[β1(N)SI10β2(N)SI200],V=[k1(1-r1)E1+μE1-pr2I1r2I1+(μ+d1)I1-k1(1-r1)E1(μ+k2)E2-(1-p)r2I1(μ+d2)I2-k2E2μS+β1(N)SI1+β2(N)SI2-Λ].
The infected compartments are E1, I1, E2, and I2. Thus(3.3)F=[0β1(S0)S0000000000β1(S0)S00000],V=[μ+k1(1-r1)-pr200-k1(1-r1)μ+d1+r2000-(1-p)r2μ+k2000-k2μ+d2].
The dominant eigenvalues of FV-1 are given by (3.4)R1=k1(1-r1)β1(S0)S0(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1),R2=k2β2(S0)S0(μ+k2)(μ+d2),
where R1 and R2 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 FV-1 is(3.5)ρ(FV-1)=max⁡{R1,R2}.
If R1<1, and R2<1, then ρ(FV-1)<1. By Theorem 3.2 in van den Driessche and Watmough [27], we know that the disease-free steady state P0 is locally asymptotically stable, P0 is unstable if R1>1 or R2>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,d1=0, d2=0. Consequently, the total population size N(t) satisfies the equation dN/dt=Λ-μN and N(t)→Λ/μ as t→∞. 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 N=S0=Λ/μ: (3.6)E1′=β1(S0)I1(S0-E1-I1-E2-I2)+pr2I1-(μ+k1(1-r1))E1,I1′=k1(1-r1)E1-(r2+μ)I1,E2′=β2(S0)I2(S0-E1-I1-E2-I2)+(1-p)r2I1-(μ+k2)E2,I2′=k2E2-μI2.
Notice that the S equation is eliminated from (3.6) and the variable S is replaced by S0-E1-I1-E2-I2. By introducing the fractions(3.7)x1=E1N,x2=I1N,y1=E2N,y2=I2N,
we obtain the equivalent limiting system from (3.6) (3.8)x1′=S0β1(S0)(1-x1-x2-y1-y2)x2+pr2x2-(μ+k1(1-r1))x1,x2′=k1(1-r1)x1-(μ+r2)x2,y1′=S0β2(S0)(1-x1-x2-y1-y2)y2+(1-p)r2x2-(μ+k2)y1,y2′=k2y1-μy2.
Obviously(3.9)0≤x1+x2+y1+y2≤1,
for all time t≥0.For a bounded real-valued functionf on [0,∞), we define(3.10)f∞=lim⁡inf⁡t→∞f(t),f∞=limsup⁡t→∞f(t).Lemma 3.1 (see [30]).
Letf:[0,∞)→∞ be bounded and continuously differentiable. Then there are sequences tn,sn→∞ with the following properties:
(3.11)f(tn)⟶f∞,f′(tn)⟶0,f(sn)⟶f∞,f′(sn)⟶0,
for n→∞.Theorem 3.2.
In the absence of disease-induced death rate, the disease-free equilibriumP0 of system (2.2) is globally asymptotically stable if R1<1, R2<1.Proof.
By Lemma3.1 and the x2′, y2′ equations in (3.8) we have
(3.12)x2∞≤k1(1-r1)μ+r2x1∞,y2∞≤k2μy1∞.
Using thex1′ equation in (3.8) and choosing tn→∞ such that(3.13)x1(tn)⟶x1∞,x1′(tn)⟶0,t⟶∞,
we get
(3.14)0≤S0β1(S0)(1-x1-x2-y1-y2)∞x2∞-(μ+k1(1-r1))x1∞≤S0β1(S0)x2∞-(μ+k1(1-r1))x1∞≤S0β1(S0)k1(1-r1)μ+r2x1∞-(μ+k1(1-r1))x1∞≤(μ+k1(1-r1))x1∞[k1(1-r1)S0β1(S0)(μ+r2)(μ+k1(1-r1))-1]≤(μ+k1(1-r1))x1∞[k1(1-r1)S0β1(S0)(μ+r2)(μ+k1(1-r1))-k1pr2(1-r1)-1]=(μ+k1(1-r1))x1∞(R1-1).
It is shown that x1∞≤0 since R1<1, as x1∞≥0, we have that x1∞=0. The inequalities in (3.12) also imply that x2∞=0.
Using they1′ equation in (3.8) and choosing sn→∞ such that(3.15)y1(sn)⟶y1∞,y1′(sn)⟶0,t→∞,
we get
(3.16)0≤S0β2(S0)(1-x1-x2-y1-y2)∞y2∞-(μ+k2)y1∞≤S0β2(S0)y2∞-(μ+k2)y1∞≤S0β2(S0)k2μy1∞-(μ+k2)y1∞≤(μ+k2)y1∞[k2S0β2(S0)μ(μ+k2)-1]=(μ+k2)y1∞(R2-1).
It is shown that y1∞≤0 since R2<1, as y1∞≥0, we have that y1∞=0. The inequalities in (3.12) also imply that y2∞=0.
Hence(3.17)lim⁡t→∞x1(t)=lim⁡t→∞x2(t)=lim⁡t→∞y1(t)=lim⁡t→∞y2(t)=0.
Therefore, the disease-free equilibrium P0 is globally asymptotically stable.Next, we consider that there is disease-induced death rate, which means thatd1>0, d2>0. Then, we can get the following theorem.Theorem 3.3.
Ford1>0, d2>0, the disease-free equilibrium P0 of system (2.2) is globally asymptotically stable if R1<1-k2r2(1-p)/((μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)), R2<1.Proof.
Let us consider the following Lyapunov function:(3.18)V(E1,I1,E2,I2)=k1(1-r1)E1+(μ+k1(1-r1))I1+k2E2+(μ+k2)I2.
Its time derivative along the solutions of system (2.2) satisfies
(3.19)V̇=k1(1-r1)Ė1+(μ+k1(1-r1))İ1+k2Ė2+(μ+k2)İ2=k1(1-r1)[β1(N)SI1+pr2I1-(μ+k1(1-r1))E1]+(μ+k1(1-r1))[k1(1-r1)E1-(μ+d1+r2)I1]+k2[β2(N)SI2+(1-p)r2I1-(μ+k2)E2]+(μ+k2)[k2E2-(μ+d2)I2],=[k1(1-r1)β1(N)S+k1(1-r1)pr2-(μ+d1+r2)(μ+k1(1-r1))]I1+[k2β2(N)S-(μ+d2)(μ+k2)]I2+k2r2(1-p)I1.
Now, using (2.3), one has β1(N)S≤β1(S0)S0, β2(N)S≤β2(S0)S0, then
(3.20)V̇≤[(μ+d1+r2)(μ+k1(1-r1))-k1(1-r1)pr2](R1-1)I1+k2r2(1-p)I1+(μ+d2)(μ+k2)(R2-1)I2,=[[(μ+d1+r2)(μ+k1(1-r1))-k1(1-r1)pr2](R1-1)+k2r2(1-p)]I1+(μ+d2)(μ+k2)(R2-1)I2.
Since
(3.21)(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)=μ(μ+d1+r2)+k1(1-r1)(μ+d1)+k1r2(1-r1)-k1pr2(1-r1)=μ(μ+d1+r2)+k1(1-r1)(μ+d1)+k1r2(1-r1)(1-p)≥0,
and R1<1-k2r2(1-p)/((μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1))≤1, we have
(3.22)[(μ+d1+r2)(μ+k1(1-r1))-k1(1-r1)pr2](R1-1)+k2r2(1-p)≤0.
Thus, V̇≤0 if R1<1-k2r2(1-p)/((μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)), R2<1. It is obvious that V̇=0 if and only if I1=0, I2=0. Then, the largest invariant set of system (2.2) on the set {(S,E1,I1,E2,I2)∈ℝ+5, V̇(E1,I1,E2,I2)=0} is the disease-free equilibrium P0. Therefore, it follows from LaSalle's invariance principle [31, 32] that P0 is globally stable if R1<1-k2r2(1-p)/((μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)), R2<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 R1 and R2.LetP1=(Ŝ,0,0,Ê2,Î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 (3.23)Λ-β2(N̂)ŜÎ2-μŜ=0,β2(N̂)ŜÎ2-(μ+k2)Ê2=0,k2Ê2-(μ+d2)Î2=0,
adding the above three equations, we have(3.24)Λ-μN̂-d2Î2=0.Using (3.24), the first and second equation of (3.23), we can easily express Ŝ, Ê2 and Î2 in terms of N̂ in the form:
(3.25)Ŝ=Λd2μd2+β2(N̂)(Λ-μN̂),Ê2=Λβ2(N̂)(Λ-μN̂)(μ+k2)[μd2+β2(N̂)(Λ-μN̂)],Î2=Λ-μN̂d2.SubstitutingŜ, Ê2, Î2 in the third equations of (3.23) yields(3.26)(Λ-μN̂)F(N̂)=0,
where(3.27)F(N̂)=Λk2d2β2(N̂)-(μ+k2)(μ+d2)[μd2+β2(N̂)(Λ-μN̂)].Clearly,Λ-μN̂=0 is a fixed point of (3.23), which corresponds to the disease-free equilibrium P0. Since N̂∈[0,S0], we get(3.28)F(0)=Λk2d2β2(0)-(μ+k2)(μ+d2)[μd2+Λβ2(0)]=Λk2d2β2(0)-[μ2+μ(d2+k2)+d2k2][μd2+Λβ2(0)]=Λk2d2β2(0)-[μ2+μ(d2+k2)][μd2+Λβ2(0)]-Λd2k2β2(0)-μk2d22=-[μ2+μ(k2+d2)][μd2+Λβ2(0)]-μk2d22,F(S0)=Λk2d2β2(S0)-μd2(μ+k2)(μ+d2)=μ[k2d2β2(S0)S0-d2(μ+k2)(μ+d2)]=μd2(μ+k2)(μ+d2)[k2β2(S0)S0(μ+k2)(μ+d2)-1]=μd2(μ+k2)(μ+d2)(R2-1).It appears thatF(0)<0 and F(S0)>0 when R2>1. The existence of fixed point follows from the intermediate value-theorem. Now, F(N̂) is monotone increasing, so that F(N̂)=0 has only one positive root in the interval [0,S0].Thus, we have established the following resultLemma 3.4.
WhenR2>1, the system (2.2) has a unique drug-resistant-TB-strain-only equilibrium P1=(Ŝ,0,0,Ê2,Î2) with Ŝ, Ê2, and Î2 all nonnegative.


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

Theorem 3.5.
IfR1<1<R2, the drug-resistant-TB-strain-only equilibrium P1 exists and is locally asymptotically stable.Proof.
Reorder the equilibriumP1 as (0,0,Ê2,Î2,Ŝ). Similarly as in the proof of locally stability of the disease-free steady state P0, we have
(3.29)F1=[0β1(N̂)Ŝ00],V1=[μ+k1(1-r1)-pr2-k1(1-r1)μ+d1+r2],F1V1-1=A[k1(1-r1)β1(N̂)Ŝ(μ+k1(1-r1))β1(N̂)Ŝ00],
where
(3.30)A=1(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1).
The spectral radius of F1V1-1 is given by
(3.31)ρ(F1V1-1)=k1(1-r1)β1(N̂)Ŝ(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1),≤k1(1-r1)β1(N̂)N̂(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1).
Since (β1(N)N)′≥0, and N̂≤S0, then
(3.32)β1(N̂)N̂≤β1(S0)S0,
Thus
(3.33)ρ(F1V1-1)≤k1(1-r1)β1(S0)S0(μ+d1+r2)(μ+k1(1-r1))-k1pr2(1-r1)=R1.
As R1<1, we have ρ(F1V1-1)<1, which implies that the drug-resistant-TB-strain-only equilibrium P1 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 equilibriumP1 is global asymptotically stable when R1<1<R2.Proof.
IfR1<1<R2, from Theorem 3.5, P1 is locally asymptotically stable. In the following, we only need to prove that the P1 is a global attractor.
There is no disease-induced death rate, which meansd1=d2=0, using a similar argument as in the proof of Theorem 3.2, we can show that if R1<1,(3.34)lim⁡t→∞x1(t)=lim⁡t→∞x2(t)=0.
Then the equivalent limiting system of (3.6) is
(3.35)E2′=β2(S0)I2(S0-E2-I2)-(μ+k2)E2,I2′=k2E2-μI2.
Let(P,Q) be the vector field defined by system (3.35). Then for the Dulac function D(E2,I2)=1/E2I2, there holds(3.36)∂E2(DP)+∂I2(DQ)=-β2(S0)(S0-I2)E22-k2I22<0.
From Dulac's criterion, we can see equivalent limiting system (3.35) does not have a limit cycle. Therefore, the local stability of P1 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(W1)
limN→0(1/β1(N))=+∞,(W2)
β1(N)/β2(N)=β1/β2,then, existence, uniqueness and global stability of interior equilibrium can be obtained. LetP2=(S*,E1*,I1*,E2*,I2*) be interior equilibrium of system (2.2). Then, the interior equilibrium (steady state with I1>0, I2>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 (3.37)Λ-β1(N*)S*I1*-β2(N*)S*I2*-μS*=0,β1(N*)S*I1*+pr2I1*-k1(1-r1)E1*-μE1*=0,k1(1-r1)E1*-(r2+μ+d1)I1*=0,β2(N*)S*I2*+(1-p)r2I1*-(μ+k2)E2*=0,k2E2*-(μ+d2)I2*=0,
adding the above five equations, we have(3.38)Λ-μN*-d1I1*-d2I2*=0.Using (3.37) and (3.38), we get the following equivalent equations set (3.39)Λ-β1(N*)S*I1*-β2(N*)S*I2*-μS*=0,β1(N*)S*I1*+pr2I1*-k1(1-r1)E1*-μE1*=0,k1(1-r1)E1*-(r2+μ+d1)I1*=0,k2E2*-(μ+d2)I2*=0,Λ-μN*-d1I1*-d2I2*=0.
From (3.39), we have (3.40)S*=β(S0)S0R1β1(N*),I1*=d2(Λ-μS*)-S*β2(N*)(Λ-μN*)d2S*β1(N*)-d1S*β2(N*),E1*=μ+d1+r2k1(1-r1)I1*,I2*=(Λ-μN*)d2-d1d2I1*,E2*=(μ+d2)(Λ-μN*)k2d2-d1(μ+d2)k2d2I1*.Substituting (3.40) into β2(N*)S*I2*+(1-p)r2I1*-(μ+k2)E2*=0, N*satisfies the following equation:(3.41)G(N*)=0,
where (3.42)G(N*)=β1(S0)S0(Λ-μN*)R1[k2d2S0β1(S0)β2(N*)R1β1(N*)-d2μ1-(1-p)k2d2r2β2(N*)β1(N*)]+d2Λ(1-β1(S0)β1(N*)R1)[d1μ1+(1-p)k2d2r2-d1k2S0β1(S0)β2(N*)R1β1(N*)],μ1=(μ+k2)(μ+d2).
Since β1(N)/β2(N)=β1/β2, limN→0(1/β1(N))=+∞, if R1>d1k2β2S0β1(S0)/β1[d1μ1+(1-p)k2d2r2], then(3.43)lim⁡N*→0G(N*)=-∞.
It means that G(N*)=0 has only one positive root in the interval [0,S0] if G(N*) is monotone increasing and G(S0)≥0, which is equivalent to (3.44)R1>H1+H2H3+H4,R2<[1+(1-p)d2k2r2d1μ1]R1,
where (3.45)H1=k2d2β2S0β1(S0)β12(N*),H2=d1d2k2β2S0β1(S0)β1′(N*),H3=β12(N*)[d2β1μ1+d2k2β2r2(1-p)],H4=d2β1′(N*)[d1β1μ1+d2k2β1r2(1-p)].
Let (3.46)H=max⁡{1,H1+H2H3+H4,d1k2β2S0β1(S0)β1[d1μ1+k2d2r2(1-p)]},f(R1)=[1+(1-p)d2k2r2d1μ1]R1.
When(3.47)R1>H,R2<f(R1),
it shows that G(S0)≥0, and G(N*) is monotone increasing, by the intermediate value-theorem, G(N*)=0 has only one positive root in the interval [0,S0], which signify that P2 is unique. Obviously, S*>0. In order to make certain that E1*>0, I1*>0, E2*>0, I2*>0, we obtain that E1*>0, I1*>0, E2*>0, I2*>0 if and only if(3.48)0<I1*<Λ-μN*d1.
Due to d2>d1>0, introducing (3.40) into (3.48), we acquire(3.49)B1<R1<B2,
where (3.50)B1=S0β1(S0)[β2(N*)(Λ-μN*)+μd2]Λd2β1(N*),B2=S0β1(S0)[β1(N*)(Λ-μN*)+μd1]Λd1β1(N*),
and B1<B2.Now, we can get the following conclusion.Lemma 3.7.
When conditions (W1), (W2), (3.47) and (3.49) are satisfied, the system (2.2) has a unique interior equilibrium P2=(S*,E1*,I1*,E2*,I2*) with S*, E1*, I1*, E2*, and I2* all nonnegative.Remark 3.8.
Forβ1(N)=β1 and β1(N)=β1/N, condition (W1) is not satisfied. However, at least when β1(N)=β1C(N), condition (W1) 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 equilibriumP2 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(N)=β2(N)=β(N), then, we get the following result.Theorem 3.9.
If conditions (W1), (W2), (3.47) and (3.49) are satisfied, interior equilibrium P2 of system (2.2) is is globally asymptotically stable when
(3.51)SS*=E1E1*=E2E2*≥1,I2I2*≥I1I1*orSS*=E1E1*=E2E2*≤1,I2I2*≤I1I1*.Proof.
Following [33], we consider the Lyapunov function:
(3.52)U(S,E1,I1,E2,I2)=(S-S*ln⁡S)+(E1-E1*ln⁡E1*)+C(I1-I1*ln⁡I1*)+(E2-E2*ln⁡E2*)+D(I2-I2*ln⁡I2*),
where C and D are positive constants to be determined later. Differentiating this function with respect to time along the solutions of system (2.2) yields
(3.53)U̇=(1-S*S)Ṡ+(1-E1*E1)Ė1+C(1-I1*I1)İ1+(1-E2*E2)Ė2+D(1-I2*I2)İ2=(1-S*S)[Λ-β(N)SI1-β(N)SI2-μS]+(1-E1*E1)[β(N)SI1+pr2I1-k1(1-r1)E1-μE1]+C(1-I1*I1)[k1(1-r1)E1-(r2+μ+d1)I1]+(1-E2*E2)[β(N)SI2+(1-p)r2I1-(μ+k2)E2]+D(1-I2*I2)[k2E2-(μ+d2)I2].
Considering (3.37), we can deduce that
(3.54)Λ=β(N*)S*I1*+β(N*)S*I2*+μS*,μ+k1(1-r1)=β(N*)S*I1*E1*+pr2I1*E1*,μ+d1+r2=k1(1-r1)E1*I1*,μ+k2=β(N*)S*I2*E2*+(1-p)r2I1*E2*,μ+d2=k2E2*I2*.
Then, (3.53) becomes
(3.55)U̇=(1-S*S)[β(N*)S*I1*+β(N*)S*I2*+μS*-β(N)SI1-β(N)SI2-μS]+(1-E1*E1)[β(N)SI1+pr2I1-β(N*)S*I1*E1E1*-pr2I1*E1E1*]+C(1-I1*I1)[k1(1-r1)E1-k1(1-r1)E1*I1I1*]+(1-E2*E2)[β(N)SI2+(1-p)r2I1-β(N*)S*I2*E2E2*-(1-p)r2I1*E2E2*]+D(1-I2*I2)[k2E2-k2E2*I2I2*]=-μ(S-S*)2S+β(N*)S*I1*(1-β(N)SI1β(N*)S*I1*)(1-S*S)+β(N*)S*I2*(1-β(N)SI2β(N*)S*I2*)(1-S*S)+(1-E1*E1)[β(N*)S*I1*(β(N)SI1β(N*)S*I1*-E1E1*)+pr2I1*(I1I1*-E1E1*)]+C(1-I1*I1)[k1(1-r1)E1*(E1E1*-I1I1*)]+(1-E2*E2)[β(N*)S*I2*(β(N)SI2β(N*)S*I2*-E2E2*)+(1-p)r2I1*(I1I1*-E2E2*)]+D(1-I2*I2)[k2E2*(E2E2*-I2I2*)].
Now, using (3.37), we have
(3.56)k1(1-r1)E1*=(r2+μ+d1)I1*,k2E2*=(μ+d2)I2*.
Then, (3.55) can be rewritten as follows:
(3.57)U̇=-μ(S-S*)2S+β(N*)S*I1*[(1-β(N)SI1β(N*)S*I1*)(1-S*S)+(1-E1*E1)(β(N)SI1β(N*)S*I1*-E1E1*)]+pr2I1*[(1-E1*E1)(I1I1*-E1E1*)+C(μ+r2+d1)pr2(1-I1*I1)(E1E1*-I1I1*)]+β(N*)S*I2*[(1-β(N)SI2β(N*)S*I2*)(1-S*S)+(β(N)SI2β(N*)S*I2*-E2E2*)(1-E2*E2)]+(1-p)r2I1*[(I1I1*-E2E2*)(1-E2*E2)+D(μ+d2)I2*(1-p)r2I1*(1-I2*I2)(E2E2*-I2I2*)].
Now, let (x,w1,z1,w2,z2)=(S/S*,E1/E1*,I1/I1*,E2/E2*,I2/I2*) and g(N)=β(N)/β(N*), then we get
(3.58)U̇=-μ(S-S*)2S+β(N*)S*I1*[(1-g(N)xz1)(1-1x)+(1-1w1)(g(N)xz1-w1)]+pr2I1*[(1-1w1)(z1-w1)+C(μ+r2+d1)pr2(1-1z1)(w1-z1)]+β(N*)S*I2*[(1-g(N)xz2)(1-1x)+(g(N)xz2-w2)(1-1w2)]+(1-p)r2I1*[(z1-w2)(1-1w2)+D(μ+d2)I2*(1-p)r2I1*(1-1z2)(w2-z2)],=-μ(S-S*)2S+f(x,w1,z1,w2,z2),
where
(3.59)f(x,w1,z1,w2,z2)=β(N*)S*I1*f1(x,w1,z1)+pr2I1*f2(w1,z1)+β(N*)S*I2*f3(x,w2,z2)+(1-p)r2I1*f4(z1,w2,z2),f1(x,w1,z1)=(1-g(N)xz1)(1-1x)+(1-1w1)(g(N)xz1-w1),f2(w1,z1)=(1-1w1)(z1-w1)+C(μ+r2+d1)pr2(1-1z1)(w1-z1),f3(x,w2,z2)=(1-g(N)xz2)(1-1x)+(g(N)xz2-w2)(1-1w2),f4(z1,w2,z2)=(z1-w2)(1-1w2)+D(μ+d2)I2*(1-p)r2I1*(1-1z2)(w2-z2).
We can choose
(3.60)C=pr2(μ+r2+d1),D=(1-p)r2I1*(μ+d2)I2*.
For x=w1=w2, we have
(3.61)f1(x,w1,z1)=(1-g(N)xz1)(1-1x)+(1-1w1)(g(N)xz1-w1)=g(N)xz1(1x-1+1-1w1)+1-1x-w1(1-1w1)=1-1x-w1+1=2-1x-x,
which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if x=1(3.62)f2(w1,z1)=(1-1w1)(z1-w1)+(1-1z1)(w1-z1)=(z1-w1)(1-1w1-1+1z1)=(z1-w1)(1z1-1w1)=2-w1z1-z1w1,
which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if w1=z1(3.63)f3(x,w2,z2)=(1-g(N)xz2)(1-1x)+(g(N)xz2-w2)(1-1w2)=g(N)xz2(1-1w2+1x-1)+1-1x-w2+1=1-1x-w2+1=2-1x-x,
which is less than or equal to zero by the arithmetic-geometric mean inequality, with the following equality if and only if x=1.
Forw2≥1, z2≥z1 or w2≤1, z2≤z1, we have(3.64)f4(z1,w2,z2)=(z1-w2)(1-1w2)+(1-1z2)(w2-z2)≤(z2-w2)(1-1w2-1+1z2)=(z2-w2)(1z2-1w2)=2-w2z2-z2w2,
which is less than or equal to zero by the arithmetic-geometric mean inequality, with equality if and only if w2=z2.
Thus,U̇(S,E1,I1,E2,I2) is less or equal to zero with equality only if x=1,w1=z1,w2=z2. Since x=w1=w2, the largest invariant set of system (2.2) on the set {(S,E1,I1,E2,I2)∈ℝ+5, U̇(S,E1,I1,E2,I2)=0} is the endemic equilibrium point P2. By LaSalle's invariance principle [31, 32], we can conclude that the interior equilibrium P2 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 Table1. Numerical results are displayed in the following figures.Table 1
Description and estimation of parameters.
ParameterDescriptionEstimated valueΛRecruitment rate of susceptible individuals into the community2 year−1β1(N)Transmission coefficient of sensitive strainVariableβ2(N)Transmission coefficient of resistant strainVariableμNaturally death rate0.01 year−1k1Rate of progression from infection to infectious for sensitive strain0.005 year−1k2Rate of progression from infection to infectious for resistant strain0.005 year−1r1Rate of effective chemoprophylaxis to sensitive strain0r2Rate of effective therapy to sensitive strain0.4 year−1d1Disease-induced mortality rate of sensitive strainVariabled2Disease-induced mortality rate of resistant strainVariablepProportion of sensitive TB individuals received effective therapy that have progression to theE1 class0.9 year−1Taking no account of disease-induced death rate, we assume thatd1=d2=0, β1=0.001, β2=0.0001 (so that R1<1, R2<1). For β1(N)=β1, β2(N)=β2, β1(N)=β1/N, β2(N)=β2/N and β1(N)=β1C(N), β2(N)=β2C(N), C(N)=bN/(1+bN+1+2bN), b=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(N)=β1, β2(N)=β2, β1=0.001, β2=0.0001, d1=d2=0, so that R1<1, R2<1. It shows that the disease-free equilibrium P0 is stable.Figure 3
Trajectories of system (2.2) when β1(N)=β1/N, β2(N)=β2/N, β1=0.001, β2=0.0001, d1=d2=0, so that R1<1, R2<1. It shows that the disease-free equilibrium P0 is stable.Figure 4
Trajectories of system (2.2) when β1(N)=β1C(N), β2(N)=β2C(N), C(N)=N/(1+N+1+2N), β1=0.001, β2=0.0001, d1=d2=0, so that R1<1, R2<1. It shows that the disease-free equilibrium P0 is stable.When there is disease-induced death rate, which means thatd1>0, d2>0, we take d1=0.02, d2=0.07. Firstly, we consider system (2.2) with β1(N)=β1, β2(N)=β2 and choose β1=β2=0.001, so that R1<1, R2<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.001, β2=1, so that R1<1, R2>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.25, so that R1>H, R2<f(R1), 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(N)=β1, β2(N)=β2, β1=β2=0.001, d1=0.02, d2=0.07, so that R1<1, R2<1. It shows that the disease-free equilibrium P0 is stable.Figure 6
Trajectories of system (2.2) when β1(N)=β1, β2(N)=β2, β1=0.001, β2=1, d1=0.02, d2=0.07, so that R1<1, R2>1. It shows that the boundary equilibrium P1 is stable.Figure 7
Trajectories of system (2.2) when β1(N)=β1, β2(N)=β2, β1=1, β2=0.25, d1=0.02, d2=0.07, so that R1>H, R2<f(R1). It shows that the interior equilibrium P2 is stable.Next, we need to think the case whenβ1(N)=β1/N, β2(N)=β2/N. Figures 8, 9 and 10 show the trajectories plot and its plane figure of system (2.2) when β1=β2=0.001 (so that R1<1, R2<1), β1=0.001, β2=1 (so that R1<1, R2>1), and β1=1, β2=0.25 (so that R1>H, R2<f(R1)), respectively. As the previous cases, we can find that the trajectories of system (2.2) converge to the disease-free equilibrium when R1<1, R2<1 (see Figure 8), the drug-resistant-TB-strain-only equilibrium is global asymptotically stable when R1<1, R2>1 (see Figure 9), while the trajectories of system (2.2) converge to the unique interior equilibrium when R1>1, R2<f(R1) (see Figure 10).Figure 8
Trajectories of system (2.2) when β1(N)=β1/N, β2(N)=β2/N, β1=β2=0.001, d1=0.02, d2=0.07, so that R1<1, R2<1. It shows that the disease-free equilibrium P0 is stable.Figure 9
Trajectories of system (2.2) when β1(N)=β1/N, β2(N)=β2/N, β1=0.001, β2=1, d1=0.02, d2=0.07, so that R1<1, R2>1. It shows that the boundary equilibrium P1 is stable.Figure 10
Trajectories of system (2.2) when β1(N)=β1/N, β2(N)=β2/N, β1=1, β2=0.25, d1=0.02, d2=0.07, so that R1>H, R2<f(R1). It shows that the interior equilibrium P2 is stable.Eventually, we assume thatβ1(N)=β1C(N), β2(N)=β2C(N), C(N)=bN/(1+bN+1+2bN) in system (2.2), of which b=1. Numerical results are depicted in Figures 11, 12, and 13. Figure 11 presents the trajectories and its plane figure when β1=β2=0.001 (so that R1<1, R2<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.001, β2=1 (so that R1<1, R2>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.25 (so that R1>H, R2<f(R1)). 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(N)=β1C(N), β2(N)=β2C(N), C(N)=N/(1+N+1+2N), β1=β2=0.001, d1=0.02, d2=0.07, so that R1<1, R2<1. It shows that the disease-free equilibrium P0 is stable.Figure 12
Trajectories of system (2.2) when β1(N)=β1C(N), β2(N)=β2C(N), C(N)=N/(1+N+1+2N), β1=0.001, β2=1, d1=0.02, d2=0.07, so that R1<1, R2>1. It shows that the boundary equilibrium P1 is stable.Figure 13
Trajectories of system (2.2) when β1(N)=β1C(N), β2(N)=β2C(N), C(N)=N/(1+N+1+2N), β1=1, β2=0.25, d1=0.02, d2=0.07, so that R1>H, R2<f(R1). It shows that the interior equilibrium P2 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βi(N), i=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 R1, R2 provide the threshold conditions which determine the competitive outcomes of the two strains. When R1<1, R2<1, the disease-free equilibrium P0 is globally stable and the disease will die out. When R1<1, R2>1, a unique boundary equilibrium P1 exists and is globally stable. While when R1>H, R2<f(R1), there exists a unique interior equilibrium P2 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.



---

*Source: 293747-2011-02-21.xml*