A two-strain TB model with multiple latent stages.

A two-strain tuberculosis (TB) transmission model incorporating antibiotic-generated TB resistant strains and long and variable waiting periods within the latently infected class is introduced. The mathematical analysis is carried out when the waiting periods are modeled via parametrically friendly gamma distributions, a reasonable alternative to the use of exponential distributed waiting periods or to integral equations involving ``arbitrary'' distributions. The model supports a globally-asymptotically stable disease-free equilibrium when the reproduction number is less than one and an endemic equilibriums, shown to be locally asymptotically stable, or l.a.s., whenever the basic reproduction number is greater than one. Conditions for the existence and maintenance of TB resistant strains are discussed. The possibility of exogenous re-infection is added and shown to be capable of supporting multiple equilibria; a situation that increases the challenges faced by public health experts. We show that exogenous re-infection may help established resilient communities of actively-TB infected individuals that cannot be eliminated using approaches based exclusively on the ability to bring the control reproductive number just below 1.

1. TB reservoirs. Tuberculosis (TB) is a bacterial disease caused by the etiological agent Mycobacterium tuberculosis. The World Health Organization (WHO) reports roughly 9.4 million new cases (incidence) per year, an active-TB prevalence of 14 million, and 1.6 to 1.9 million deaths per year, a number that includes 400,000 deaths coming from HIV positive individuals each year. Most active-TB cases are concentrated in South East Asia, African and Western Pacific regions (page 1; [43]). Approximately 250,000 multi-drug TB resistant (MDR TB) cases (including increasing cases of extremely drug resistant or XDR-MDR TB), most involving undiagnosed (nearly ninety percent) individuals, were reported by 2009 and the number seems to be increasing. Funding to control TB has passed the 11 billion mark (US dollars) (page 1; [43]). Why has it been so difficult to control TB?
"Latent tuberculosis infection (LTBI) ... a state of persistent immune response to stimulation by Mycobacterium tuberculosis antigens without evidence of clinically manifested active TB" (page 10; [19]) is extremely prevalent, a huge TB reservoir. "... One third of the world's population is estimated to be infected with M. tuberculosis ... are at risk for developing active TB disease and becoming infectious. The lifetime risk of [TB reactivation] for a person with documented LTBI is estimated to be 5 − 10%, with the majority developing TB disease within the first five years after initial infection ... the risk of developing TB disease following infection depends on several factors, the most important one being the immunological status of the host" (page 10; [19]) Tuberculosis (TB) has often simplistically but effectively been classified as a fast/slow disease since there is a small window in time, following infection, that may give rise to active TB cases. However, in general, TB-activation (and re-activation) is slow and rare [28]. Here, it is assumed that TB may be activated at "all" speeds, with most active-TB cases the result of extremely slow rates of activation at this time in history. TB progression has, over the past centuries, moved from fast (and common) to primarily slow (and rare) progressing disease. Hence, the use of flexible distributions to model TB latency periods is natural and necessary. Today, active-TB (infectious) individuals (8−14 million per year globally) remain so for relatively short periods of time due to the generalized availability of antibiotics. Yet, active TB disease remains a serious global health threat, the result of a lack of universal access to antibiotics and the high levels of non-compliance with antibiotic treatment when available. Further, the emergence and/or re-emergence of diseases that have a strong deleterious impact on the immunological status of individuals (HIV and Ebola) have increased the importance that of global threat posed by TB. In short, TB progression, the transfer or movement of infected latent individuals to the active stage, is today relatively uncommon (given the size of the reservoir), taking place on time scales not too different than the average life span of human hosts albeit immune-compromising diseases are known to facilitate the reactivation of latent TB, with malnutrition playing a central role for decades. [2,3,9,10,12,32,36,38,42].
In short, there is a huge TB-latent human reservoir, for Mycobacterium tuberculosis [6]. This reservoir has increased its average probability of re-activation due to the emergence and growth HIV and TB drug-resistant strains; data supporting the premises of this article, which can be summarized in two questions. What is the role of long and variable latency periods and treatment on the co-evolutionary dynamics of MDR-TB and TB? and what is the role if exogenous re-infection? 2. Drug resistance. One of the biggest health challenges facing the world is tied in to the dramatic increases in the levels of drug resistance, particular in hospital settings [18,31]. The case of Tuberculosis, a disease heavily dependent on multidrug treatment that must be faithfully carried out over 6-9 months [10], is a source of concern due the fact that a number of TB-active individuals do not complete treatment giving rise to the emergence of drug resistance TB strains [11,23,27]. In this manuscript, we study the dynamics of TB-drug resistance in a setting that expands earlier works [7,11,12,28]. In this manuscript, the latency period distributions for the sensitive and resistant TB strains are modeled using gamma distribution to better account for the length and variability of these periods on TB dynamics. This paper is organized as follows: Section 3 introduces the general model, identifies some key epidemiological parameters and averages; Section 4 carries the basic analysis of special sub-models where the latency periods for the drug sensitive and resistant strains have different distributions (not the same gamma distribution), conditions for co-existence and competitive exclusion are identified using appropriately defined and identified threshold conditions; Section 5 proceeds with the full analysis of the general model under the assumption that drug sensitive and drug resistant strains are given by distinct gamma distributions; numerical simulations are used to highlight the effects of using distinct gamma distributions; Section 6 looks at the impact of exogenous re-infection in two-strain models including the implications for public health policy while Section 7 collects conclusions and discusses the implications of the models' results.
3. Model formulation. The total population at time t, denoted by N (t), is subdivided into five disjoint classes of susceptible (S(t)), latent (L(t); with n latent stages), infectious (I(t)), drug-resistant latent (M(t); with m latent stages), and drug-resistant infectious (J(t)) individuals, so that It is assumed that susceptible humans are recruited into the population at a constant rate Λ while susceptible individuals acquire TB infections at the per-capita rate and drug-resistant infections at the per-capita where β 1 and β 2 denote susceptibility to TB and drug-resistant TB infections, respectively. It is assumed that latent individuals do not transmit infection and that the per-capita natural death rate in each class is given by the constant µ. Thus, the rate of change of the susceptible population is given by The population of latent individuals in stage 1 (L 1 ) is generated via the infection of susceptible individuals (at the rate λ 1 ) and via the fraction q of the rate of treated infectious individuals, with r denoting the per-capita treatment rate for infectious individuals. This latent population progresses to the next latent stage (L 2 ; at a rate a 1 σ) while experiencing the per-capita natural death µ, and so, In general, population of latent individuals in stage i (with 2 ≤ i ≤ n) come from the progression of individuals from the previous stage. That is, those in stage L i−1 move into the stage i at the per-capita rate a i−1 σ. This population looses individuals via progression to the next latent stage at the per-capita rate a i σ and by natural death at the per-capita rate µ, so that The population of infectious individuals comes from the latent individuals in the final stage (n). It is assumed that they develop symptoms at this stage at the percapita rate a n σ while loosing individuals via antibiotic treatment at the per-capita rate r, or via natural death at the per-capita death rate µ, or via disease-induced death at the per-capita rate d 1 . Thus, we have, dI dt = a n σL n − (r + µ + d 1 )I.
The population of drug-resistant latent individuals in stage 1 (M 1 ) come from the population of susceptible individuals acquiring drug-resistant infections at the per susceptible and per infective (with resistant strain) rate λ 2 . It looses individual via progression to the next drug-resistant latent stage M 2 ; at the per-capita rate b 1 η while dying from natural causes at the per-capita death µ, so that The population of drug-resistant latent individuals in stage k (with 2 ≤ k ≤ m) comes from the progression of individuals in stage k − 1 at the per-capita rate b k−1 η. It losses individuals via progression to the next drug-resistant latent stage at per-capita rate b k η or by natural death at the per-capita rate µ. Thus, Finally, the population of drug-resistant infectious individuals is generated when drug-resistant individuals in the final stage (m) develop symptoms at the per-capita rate b m η with the fraction (1 − q) of treated infectious individuals developing drugresistant TB (due to non-compliance with treatment). This class looses individuals by natural death at the per-capita rate µ and via disease-induced death at the per-capita rate d 2 , so that In the above formulation, a i , b k (i = 1, 2, · · · , n; k = 1, 2, · · · , m) are assumed to be constants and so the distributions of latent periods at each stage is exponential. That is, P Li = a i σe −aiσS , i = 1, 2, · · · , n, In Eq. (3), T Li = 1 aiσ and T M k = 1 b k η denote the mean stage latent periods. Using the definitions in Eq.
representing the mean time spent in the latency stage (e.g., 1 σ for the latent class, L), which is "shared" among the various stages with the time period 1 σ distributed equally (if a 1 = a 2 = · · · = a n = n) or unequally (if a 1 = a 2 = · · · = a n = n) among all the L i (i = 1, 2, · · · , n) stages. Letting, and making use of Eqs. (3) and (4), it follows that the time spent in the the classes L and M are gamma-distributed. Specifically, their corresponding respective distributions are with the mean associated latent periods given, respectively, by Putting the definitions, formulations and assumptions together, it follows that a model for the transmission dynamics of TB, an infectious disease that involves a wild and resistant (antibiotic generated) strains, must involve latent, infectious, drug-resistant latent and drug-resistant infectious periods, that when modeled via gamma distributed sojourn periods, give rise to the following non-linear system of nonlinear differential equations (model flow diagram is given in Figure 1 and the associated variables and parameters are tabulated in Table 1).

Symbol Description
Variables Population of latent individuals Population of treated individuals who developed drug resistance Parameters Λ Recruitment rate µ Natural death rate β 1 , β 2 Effective contact rate needed to acquire a TB and a drug-resistant TB infections, respectively r Per-capita treatment rate for infectious individuals q Fraction of treated infectious individuals who complete treatment a i σ Progression rate from latent stage i to stage i + 1(i = 1, · · · , n − 1) a n σ Progression rate of latent individuals in stage n to infectious class b k η Progression rate from drug-resistant latent stage k to stage k + 1(k = 1, · · · , m − 1) b m η Progression rate of drug-resistant latent individuals in stage m to drug-resistant infectious class d 1 , d 2 TB-induced mortality rate for individuals in I and J classes, respectively Table 1. Description of variables and parameters of the model (5).
3.1. Basic model properties. In this section we show that our model is "biologically" well posed.
Proof. Let It follows from the first equation of Model (5) that which can be written as, Thus, so that, Similarly, it can be shown that L i (t) ≥ 0, I(t) ≥ 0, M k (t) ≥ 0, and J(t) ≥ 0 (with i = 1, · · · , n, k = 1, · · · , m )for all time t > 0. Hence, all solutions of Model (5) remain positive for all non-negative initial conditions, as required. Proof. Consider the biologically-feasible region, D and observe that the rate of change of the total population, obtained by adding all the equations of the model (5), is given by It follows that dN dt < 0 whenever N > Λ µ . Further, since dN dt ≤ Λ − µN , it is clear that, using a standard comparison theorem [24] that .
Hence, all solutions of the model with initial conditions in D do remain in D for all t > 0 (i.e., the ω-limits sets of System (5) are contained in D), that is, the set D is positively-invariant and attracting.
Since the region D is positively-invariant, the unique solution of the model (5) exists and depends continuously on the initial data of the model (hence, it is sufficient to study its asymptotic dynamics in the region D [21]). 4. Special model case: n = 2 and m = 1. First, we focus on the analysis of the special case of Model (5) obtained when n = 2 and m = 1, namely, Making use of the next generation operator method [14], we compute R C and address the conditions that guarantee the local stability of ε 0 . The non-negative matrix, F C , of the new infection terms, and the M -matrix, V C , of the transition terms associated with the model (6), are given, respectively, by It follows, from Theorem 2 in [14], that the control reproduction numbers of Model (6), defined by R C = ρ(F C V −1 C ) (where ρ denotes the spectral radius of the next generation matrix F C V −1 C ) are given by: and R C1 and R C2 are control reproduction numbers for the TB wild strain only and for the drug-resistant TB strain, respectively. Thus, the spectral radius of the matrix survives and proceeds to stage I. A fraction h 2 = qr r+µ+d1 of infectious individuals re-enters compartment L 1 . Hence, a fraction h 1 of latent individuals passes through compartment I at least once, a fraction h 2 1 h 2 passes through at least twice, and a fraction h k 1 h k−1 2 passes through at least k times. The expected number of times an latent individual passes through compartment I is Thus, the expected total time that a newly-infected individual spends in compartment I is given by Thus, the generation of secondary TB-strain infections generated by a "typical" infectious individual when an appropriate fraction of population is receiving treatment is given by

4.2.2.
Terms in the expression for R C2 . Similarly, the average duration in the class J is given by T J = 1 µ+d2 , the newly-infected individuals entering compartment M with a fraction b1η (b1η+µ) surviving and proceeding to stage J. Thus, the expected total time that a newly-infected individual spends in compartment J is given by .
Thus, the generation of secondary TB resistant-strain infections generated by a "typical" resistant-strain infected individual given that an appropriate fraction of the population is receiving treatment, is given by  (6) given by (7), is locally-asymptotically The epidemiological implications associated with Theorem 4.1 leads to the conclusion that TB can be controlled in a community provided R C < 1 and the initial sizes of the sub-populations involved in Model (6) are in the basin of attraction of the DFE. The dimensionless threshold quantity or ratio, R C , represents the average number of new TB infections generated by a single infectious individual in a susceptible population in the presence of treatment [14].
The proof of Theorem 4.2, is given in Appendix A.

4.3.
Existence of an endemic equilibrium point (EEP). Computing endemic states are important, particularly in the case of TB, since roughly a third of the world population live in a TB-latent stage. We divide the analyses in various cases as follows: Case (i): The case q = 1, that is, treatment is 100% effective. It is assumed, that at the time (t = 0), a subpopulation of drug-resistant infectious individuals exists. That is, J(0) > 0 (due to immigration or due to the inadequate prior use of the drugs used for treatment). We see that when q = 1, two possible endemic equilibria for the System (6) can be identified, namely, two boundary equilibria ε 1 (when only the first strain is present) and ε 2 (when only the second strain is present). There is no interior equilibrium (where both strains are present).

4.3.1.
The TB-strain only equilibrium. This is obtained by setting M = J = 0. The TB only equilibrium in terms of the equilibrium value of the force of infection λ * * 1 is given as where, The expression for λ 1 , defined in Eq. (1), at the endemic steady-state, denoted by λ * * 1 , is given by where N (t) is now replaced by its limiting value, N * * = Λ µ . Substituting equations in Eq. (12) into Eq. (13) and simplifying leads to the following equation in terms of λ * * 1 , Thus, The components of the unique endemic equilibrium (ε 1 ) can be obtained by substituting the unique value of λ * * 1 , given in Eq. (14), into the expressions in (12). Thus, the following result has been established. Lemma 4.3. Model (6) has a unique TB-strain-only equilibrium, given by The proof of Theorem 4.4 is given in Appendix B.

4.3.2.
The drug resistant TB-strain only equilibrium. This is obtained by setting L 1 = L 2 = I = 0. The drug resistant TB only equilibrium in terms of the equilibrium value of the force of infection λ * * 2 is given as where, The expression for λ 2 , defined in Eq. (2), at the endemic steady-state, denoted by λ * * 2 , is given by Substituting the equations in (16) into Eq. (17), and simplifying, we have (in Thus, λ * * 2 > 0, whenever R C2 > 1. Therefore, the following result has been established. Lemma 4.5. Model (6) has a unique drug resistant TB-strain-only equilibrium, given by (ε 2 ) whenever R C2 > 1 > R C1 .
The proof of Theorem 4.6 is given in Appendix C. Figure 2 captures the bifurcation possibilities for the case q = 1. There are three regions I, II, and III in the parameter space (R C1 , R C2 ). In the region I, ε 0 is a global attractor and other equilibria are unstable when they exist. In regions II and III, ε 3 does not exist and ε 1 and ε 2 are l.a.s., respectively. Case (ii): 0 < q < 1 treatment is not 100% efficient due to non-compliance There are two possible endemic equilibria for System (6) when 0 < q < 1, namely, ε 2 (when only the second strain is present), and the interior equilibrium point ε 3 (when both strains exist). In the last case, there is no boundary equilibrium ε 1 , that is, there is no equilibrium where only the first strain is present.

4.3.3.
The drug resistant TB-strain only equilibrium. Analysis of ε 2 in the case 0 < q < 1 is identical in the case q = 1, but it is not shown here to avoid repetition. Thus, the following results have been established.

4.3.4.
Interior equilibrium. The endemic equilibrium where both TB strains co-exist is denoted by where, The expression for λ 1 , defined in (1), at the endemic steady-state, denoted by λ * * 1 , is given by Substituting the equations in Eq. (20) into Eq. (21), and simplifying, gives the following equation (in terms of λ * * 1 and λ * * 2 ) Since the left-hand side of Eq. (22) is always positive, it is necessary that R C1 > 1. The expression for λ 2 , defined in Eq. (2), at the endemic steadystate, denoted by λ * * 2 , is given by Substituting the equations in Eq. in (20) into Eq. (23), and simplifying, gives the following equation (in terms of λ * * 1 and λ * * 2 ) so that, Further, for (ε 3 ) to exist, it is necessary that the TB and resistant strains exist (i.e., R C1 > 1 and R C2 > 1). Thus, we have established the following result: Lemma 4.9. The model (6) has a unique interior equilibrium, given by ε 3 whenever R C1 > R C2 > 1.
The proof of Theorem 4.10 is given in Appendix D. Figure (3) gives a bifurcation diagram for the case 0 < q < 1. There are three regions I, III and IV in parameter space (R C1 , R C2 ) (ε 1 does not exist.). The regions correspond to zones of stability for ε 0 , ε 2 , and ε 3 , respectively. When q → 0, the region IV becomes wider (i.e. more drug resistant individuals exist there and so "more" coexistence). Whenever q → 1, the region IV narrows (i.e. "less" drug-resistant individuals exist there).

5.
Analysis of the full model. In this section, System (5) (i.e., the full model) is analyzed, that is, we look at its equilibria and its stability properties.

5.1.
Local stability of disease-free equilibrium (DFE). Model (5) has a unique disease-free equilibrium, which is obtained by setting the right-hand sides of the equations in the model (5) equal to zero and solving and so, we get that,  (6) in the case 0 < q < 1.
We proceed to compute the control reproduction number via the next generation operator method; and so, we make use of the associated matrices F C , for the new infection terms, and V C , for the remaining transition terms. They are given by and where, C F is a (m + 1) × (n + 1) zero matrix, B F , B V are (n + 1) × (m + 1) zero matrices. Further, A F is a (n + 1) × (n + 1) matrix and D F is a (m + 1) × (m + 1) matrix given by, · · · −a n−1 σ a n σ + µ 0 0 0 · · · 0 −a n σ r

given by
It follows that where R C1 and R C2 are reproduction numbers for TB strain only and drug-resistant TB strain only, respectively. Hence, that the spectral radius of the matrix

5.2.
Interpretation of R C .

5.2.1.
Terms in the expression for R C1 . The average duration that an individual spends in class I per visit to this class is given by T I = 1 r+µ+d1 . Of the newlyinfected individuals entering into compartments L i , i = 1, · · · , n, the fraction survives and proceeds to stage I. The fraction h 2 = qr r+µ+d1 of infectious individuals re-enters compartment L 1 . Hence, the fraction h 1 of latent individuals passes through compartment I at least once, the fraction h 2 1 h 2 passes through at least twice, and the fraction h k 1 h k−1 2 passes through at least k times. Hence, the expected number of times that a latent individual passes through compartment I is Thus, the expected total time that a newly-infected individual spends in compartment I, is given by Hence, the number of secondary new TB-strain infections generated by a TB-active individual is given by

5.2.2.
Terms in the expression for R C2 . Similarly, the average duration in the class J is given by survives and proceeds to stage J. Thus, the expected total time that a newlyinfected individual spends in compartment J is given by Hence, the number of secondary new TB resistant-strain infections generated by a TB-resistant active case individual is given by Therefore, we have arrived at the following results: Theorem 5.1. The DFE, ε 0 , of Model (5), given by (25), is l.a.s. if R C1 < 1 and Theorem 5.2. The DFE, ε 0 , of the model (5) is g.a.s. in D whenever R C1 < 1 and R C2 < 1.
The proof of Theorem 5.2 is given in Appendix E.

Existence of endemic equilibrium point (EEP).
Case (i): q = 1, that treatment is 100% effective It is assumed that at the initial time (t = 0) there are drug-resistant infectious individuals in the population, that is, J(0) > 0. There are two possible endemic equilibria for System (6) when q = 1, two boundary equilibria ε 1 (when only the first strain is present) and ε 2 (when only the second strain is present). There is no interior equilibrium (where both strains are present).

5.3.1.
The TB-strain only equilibrium. This is computed by setting M 1 = · · · = M m = J = 0 and observing that the TB-only equilibrium is a function of the force of infection λ * * 1 . Specifically, it is given by where, L * * n = r + µ + d 1 a n σ I * * , The expression for λ 1 , defined in (1), at the endemic steady-state, here denoted by λ * * 1 , which is given by where N (t) is now replaced by its limiting value, N * * = Λ µ . Substituting the equations in Eq. (28) into (29) and simplifying, leads to the following equation in terms of λ * * 1 : The components of the unique endemic equilibrium (ε 1 ) are obtained by substituting the unique value of λ * * 1 , given in (30), into the expressions in (28). Hence, the following result has been established. Lemma 5.3. Model (5) has a unique TB-strain-only equilibrium given by ε 1 whenever R C1 > 1 > R C2 .
The proof of Theorem 5.4 is given in Appendix F.

5.3.2.
Drug resistant TB-strain equilibrium. An expression of the drug resistant TB-strain equilibrium is obtained by setting L 1 = · · · = L n = I = 0 from where this equilibrium can be expressed in terms of the equilibrium value of the force of infection λ * * 2 . It is given by The expression for λ 2 , defined in Eq. (2), at the endemic steady-state, denoted by λ * * 2 , is given by Substituting the equations in (32) into Eq. (32) and simplifying, leads to the following equation in terms of λ * * 2 : So that, Thus, the following result has been established: Lemma 5.5. Model (5) has a unique drug resistant TB-strain-only equilibrium, which it is given by ε 2 whenever R C1 < 1 < R C2 .
The proof of Theorem 5.6 is given in Appendix G. Case (ii): 0 < q < 1, that is, treatment is not 100% effective due to noncompliance.
There are two possible endemic equilibria for model system (6) when 0 < q < 1: ε 2 (when only the second strain is present), and the interior equilibrium point ε 3 (when both strains exist). In this case there is no second boundary equilibrium ε 1 (where only the first strain is present).

5.3.3.
The drug resistant TB-strain only equilibrium. Analysis of ε 2 in the case 0 < q < 1 is identical to the case q = 1 and so it is not included. Thus, the following result is true.

5.3.4.
Interior equilibrium. The endemic equilibrium where both TB strains exist is denoted by A TWO-STRAIN TB MODEL 759 L * * n = r + µ + d 1 a n σ I * * , The expression for λ 1 , defined in (5), at the endemic steady-state, is denoted by λ * * 1 . It is given by Substituting the equations in (34) into Eq. (35) and simplifying, leads to the following equation, in terms of λ * * 1 and λ * * 2 : Since the left-hand side of Eq. (36) is always positive, it is necessary that R C1 > 1.
The proof of Theorem 5.10 is given in Appendix H.
with N = 25000. The parameters of model (6) take the values of Table 3.
The numerical simulations of Model (6) when R C < 1(R C1 < 1 and R C2 < 1) (Figure 8), suggest that the associated disease-free equilibrium (ε 0 ) is stable when it exists. We took β 1 = 0.2 and β 2 = 0.1, corresponding to R C1 = 0.5752 and R C2 = 0.4452, with the rest of the parameters taking the values in Table 3. Figure 9 shows that, for R C1 > 1 and R C2 < 1, the TB-strain-only equilibrium (ε 1 ) exists. We considered q = 1, β 1 = 0.8 and β 2 = 0.1, which corresponds to R C1 = 2.3007 and R C2 = 0.4452. The rest of the parameters take the values in Table 1. We observe that the state variables converge to (ε 1 ) when t → ∞. Figure 10 shows that, for R C1 < 1 and R C2 > 1, the drug resistant TB-strainonly equilibrium (ε 2 ) exists. We considered q = 0.5, β 1 = 0.2 and β 2 = 0.6, corresponding to R C1 = 0.1484 and R C2 = 2.6709, and the rest of the parameters take the values in Table 1. We observe that the state variables converge to (ε 2 ) when t → ∞. Figure 11 shows that the interior equilibrium (ε 2 ) exists when R C1 > 1 and R C2 > 1. We took q = 0.5, β 1 = 2 and β 2 = 0.6, values that correspond to R C1 = 1.4840 and R C2 = 2.6709. The rest of the parameters take the values in Table 1. We observe that the state variables converge to (ε 2 ) when t → ∞.
The cumulative number of new cases, as a function of time, for various values of the effective contact rate for TB infection (β 1 ), it is shown in Figure 12. It follows from these figures that, as expected, an increase of effective TB contact rates significantly increase cumulative disease incidence. The effect prevalence of the disease is depicted in Figure 14.
The effect of the effective contact rate for drug-resistant TB infection (β 2 ) on the cumulative incidence and prevalence of the disease is depicted in Figures 13 and  15, respectively.

5.5.
On models with parametrically friendly sojourn distributions. Sections 3, 4, and 5 introduced two-strain TB models that included long and variable periods of latency, which are characteristic of TB. These models expand on those first introduced in [10,11,16]. The results, in the context of single strain SEIR models, are those found in [17], here generalized two include two-strains. The results in this paper combined the results for two strains found in [10] and single strain found in [17]. The formulae includes explicit expression for the basic reproduction number for generalized strain-specific gamma distributions and show that the results carry out in the case of two strains in the presence of antibiotic resistance due to non-treatment compliance. The modeling and the analysis has been used in similar settings before. A good reference that includes a detailed analysis of quite general models is that of [22]. 6. Exogenous reinfections. Exogenous reinfection has been identified as a possibly significant mechanism in the transmission dynamics of TB, particularly in environments where tuberculosis infectivity and prevalence are high. Mathematical models that focus on the impact of exogenous reinfection on TB dynamics have been developed and analyzed [15,41,42].
Some experience TB recurrence, that is, the return of symptoms, active TB, after the effect of treatment is over. In fact, endogenous infection and exogenous reinfection are well-documented paths that lead to active tuberculosis among individuals who have experienced prior active-TB infections. Whether or not the major mechanisms behind the recurrence of TB include endogenous reactivation (exacerbation of an old infection) or exogenous reinfection (recurrent infection by a different strain) has been debated for decades [25]. What does activate TB? Endogenous reactivation is the leading suspect in the search for significant mechanisms responsible for tuberculosis progression since the 1960s [39]. TB models have incorporated endogenous reactivation (relapse) via the generation of secondary cases from contacts between treated individuals and those with active infections [1,4,7,15,40,41]. However, the evidence that exogenous reinfection is also an important source of TB recurrence among "recovered" is strong [2,4,13,15,33,41]. There are studies supporting this mode of infection. For example, a recent Shanghai study found that 61.5% of TB recurrent cases from 1999 through 2004could be attributed to exogenous reinfection [33] (see also exogenous reinfection as reported in [5,26,29,30, 34]) 6.1. Model equations. Exogenous re-infection can be achieved with minimal second order progression interactions between TB-exposed and TB-active. The impact of this second order interactions are rather surprising and unique [15,41,42] and so, the goal here is to highlight its dynamical role in a model that includes sensitive and resistant TB strains, the first time that this has been explored in two strain models. And so, two latent classes are considered, drug-sensitive and drug-resistant strains within a highly simplified version of the general model proposed in the prior sections, since the primary goal, is that of highlighting the impact of second order interactions, exogenous re-infection to be precise, have on the dynamical richness of the TB models just studied. Members in the L 1 latent class have a higher risk of developing active TB at the rate α 0 L 1 (fast route to active TB) but in general they will move at the rate a 1 σL 1 to the next latent latent stage. It is assumed that exogenous reinfection does not act in a substantive way on the L 1 class because the time spent in this class is very small, and so this effect is totally neglected in this simplified model. And so, since individuals tend to spend long periods of time in the latent state, it is assumed that the average residence time in L 2 is significant. Consequently, individuals in this class may be reinfected, moving into the TB-active class, via the exogenous reinfection route. We also assume that the drugsensitive strain will not play a role in the process of exogenous reinfection for the drug-resistant strain (another simplifying assumption). Time scale considerations makes it possible to assume that successfully treated individuals will become members of the L 2 class at the rate qrI. It is assumed that α 0 << a 1 σ since only 10% of TB infected individuals move into the active TB cases class. In addition to new the parameter α 0 , we add two more parameters β 3 and β 4 to account for the risk of exogenous reinfections under the just described simplifying assumptions. The model equations with exogenous reinfections are therefore given by the following nonlinear system, A flow diagram of Model (39) is found in Figure 4.  (39) 6.2. Computation of the basic reproduction number. The next-generation matrix approach [14] is used again to compute the basic reproduction number.
To cluster blocks of zeroes in the linearization matrices are used to rearrange the equations in Model (39) as follows.
In block form

Using the invertible matrices
we can express V −1 as The expression r c = r (1−q)a2σ+µ a2σ+µ gives the effect of treatment on drug-sensitive TB infections. Since successfully treated individuals go back to the chronically latentlyinfected class L 2 class rather than the susceptible class), the effect of treatment rate needs to be discounted by the factor of (1−q)a2σ+µ a2σ+µ < 1. We find that We introduce a new aggregate constants C = (µ + d 2 )(r c + µ + d 1 ) to simplify the expressions. The explicitly expression for v 21 , a crucial component of the inverse of V, is given by Since FV −1 = βv 21 βv 22 0 0 the dominant eigenvalue of FV −1 is the same as the dominant eigenvalue of the sub-matrix given by the first two rows and first two columns of βv 21 , that is, Hence, we can define and And so, the basic reproduction number is In order that we gain a better understanding of the basic reproduction number, we rearrange R 1 0 and R 2 0 as follows: R 1 0 is the basic reproduction number for the drug-sensitive strain and R 2 0 for drugresistant strain. Each term in R 1 0 and R 2 0 has a clear biological meaning. We provide a detailed interpretation of R 1 0 for illustrative purposes. Its first part is β1 rc+µ+d1 α0 a1σ+α0+µ giving the secondary cases generated via the fast route. The term 1 rc+µ+d1 stands for the adjusted average infectious period of drug-sensitive strain TB-active individuals. The term α0 a1σ+α0+µ gives the "probability" of progressing to the TB-active state before moving into the second latent stage or dying. Its second part, β1 rc+µ+d1 a1σ a1σ+α0+µ a2σ a2σ+µ gives the secondary cases generated via the slow route. The infected in L 1 class move to L 2 class with "probability" a1σ a1σ+α0+µ and the infected in the L 2 class progresses into TB-active class with "probability" a2σ a2σ+µ . The average infectious period is assumed to be same as in the fast route, namely, 1 rc+µ+d1 .
A TWO-STRAIN TB MODEL 765 6.3. Bifurcation classification at R 0 = 1. This section classifies the bifurcation of Model (39) that originates at R 0 = 1. The following theorem (see [12]) helps to identify the bifurcation type.
Theorem 6.1. Consider a system of ordinary differential equations with a parameter φ.
This theorem can be used in a straightforward way after we relabeled the variables and rewrite System (39) into the following way, The Jacobian matrix of System (47) around the DFE (Λ/µ, 0, 0, 0, 0, 0) is given by We shall deal with the case R 1 0 ≥ R 2 0 or R 2 0 ≥ R 1 0 in exactly the same way. . It can be checked that is a right eigenvector of J(R 0 = 1) corresponding to eigenvalue 0. With little longer algebraic management then finding the right eigenvector, we can find a left eigenvector of J(R 0 = 1) corresponding to eigenvalue 0 is given by We compute a using (45) as follows.
In processing these computations, equation β2 µ+d2 b1η b1η+µ = 1 was used. Choosing β 2 as the bifurcation parameter, we can check that Therefore, we have established the following theorem. Comparing Model (39) with the earlier models discussed in this paper we find that the bifurcation at R 0 = 1 has changed. Theorem 6.2 essentially concludes that whenever exogenous reinfections are incorporated, the basic reproduction number alone is not enough to determine dynamical outcomes. The initial values must now be considered since there are a pair of stable equilibria. 6.3.2. The case of R 1 0 ≥ R 2 0 . We find that if R 1 0 ≥ R 2 0 then a backward bifurcation is possible. In this case R 0 = R 1 0 > R 2 0 . Now we switch the bifurcation parameter from β 2 to β 1 . One can verify that R 0 = 1 if only if be a left eigenvector of J corresponding to the eigenvalue 0 with . and observe that a right eigenvector of J corresponding to the eigenvalue 0 is W = [w 1 , w 2 , w 3 , w 4 , w 5 , w 6 ] with Using W and V , we are able to compute a in (45) Since we only concerned here with the sign of a, we focused on a positive multiplier of a instead, namely Noticing that w4 w5 = a1σ(rc+µ+d1) a1σ2σ+(a2σ+µ)α0 + qr a2σ+µ , it can be checked that β 1 v 2 = r c + µ + d 1 using R 0 = R 1 0 = 1. Thus, aΛ Straightforward computations lead to the expressions Simplifying the summation term gives Since β 1 is the selected the bifurcation parameter, we can check the sign of b via the expression Parallel to Theorem 6.2, we arrive to the following theorem: then Model (39) undergoes a backward bifurcation at R 0 = 1, while if Model (39) undergoes a forward bifurcation at R 0 = 1. Theorem 6.2 and Theorem 6.3 show that if exogenous reinfection forces (β 3 or β 4 ) reach some critical values then a backward bifurcation can happen, creating an endemic equilibrium when R 0 < 1. These cases are numerically illustrated in Figures 5 and 6, respectively, by choosing not necessarily realistic parameter values. In such cases, the disease will not go away when R 0 < 1. Definitely, latently-infected individuals must avoid becoming infected again, killing the possible occurrence of a backward bifurcation. 6.4. Multiple endemic equilibria. The analysis of the case of R 0 > 1 gets us into pretty rough algebra and so, we investigate the model numerically, particularly documenting the appearance of multiple endemic equilibria since the goal is just to show that the introduction of the possibility of exogenous reinfection within drug-sensitive and drug-resistant strains, in highly simplified circumstances, can indeed support complex dynamical outcomes, particularly the possibility of multiple endemic states . In order to just illustrate these possibilities, we choose the following (unrealistic) parameter values, d 1 = .01, d 2 = .01, Λ = 1, µ = .14, β 3 = 0.3, β 4 = 0.3, q = 0.2, β 1 = 2, a 1 = 1, a 2 = 2, σ = 0.04, η = 0.03, b 1 = .01, r = 0.01, and α 0 = 0.01. Hence, we let β 2 vary from 1.5 to 20 as we generate the bifurcation curves in Figure 7. For this set of parameters we have that R 0 = 1.6222 > 1 and yet as it can be seen from Figure 7, there are four equilibria, two of which are stable.  In this Section 6, we have carried out a bifurcation analysis of a highly simplified TB-model, with the only objective, not unexpected given the work in [41,42], in order to highlight the possibility of multiple steady states involving drug-resistant and drug-sensitive TB. The model is not being used to fit any specific outbreak or a particular endemic situation, it is being introduced just to highlight, like in [41,42], that multiple steady states are indeed possible.

7.
Conclusions. Tuberculosis (TB) has colonized 2 billion humans. Fortunately, despite dramatic increases in population density, the growth of mega-cities, intense population mobility and higher contact rates (mass public and air transportation), the likelihood of developing active-TB decreased significantly before the introduction of antibiotics in Europe and the USA. Today, we see "only" about 9 million new cases of active TB per year, throughout the world, that result in roughly 3 million deaths per year. Yet, the fact remains, that one out of three humans is the host of a TB strain. A fact, that should not be ignored as the potential for massive re-activation cannot be totally ruled out. And so, solid understanding of the processes of re-infection and/or re-activation as well as the population-level consequences of such mechanisms for reactivation are critically important. The HIV epidemics increased the rates of TB re-activation ( [32]) for some time and so, we must be prepared. In addition, the lack of compliance with multi-drug TB treatment protocols has also been contributing to the increases in the number of cases of DR-and XDR-TB.
Models have shown that TB exogenous re-infection can support multiple endemic states and bi-stability ( [41,42], a situation that we have shown, in this paper to be possible, in the presence of drug-sensitive and drug-resistant TB strains. It should be intuitively clear, after reading Section 6 that the impact of exogenous reinfection is likely to be altered by changes in HIV prevalence, deteriorating economic conditions (famine, wars, climate change) and the growth of DR-and XDR-TB cases. In this paper, we have just shown that it is dynamically possible in a highly simplified model, without any attempts to connect our results to any specific realistic system, that exogenous re-infection is indeed worth studying.
Looking at policies that bring R 0 to a value less than one may therefore not be sufficient if reactivation or re-infection are strong. In this manuscript, we have looked via rather simple models that ignore epidemiological, social and socioeconomic factors, host's heterogeneity, population structure and "mobility," at the possibility of multiple steady states and bi-stability. In our analysis, we did not incorporate HIV co-infections, or DR-TB, or XDR-TB. The addition of such levels of complexity would more likely strengthen the possibility that multiple steady states are supported by TB dynamics. Reactivation or re-infection mechanisms are therefore extremely likely to play a central role, in the context of complex social environments in highly heterogeneous communities, on TB dynamics.
Prior studies on the role of long and variable periods of latency have been conducted in general [10,12,16,17,22] and so, those interested, should look at their approaches. Here, we also provide specific formulae and a detailed careful analysis in situations when drug-sensitive and drug-resistant TB strains are present. Formulae are provided when a gamma or generalized gamma distributions of sojourn times in the latent class are considered. The formulae are useful enough to those interested in connecting models to data and those that must focus on the study of the transient dynamics of epidemic models-a problem of particular interest given that antibiotic resistance appears to be invading an increasing number of communities. Success due to the relatively easy access that some countries provide to antibiotics and a lack of understanding of the problems generated when those undergoing treatment do not comply rigorously with treatment protocols has been difficult to reach and so, the question remains: How long will we have before drug resistance becomes pervasive?
The models and analysis presented here may help practitioners and public health officials address the above issues via easily carried out simulations. The last section focused on the impact of incorporating of exogenous infection within a simplified two strain model that includes sensitive and drug resistants strains; with drug resistant arising due to a lack of compliance with treatment has the sole purpose of highlighting the possibility of multiple steady states in our multi-strain model. TB-dynamical outcomes are enhanced when exogenous reinfection is included and it seems more so in a two-strain model. It is not surprising to see that the incorporation of exogenous reinfection is not yet a popular alternative among those who believe that R 0 is the key to public policy formulation. Although, the use of R 0 has been extraordinarily useful, providing overall quantitative perspectives on what may be needed to ameliorate the impact of a disease, the fact remains that the preponderance of its use under the limited assumptions of homogenous mixing, needs to be study in deeper ways. Models that incorporate gender differences, variability in susceptibility, age structure, socio-economic factors, mobility and behavioral responses are being systematically considered by computational epidemiologists. The conclusions tend to be that while R 0 is still extremely useful, it is certainly not sufficient for the development of public health policy that accounts for the role of initial conditions [8]. In fact, here, it is shown that exogenous re-infection makes life difficult, particularly to those involved in making health policy decisions. The analysis shows that in general initial conditions do matter, that is, variations in initial conditions, may in fact lead to quite distinct epidemic outbreaks. where, , and, Thus, J is a non-negative matrix since Hence, it follows from (50) that Since R C = max{R C1 , R C2 } ≤ 1 (or, equivalently, the eigenvalues of the matrix F − V all have negative real parts), it follows that the linearized differential inequality system (51) is stable whenever R C1 < 1 and R C2 < 1. Thus, it follows, by Comparison Theorem [24], that lim t→∞ (L 1 (t), L 2 (t), I(t), M (t), J(t)) → (0, 0, 0, 0, 0). (6) showes that S(t) → S * as t → ∞ (for R C1 < 1 and R C2 < 1). Therefore, Hence, the DFE (ε 0 ) of the model (6) is g.a.s. in D whenever R C1 < 1 and R C2 < 1.

B. Proof of Theorem 4.4.
Proof. In this case, M = J = 0 in system (6), that is, The components of the unique endemic equilibrium ε 1 can then be obtained by substituting the unique value of λ * * 1 , given in (14), into the expressions in (12) as following , The Jacobian matrix of model system (52) at ε 1 is given by When R C1 > 1, hence, J (ε T 1 ) is a strictly column diagonally dominant matrix. Also all diagonal elements of J (ε T 1 ) are negative. Therefore, using the Gershgorin circle theorem, the radius of the disc, R i = j =i |a ij |, can easily be defined as R j = j =i |a ji | 776 A. JABBARI, C. CASTILLO-CHAVEZ, F. NAZARI, B. SONG AND H. KHEIRI observing that the TB-strain only equilibrium ε 1 is l.a.s whenever R C1 > 1 > R C2 .
C. Proof of Theorem 4.6.
Proof. In this case, L 1 = L 2 = I = 0 in system (6), that is, The components of the unique endemic equilibrium ε 2 can then be obtained by substituting the unique value of λ * * 2 , given in (18), into the expressions in (16) as following The Jacobian matrix of model system (54) at ε 2 is given by When R C2 > 1, hence, J (ε T 2 ) is a strictly column diagonally dominant matrix. Also all diagonal elements of J (ε T 2 ) are negative, so, using the Gershgorin circle theorem we observe that the drug-resistant TB-strain only equilibrium ε 2 is l.a.s. whenever R C2 > 1 > R C1 .
D. Proof of Theorem 4.10.
Proof. The components of the unique endemic equilibrium ε 3 can be obtained by substituting the unique value of λ * * 1 and λ * * 2 given in (22) and (24), into the expressions in (20) The Jacobian matrix at ε 3 is given by with, .
hence, J (ε T 3 ) is a strictly column diagonally dominant matrix. Also all diagonal elements of J (ε T 3 ) are negative, so, using the Gershgorin circle theorem we observe that interior equilibrium ε 3 is l.a.s. whenever R C1 > R C2 > 1.

E. Proof of Theorem 5.2.
Proof. It should, first of all, be mentioned that the system (5) satisfies the Type K condition [35] (hence, Comparison Theorem can be used [24]). where, The matrices F and V , defined in (26) and (27), respectively.

F. Proof of Theorem 5.4.
Proof. ε 1 exists and is unique if and only if R C1 > 1 and ε 1 to exist alone, it is necessary that the resistant strain does not exist (i.e., R C2 < 1). The components of the unique endemic equilibrium ε 1 can then be obtained by substituting the unique value of λ * * 1 , given in (30), into the expressions in (28) as following , · · · L * * n = κ 1 (R C1 − 1)Λ β 1 a n σ , The Jacobian matrix at ε 1 is given by When R C1 > 1, hence, J (ε T 1 ) is a strictly column diagonally dominant matrix. Also all diagonal elements of J (ε T 1 ) are negative, so, using the Gershgorin circle theorem we observe that the the TB-strain only equilibrium ε 1 is l.a.s. whenever R C1 > 1 > R 2 .
G. Proof of Theorem 5.6.
Proof. ε 2 to exist alone, it is necessary that the resistant strain does not exist (i.e., R C1 < 1). The components of the unique endemic equilibrium ε 2 can then be obtained by substituting the unique value of λ * * 2 , given in (33), into the expressions in (31) as following The Jacobian matrix at ε 2 is given by When R C2 > 1, hence, J (ε T 2 ) is a strictly column diagonally dominant matrix. Also all diagonal elements of J (ε T 2 ) are negative, so, using the Gershgorin circle theorem we observe that the the drug resistant TB-strain-only equilibrium ε 2 is l.a.s whenever R C2 > 1 > R C1 .
H. Proof of Theorem 5.10.
hence, J (ε T 3 ) is a strictly column diagonally dominant matrix. Also all diagonal elements of J (ε T 3 ) are negative, so, using the Gershgorin circle theorem we observe that the the interior equilibrium ε 3 is l.a.s. whenever R C1 > R C2 > 1.   Table 3. Values of a i and b i (i = 1, 2, 3, 4) for various number of disease stages (m and n).      Table  3. Parameter values used are as given in Table 3.   Table 3.