# Modeling the Spread of Tuberculosis in Semiclosed Communities

Authors: Mauricio Herrera; Paul Bosch; Manuel Nájera; Ximena Aguilera
Journal: Computational and Mathematical Methods in Medicine
 (2013)

---

## Abstract

We address the problem of long-term dynamics of tuberculosis (TB) and latent tuberculosis (LTB) in semiclosed communities. These communities are congregate settings with the potential for sustained daily contact for weeks, months, and even years between their members. Basic examples of these communities are prisons, but certain urban/rural communities, some schools, among others could possibly fit well into this definition. These communities present a sort of ideal conditions for TB spread. In order to describe key relevant dynamics of the disease in these communities, we consider a five compartments SEIR model with five possible routes toward TB infection: primary infection after a contact with infected and infectious individuals (fast TB), endogenous reactivation after a period of latency (slow TB), relapse by natural causes after a cure, exogenous reinfection of latently infected, and exogenous reinfection of recovered individuals. We discuss the possible existence of multiple endemic equilibrium states and the role that the two types of exogenous reinfections in the long-term dynamics of the disease could play.


---

## Body

## 1. Introduction

Dynamics of tuberculosis (TB) spread has been the subject of a considerable body of theoretical and mathematical work. For review, see, for example, [1, 2] and references therein. The choice of a particular model is strongly connected to the questions we want to answer, and in the present work we will address the problem of long-term dynamics of tuberculosis and latent tuberculosis (LTB) in semiclosed communities.For semiclosed communities we mean not strictly closed communities with certain mobility of their members out or into the community, with a recruitment of new members and departure of others. But, essentially these communities are congregate settings with the potential for sustained daily contact for weeks, months, and even years between community members. Basic examples of these communities are prisons, but certain urban/rural communities, schools, among others could possibly fit well into this general definition. These communities present a sort of ideal conditions for frequent TB outbreaks, enhanced TB transmission, and accelerated spread of the disease.The basic characteristics of such settings including the possibility of high concentrations of infectious individuals and immunodeficient hosts, improper precautions taken for protection, delay in diagnosis, sustained contact with the index case, and inadequate ventilation and/or overcrowding make them well suited for TB transmission, creating this way genuine high transmission pockets of TB inserted in the general population [3, 4].In fact, prisons are especially high burden communities, in which incidence and prevalence of TB are very high, and consequently the frequency of infections and reinfections considerably increases in comparison with population at large; see the works by Chiang and Riley [5] and by Baussano et al. [6].Studying the dynamics of the TB spread in semiclosed communities is an interesting and significant topic by itself; however, there is an important phenomenon due to which the study of these types of communities is essential in the context of TB spread. This phenomenon has been called theReservoir Effect [6, 7]. Indeed, semiclosed communities such as prisons represent a reservoir for disease transmission to the population at large and should be a source of public concern. For example, TB infection may spread into the general population through prison staff, visitors, and close contacts of released prisoners. The transmission dynamics between prisoners and the general population [6], together with immigration from developing countries with high prevalence of TB [8, 9], has been hypothesized to play a key role in driving overall population-level TB incidence, prevalence, and mortality rates.In a recent work [4] the authors have even gone further in relation to this effect and have named these communities Institutional Amplifiers of TB Propagation. Some examples of communities given by these authors are poor hospitals in which dozens of patients share poorly ventilated communal rooms, crowded prison cell blocks, and mining barracks among others.The transmission and progression of TB infection has been relatively well understood on a population scale. Generally, it is assumed that once an individual is infected with TB, he or she is immune from further infection events. Moreover, it was proposed what came to be known as theunitary concept of pathogenesis [10], which states that TB always begins with primary infection, and subsequent episodes of active TB are due to reactivation of dormant bacilli from this primary infection. However, a persistent evidence has recently been shown (see [5] for a review) that the paths to TB infection are not as linear as was suggested by the unitary concept of pathogenesis. The availability of individual, strain-specific infection histories (see, e.g., [11–13]) has made it clear that exogenous reinfection in people with previously documented TB infection does occur. The important question is whether reinfection occurs commonly enough to have an effect on the overall infection dynamics of the population [14].The relative importance of these pathways to the development of active disease has significant implications for treatment and control strategies, most notably in deciding whether latently infected and treated individuals are at risk of reinfection [15].Several authors [15–20] have declared that exogenous reinfection plays an important role in the disease progression and that the inhalation of tubercle bacilli by persons who have had a primary TB infection previously for more than five years represents an increasing risk to develop active TB soon after reinfection. A study from South Africa [21] has demonstrated that the rate of reinfection by TB after successful treatment could be higher than the rate of new TB infections. In this study the reinfection rate after successful treatment was estimated at 2.2 per 100 person-years, which was approximately seven times the crude incidence rate (313 per 100 000 population per year) and approximately four times the age-adjusted incidence rate of new TB (515 per 100 000 population per year). So, ignoring exogenous reinfection when modeling TB spread in high-incidence and high-prevalence community setting such as semiclosed communities has been seen to be inappropriate. (Henao-Tamayo et al. in [22] recently published a mouse model of TB reinfection that could help to explain immunological aspects of reinfection risk in high-incidence areas.)We will use an SEIR standard compartmental model; see for example the works by Blower et al. [23] and more recently by Liao et al. [24] with some modifications explained bellow that turn out to be quite useful in the study of the particularities of TB spread at this type of communities. This model assumes that the population in the community is homogeneous that it does not consider the heterogeneities in the social structure between community members, and it is based on the so-called mass action or fully mixing approximation. This means that individuals with whom a susceptible individual has contact are chosen at random from the whole community. It is also assumed that all individuals have approximately the same number of contacts in the same time and that all contacts transmit the disease with the same probability.The model we use in this work takes into account the following relevant facts in the context of semiclosed communities.(1)
The overcrowding in the community can increase (compared to what occurs in the population at large) the likelihood of exogenous reinfection due to repeated contacts with active infected individuals. That is, besides primary infection the model considers the possible reinfection of individuals with LTB (individuals who are assumed to be asymptomatic and noninfectious but capable of progressing to active TB) and recovered individuals (individuals who have been treated for TB in the past and been declared cured). If latently infected or recovered individuals remain in the community, they could be infected again.(2)
At present, it is not completely clear whether in all cases previous infections withMycobacterium TB with or without subsequent recovery offer some protection that could be translated into a reduced susceptibility to reinfection [5, 21, 22, 25]. So, we will be open at exploring different situations with regard to this fact in the model.(3)
Poor nutrition, immunodepression,  and other diseases increase the likelihood of accelerated progression to active TB.We will see that considering exogenous reinfection to describe TB spread produces a richer and more complex dynamics than the one observed in previous models (see e.g., [23, 25, 26]). In particular, unlike the model published by Feng et al. in [26], which uses a single parameter for exogenous reinfection, our model uses two parameters related to two possible reinfections (reinfection of latently infected and reinfection of recovered individuals).


## 2. Basic Epidemiology of TB Sources and Probability of Infection in Semiclosed Communities

The risk of infection withMycobacterium tuberculosis, the bacterium causing TB, depends mainly on two factors: first, significant exposure to a source of infection and second, the probability of getting infection if there is exposure.TB is mostly transmitted through the air; tubercle bacilli, that depends on host and agent factors, is distributed in tiny liquid droplets that are produced when someone with clinical or active TB coughs, sneezes, spits, or speaks, allowing infected individual to infect others. In closed places the bacteria are expelled into a finite volume of air unless there is ventilation, see [27]. In these conditions they may remain viable and suspended in the air for a prolonged period of time. But, the number of bacilli excreted by most persons with active pulmonary TB is relatively small [16], so the probability of TB transmission per contact, per unit of time is in general quite low. The risk of infection is very small during a single encounter with an infectious individual [28]. However, the probability of TB transmission can be enhanced by systematic and long exposure of susceptible individuals to particular infectious individuals.The risk of TB transmission is particularly high in settings with poorly ventilated areas (places with reduced air volume per occupant, with ventilation systems which recirculate the air, or with poorly filtered air exchanges) and/or closed areas in which people are in close and frequent contact. Closed regime prisons are examples of these high-risk areas. In effect, the occurrence of TB in prisons for example is usually reported to be much higher than the average levels reported for the corresponding general population [6].Although most exposed individuals develop an effective immune response to the initial infection [17], there is another factor that raises the chances of TB contagion, the fact that TB is an opportunistic disease. Indeed, infected individuals with weakened immune systems are at significant risk of developing clinical TB disease (active TB). High TB prevalence is therefore observed in individuals with HIV infection, poor nutritional status, alcoholism, drug abuse, concurrence of other pathology, and psychological stress decrease immune response levels. These conditions occur frequently in imprisoned peoples.TB is usually described as a slow disease because of its long and variable period of latency and because of its short and relatively narrow infectious period distribution. Long periods of latency (inactive TB or latent TB or LTB) imply that new cases of infection are not clinically noticeable and therefore remain unobserved for a period of time. Immune response of susceptible individuals can restrict proliferation of the bacilli leading to what seems to be long-lasting partial immunity against reinfection or a response capable of stopping the progression from LTB to active TB.Exposed individuals may remain in the latent stage for long and variable periods of time. In fact, it often happens that the host dies without ever developing active TB. The progression from latent to active TB is uncommon in the population at large. It is estimated that only about 5 to 10 percent of LTB individuals develop clinical or active TB [16], but due to the above described extreme conditions at semiclosed communities such as prisons, persons lived in these communities may be at risk of rapid progression from LTB to active TB following recent infection or reactivation of latent infection, or reinfection, see [6].Some additional known epidemiological facts to be considered for TB disease are the following.(1)
Most of the secondary infections generated by an infected individual do take place within the first months following TB activation [29].(2)
In the work by Styblo [16] it was noted that nearly 60 percent of the new cases arose during the first year following infection, while the cumulative number of cases generated over the first five years after infection accounted for nearly 95 percent of the total observed cases. People ill with TB can infect up to 10–15 other people through close contact over the course of a year [30].(3)
Case fatality among untreated pulmonary TB cases is around 66.6 percent [30].(4)
Recovered individuals, naturally or from treatment, may develop active TB again, a phenomenon known as TB relapse. (Recurrent cases (formerly relapse cases) have been treated for TB in the past and been declared successfully treated (cured/treatment completed) at the end of their treatment regimen. Recurrent cases include relapses due to the sameMycobacterium tuberculosis strain as for the previous episode as well as new episodes of TB due to reinfection.)(5)
Individuals with LTB may progress to active TB due to reexposure and reinfection. The extent to which latent tuberculosis infection could reduce the risk of progressive disease following reinfection is not known [31].


## 3. A Compartmental Model for the TB Spread

In order to describe key relevant dynamics in the study of the TB spread in semiclosed communities, we consider five compartments SEIR model represented in Figure1.Figure 1
Flow chart of TB compartmental model.The compartments are uninfected individuals (susceptible), theS class; the latent class E, that is, individuals who are assumed to be asymptomatic and noninfectious but capable of progressing to the clinical disease or active TB; the infectious class I is subdivided into two subclasses: (a) infected and infectious individuals II and (b) infected and noninfectious individuals IN; and the R class of recovered by treatment, self cure, or quarantine.Every individual in theE, II, and IN classes is considered infected. There are five possible routes toward TB infection according to this model: primary infection after a contact with infected and infectious individuals (fast TB), endogenous reactivation after a period of latency (slow TB), relapse by natural causes after a cure, exogenous reinfection of latently infected, and exogenous reinfection of recovered individuals.Thef and q are probability of developing infectious TB if one develops fast and slow TB, respectively, 2w is the relapse rate to active TB. Uninfected individuals are recruited at the rate Π, and μ is the natural mortality rate. Individuals with TB experience a death rate μT due to TB infection.After infection, a fractionp of individuals progresses to disease relatively soon (i.e., within the first two years) after infection; the remaining fraction 1-p of infected individuals become latently infected. Of newly infected individuals who thus progress quickly to disease, f represents the fraction that develops infectious disease, and 1-f represents the fraction that develops noninfectious disease.TheE class, latently infected individuals, does not shed bacilli and is not infective to others. In some latently infected individuals, the infection remains latent and it may persist for life. But, in a small minority of latently infected individuals, reactivation of the latent infection leads to the disease. The coefficients r1 and r2 denote the treatment rates for infected and infectious individuals II class and infected and noninfectious individuals IN class, respectively. The model does not consider unsuccessful treatments.The parameterβ is the primary TB transmission rate; this parameter summarizes socioeconomic and environmental factors that affect primary TB transmission. We assume that transmission rates are determined by broad demographic and social contexts, as well as by characteristics of both the transmitter and recipient (i.e., the number, viability, and virulence of the organisms within sputum droplet nuclei, immune status of the recipient, etc.)TB transmission rate in case of reinfection might be different than the transmission rate of primary infection. The quantities that take into account these differences in case of reinfection of latently infected individuals and reinfection of recovered individuals are given by the dimensionless parametersδ and η, respectively. The parameter δ is the proportion in TB transmission due to exogenous reinfection of latently infected individuals, and η is the proportion in TB transmission due to exogenous reinfection of recovered individuals. Thus, δβ is the exogenous reinfection rate of latently infected, and ηβ is the exogenous reinfection rate of recovered individuals.A conservative point of view will consider thatbiologically plausible values for the reinfection parameters δ and η are given within the intervals 0≤δ≤1, 0≤η≤1. In this case, the parameters δ and η can be interpreted as factors reducing the risk of reinfection of an individual who has previously been infected and has acquired some degree of protective immunity. However, studies on genetic predisposition [22] or in communities with cases as those reported in [21] have gathered some evidence that in certain situations there may be some increased susceptibility to reinfection. Therefore, we are willing to explore in the next sections other mathematical possibilities where the reinfection parameters can take even less usual values δ>1 and η>1.However, recurrent TB due to endogenous reactivation (relapse) and exogenous reinfection could be clinically indistinguishable [32]; they are independent events. For this reason, beside primary infection we will include in the model the possibility of endogenous reactivation and exogenous reinfection as different way toward infection. So, we have the following.(1)
TB due to the endogenous reactivation of primary infection (exacerbation of an old infection) is considered in the model by the termsqνE and (1-q)νE.(2)
TB due to reactivation of primary infection induced by exogenous reinfection is considered by the termsδqβEII and δ(1-q)βEII.(3)
Recurrent TB due to exogenous reinfection after a cure or treatment is described by the termηβIIR.The parameters of the model, its descriptions, and its units are given in Table1.Table 1
Parameters of the model, its descriptions, and its units.
Parameter
Description
Unit
β
Transmission rate
1/year
Π
Recruitment rate
1/year
c
Natural cure rate
1/year
ν
Progression rate from latent TB to active TB
1/year
μ
Natural mortality rate
1/year
μ
T
Mortality rate or fatality rate due to TB
1/year
w
Relapse rate
1/year
q
Probability to develop TB (slow case)
—
f
Probability to develop TB (fast case)
—
p
Proportion of new infections that produce active TB
—
δ
β
Exogenous reinfection rate of latent
1/year
η
β
Exogenous reinfection rate of recovered
1/year
r
1
Treatment rates forII
1/year
r
2
Treatment rates forIN
1/yearAll these considerations give us the following system of equations:(1)dSdt=Π-βSII-μS,dEdt=(1-p)βSII+ηβRII-(ν+μ)E-δβEII,dIIdt=fpβSII+qνE+wR-(μ+μT+c+r1)II+δqβEII,dINdt=(1-f)pβSII+(1-q)νE+wR-(μ+μT+c+r2)IN+δ(1-q)βIIE,dRdt=c(II+IN)-(2w+μ)R-ηβRII+r1II+r2IN.Adding all the equations in (1) together, we have
(2)dNdt=-μN-μT(II+IN)+Π,
where N=S+E+II+IN+R represents the total number of the population, and the region
(3)D={(S,E,II,IN,R)∈ℝ+5:S+E+II+IN+R≤Πμ}
is positively invariant of system (1).It is a common practice in epidemic research to introduce the basic reproduction numberR0, defined as the average number of secondary infections produced by an infected individual in a completely susceptible population, as the measure for the epidemic thresholds, if R0>1 an epidemic will arise.We have calculatedR0 for this model using the Next Generation Method [35] and it is given by
(4)R0=βΠ((hfp+(1-p)νq)(ab-mw)R0=βΠj+mw(h(1-f)p+(1-p)ν(1-q)))×(μha(ab-gw-mw))-1,
where
(5)a=μ+μT+c,b=2w+μ,h=ν+μ,g=r1+c,m=r2+c.


### 3.1. Steady-State Solutions

In order to find steady-state solutions for (1) we have to solve the following system of equations:
(6)0=Π-βSII-μS,0=(1-p)βSII+ηβRII-(ν+μ)E-δβEII,0=fpβSII+qνE+wR-(μ+μT+c+r1)II+δqβEII,0=(1-f)pβSII+(1-q)νE+wR-(μ+μT+c+r2)IN+δ(1-q)βIIE,0=c(II+IN)-(2w+μ)R-ηβRII+r1II+r2IN.Solving system (6) with respect to II we have the following equation:
(7)II(AII3+BII2+CII+D)=0.The coefficients of (7) are all expressed as functions of the parameters listed in Table 1. However, these expressions are too long to be written here. See Appendix A for explicit forms of the coefficients.


#### 3.1.1. Disease-Free Equilibrium

ForII=0 we get the disease-free steady-state solution:
(8)P0=(S0,E0,II0,IN0,R0)=(Πμ,0,0,0,0).Lemma 1.
The disease-free steady-state solutionP0 is locally asymptotic stable for R0<1.Proof.
The Jacobian matrix of system (6) evaluated in P0=((Π/μ),0,0,0,0) is
(9)J=[-μ0-βΠμ000-ν-μ(1-p)βΠμ000qνfpβΠμ-μ-μT-c-r10w0(1-q)ν(1-f)pβΠμ-μ-μT-c-r2w00r1+cr2+c-2w-μ].
The characteristic equation for J have the form
(10)(λ+μ)(λ4+a3λ3+a2λ2+a1λ+a0)=0.
Given the polynomial
(11)P(λ)=λ4+a3λ3+a2λ2+a1λ+a0=0,
in the special case when a1,a2,a3>0, 3 roots of the polynomial P(λ) have negative real part and if(i)
a0=0, the 4th root, or largest eigenvalue, is zero,(ii)
a0>0, all eigenvalues are negative,(iii)
a0<0, the largest eigenvalue has positive real part.
Thus, the stability of disease-free steady-state solution is determined solely by the sign of the constant terma0 of the polynomial P(λ) [36].
The coefficientsa0,a1,a2,a3 are all decreasing functions of β and they are linear functions with respect to the parameter β, so they all take the general form ai(β)=-Aiβ+Bi with i=0,1,2,3 and Ai,Bi>0. We can define βi=Bi/Ai. It is easy to see that for β<βi we have ai(β)>0.
For example,(12)a3(β)=-fpβΠμ+ν+3μ+μT+c+r2+2w,β3=(ν+3μ+μT+c+r2+2w)μfpΠ.
By straightforward calculations and reminding that the coefficientsai are decreasing functions of β, we found that
(13)a2(β3)<0⟹β2<β3,a1(β2)<0⟹β1<β2,a0(β1)<0⟹β0<β1.
This way,β0<β1<β2<β3 and if we take β<β0, all the coefficients ai are positive. But, for R0<1, we can see that the condition β<β0 is fulfilled. Indeed, the constant term a0 of the polynomial P(λ) can be written as
(14)a0(β)=-A0β+B0=B0(1-A0βB0)=B0(1-R0).
Using this form for the coefficienta0 we can see that if R0<1, then a0(β)>0 so β<β0.Remark 2.
ForR0>1 we have a0<0, and the disease-free steady-state solution is unstable. Indeed, if λ1,λ2,λ3, and λ4 are the roots of polynomial P(λ)=0, we have λ1λ2λ3λ4=a0<0, and it is impossible that the four roots have negative real part.


#### 3.1.2. Endemic Equilibria

From (7) we get that nontrivial solutions are possible if
(15)AII3+BII2+CII+D=0.The explicit expressions for the coefficientsA and D are
(16)A=β3δσ(μ+μT)(μ+μT+c+(1-q)r1+qr2)>0,(17)D=h((r1+a)[ab+wr2+r1(w+μ)](μ+μt)+μm(r1+a))(1-R0),
where parameters a, b, h, and m are defined as in (5).The coefficientsB and C can be written in the following general form:
(18)B=β2fB(β),C=βfC(β),
where
(19)fB(β)=-B1β+B2,fC(β)=-C1β+C2.The coefficients{B}i=1,2 and {C}i=1,2 are all positive and depend on the parameters given in Table 1. See Appendix A for the explicit form of these coefficients.Changes in the signs of the coefficientB and C as function of transmission rate β can be explained using the above defined functions fB(β) and fC(β), respectively. The functions fB(β) and fC(β) both are linear and decreasing functions of β.Consider the polynomial function(20)P(x)=Ax3+Bx2+Cx+D.From (17) we can see that for R0>1 the coefficient D is negative, so we have P(0)=D<0. On the other hand, because the coefficient A is always positive, there must be a value x* such that, for x>x*, it holds that P(x)>0. Since function P(x) is continuous, this implies the existence of solution for the equation P(x)=0.To determine how many possible endemic states arise, we consider the derivativeP′(x)=3Ax2+2Bx+C, and then we analyse the following cases.(1)
IfΔ=B2-3AC≤0, P′(x)≥0 for all x, then P(x) is monotonically increasing function and we have a unique solution, that is, a unique endemic equilibrium.(2)
IfΔ≥0, we have solutions of the equation P′(x)=0 given by
(21)x2,1=-B±B2-3AC3Aand P′(x)≥0 for all x≥x2 and x≤x1. So, we need to consider the positions of the roots x1 and x2 in the real line. We have the following possible cases.(i)
IfC≤0, then for both cases B≥0 and B<0, we have x1<0, x2>0 and P′(x)>0 for all x≥x2≥0. Given that P(0)=D<0, this implies the existence of a unique endemic equilibrium.(ii)
IfB≥0 and C≥0, then both roots x1 and x2 are negative and P′(x)>0 for all x≥0.(iii)
IfB<0 and C>0, then both roots x1 and x2 are positive and we have the possibility of multiple endemic equilibria. This is a necessary condition, but not sufficient. It must be fulfilled also that P(x1)≥0.LetβB be the value of β such that fB(βB)=0 and βC the value of β such that fC(βC)=0. Moreover, let βR0 be the value for which the basic reproduction number R0 is equal to one (the value of β such that coefficient D becomes zero).Lemma 3.
If the conditionβR0<βC<βB is met, then system (1) has a unique endemic equilibrium for all β>βR0 (Table 3).Proof.
Using similar arguments to those used in the proof of Lemma 1,  we have, given the conditionβR0<βC<βB, that for all values of β such that β<βR0, all polynomial coefficients are positive; therefore, all solutions of the polynomial are negative and there is no endemic equilibrium (positive epidemiologically meaningful solution).
ForβR0<β<βC the coefficients B and C are both positive, while the coefficient D is negative; therefore, appears only one positive solution of the polynomial (the greatest one), so we have a unique endemic equilibrium.
ForβC<β<βB, the coefficient C is negative and B is positive. According to the cases studied above we have in this situation a unique endemic equilibrium.
Finally, forβ>βB the coefficients B and C are both negative, and according to the study of cases given above we also have a unique positive solution or endemic equilibrium.Let us first consider biologically plausible values for the reinfection parametersη and δ, that is, values within the intervals 0≤δ≤1, 0≤η≤1. This means that the likelihood of both variants of reinfections is no greater than the likelihood of primary TB. So, we are considering here partial immunity after a primary TB infection.Lemma 4.
For biologically plausible values(δ,η)∈[0,1]×[0,1] system (1) fulfils the condition βR0<βC<βB.Proof.
Using straightforward but cumbersome calculations (we use a symbolic software for this task), we were able to prove that if we consider all parameters positive (as it is the case) and taking into account biologically plausible values(δ,η)∈[0,1]×[0,1], then fB(βC)>0 and D(βB)>0 and it is easy to see that these inequalities are equivalent to βR0<βC<βB.We have proven that the conditionβR0<βC<βB implies that the system can only realize two epidemiologically meaningful (nonnegative) equilibrium states. Indeed, if we consider the disease transmission rate β as a bifurcation parameter for (1), then we can see that the system experiences a transcritical bifurcation at β=βR0, that is, when R0=1 (see Figure 2). If the condition βR0<βC<βB is met, the system has a single steady-state solution, corresponding to zero prevalence and elimination of the TB epidemic for β<βR0, that is, R0<1, and two equilibrium states corresponding to endemic TB and zero prevalence when β>βR0, that is, R0>1. Moreover, according to Lemma 4 this condition is fulfilled in the biologically plausible domain for exogenous reinfection parameters (δ,η)∈[0,1]×[0,1]. This case is summarized in Table 2.Table 2
Qualitative behaviour for system (1) as a function of the disease transmission rate β, when the condition βR0<βC<βB is fulfilled.
Interval
Coefficients
Type of equilibrium
β
<
β
R
0
A
>
0, B>0, C>0, D>0
Disease-free equilibrium
β
R
0
<
β
<
β
C
A
>
0, B>0, C>0, D<0
Unique endemic equilibrium
β
C
<
β
<
β
B
A
>
0, B>0, C<0, D<0
Unique endemic equilibrium
β
>
β
B
A
>
0, B<0, C<0, D<0
Unique endemic equilibriumTable 3
Qualitative behaviour for system (1) as function of the disease transmission rate β, when the condition βB<βC<βR0 is fulfilled. Here, Δ1 is the discriminant of the cubic polynomial (20).
Interval
Coefficients
Type of equilibrium
β
<
β
B
A
>
0, B>0, C>0, D>0
Disease-free equilibrium
β
B
<
β
<
β
C
A
>
0, B<0, C>0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
C
<
β
<
β
R
0
A
>
0, B<0, C<0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
R
0
<
β
A
>
0, B<0, C<0, D<0
Unique endemic equilibriumFigure 2
Bifurcation diagram (solutionx of polynomial (20) versus β) for the condition βR0<βC<βB. βR0 is the bifurcation value. The blue branch in the graph is a stable endemic equilibrium which appears for R0>1.From Table2 we can see that although the signs of the polynomial coefficients may change, other new biologically meaningful solutions (nonnegative solutions) do not arise in this case. The system can only display the presence of two equilibrium states: disease-free or a unique endemic equilibrium.The basic reproduction numberR0 in this case explains well the appearance of the transcritical bifurcation, that is, when a unique endemic state arises and the disease-free equilibrium becomes unstable (see blue line in Figure 2).However, the change in signs of the polynomial coefficients modifies the qualitative type of the equilibria. This fact is shown in Figures5 and 7 illustrating the existence of focus or node type steady-sate solutions. These different types of equilibria as we will see in the next section cannot be explained using solely the reproduction number R0.In the next section we will explore numerically the parametric space of system (1), looking for different qualitative dynamics of TB epidemics. We will discuss in more detail how dynamics depends on the parameters given in Table 1, especially on the transmission rate β, which will be used as bifurcation parameter for the model.Let us consider here briefly two examples of parametric regimes for the model in order to illustrate the possibility to encounter a more complex dynamics, which cannot be solely explained by changes in the value of the basic reproduction numberR0.Example I. Suppose β=βR0, this implies that R0=1 and D=0; therefore, we have the equation:
(22)P(x)=Ax3+Bx2+Cx=x(Ax2+Bx2+C)=0.It is easy to see that besides zero solution, ifB<0, C>0 and B2-4AC>0,  (22) has two positive solutions x1 and x2. So, we have in this case three nonnegative equilibria for the system.The conditionB<0 for β=βR0 means fB(βR0)<0, and this in turn implies that βB<βR0. On the other hand, the condition C>0 implies fC(βR0)>0 and therefore βR0<βC.Gathering both inequalities we can conclude that ifβB<βR0<βC, then the system has the possibility of multiple equilibria.Since the coefficientsA and B are both continuous functions of β, we can always find a ϵ neighbourhood of βR0, |β-βR0|<ϵ such that the signs of these coefficients are preserved. Although in this case we do not have the solution x=0, we eventually could still have two positive solutions and consequently, multiple equilibrium states; see the green line in Figure 3.Figure 3
PolynomialP(x) for different values of β with the condition βB<βR0<βC. The graphs were obtained for values of δ=3.0 and η=2.2. The dashed black line indicates the case β=βR0. The figure shows the existence of multiple equilibria.Example II. Suppose we take numerical values for the parameters in Table 1 such that the condition βB<βC<βR0 is fulfilled.Ifβ<βB, then all coefficients of the polynomial (20) are positive and there is not nonnegative solutions. In this case, the system has only a disease-free equilibrium.ForβB<β<βC and βC<β<βR0 the signs of the coefficients of the polynomial are A>0, B<0, C>0, and D>0,A>0, B<0, C<0, D>0, respectively.In both cases the polynomial has two possibilities:(a)
three real solutions: one negative and two positive solutions forΔ1<0,(b)
one negative and two complex conjugate solutions forΔ1>0.HereΔ1 is the discriminant for the polynomial (20).In the (a) case we have the possibility of multiple endemic states for system (1). This case is illustrated in numerical simulations in the next section by Figures 8 and 9.We should note that the valueβ=βB is not a bifurcation value for the parameter β.Ifβ=βB, then A>0, B=0, C>0, and D>0. In this case we have
(23)Δ1=14D2A2+127C3A3>0.The discriminantΔ1 is a continuous function of β, for this reason this sign will be preserved in a ϵ neighbourhood of βB.We should be able to find a bifurcation value solving numerically the equation(24)Δ1(β*)=0,
where β* can be bounded by the interval βB<β*<βR0 (see Figure 4).Figure 4
Bifurcation diagram for the conditionβB<βC<βR0. β* is the bifurcation value. The blue branch in the graph is a stable endemic equilibrium which appears even for R0 < 1.Figure 5
Numerical simulation forR0=3.585422172, δ=0.9, η=0.01, and β=0.00052. The system goes toward a focus type stable stationary equilibrium.


## 3.1. Steady-State Solutions

In order to find steady-state solutions for (1) we have to solve the following system of equations:
(6)0=Π-βSII-μS,0=(1-p)βSII+ηβRII-(ν+μ)E-δβEII,0=fpβSII+qνE+wR-(μ+μT+c+r1)II+δqβEII,0=(1-f)pβSII+(1-q)νE+wR-(μ+μT+c+r2)IN+δ(1-q)βIIE,0=c(II+IN)-(2w+μ)R-ηβRII+r1II+r2IN.Solving system (6) with respect to II we have the following equation:
(7)II(AII3+BII2+CII+D)=0.The coefficients of (7) are all expressed as functions of the parameters listed in Table 1. However, these expressions are too long to be written here. See Appendix A for explicit forms of the coefficients.


### 3.1.1. Disease-Free Equilibrium

ForII=0 we get the disease-free steady-state solution:
(8)P0=(S0,E0,II0,IN0,R0)=(Πμ,0,0,0,0).Lemma 1.
The disease-free steady-state solutionP0 is locally asymptotic stable for R0<1.Proof.
The Jacobian matrix of system (6) evaluated in P0=((Π/μ),0,0,0,0) is
(9)J=[-μ0-βΠμ000-ν-μ(1-p)βΠμ000qνfpβΠμ-μ-μT-c-r10w0(1-q)ν(1-f)pβΠμ-μ-μT-c-r2w00r1+cr2+c-2w-μ].
The characteristic equation for J have the form
(10)(λ+μ)(λ4+a3λ3+a2λ2+a1λ+a0)=0.
Given the polynomial
(11)P(λ)=λ4+a3λ3+a2λ2+a1λ+a0=0,
in the special case when a1,a2,a3>0, 3 roots of the polynomial P(λ) have negative real part and if(i)
a0=0, the 4th root, or largest eigenvalue, is zero,(ii)
a0>0, all eigenvalues are negative,(iii)
a0<0, the largest eigenvalue has positive real part.
Thus, the stability of disease-free steady-state solution is determined solely by the sign of the constant terma0 of the polynomial P(λ) [36].
The coefficientsa0,a1,a2,a3 are all decreasing functions of β and they are linear functions with respect to the parameter β, so they all take the general form ai(β)=-Aiβ+Bi with i=0,1,2,3 and Ai,Bi>0. We can define βi=Bi/Ai. It is easy to see that for β<βi we have ai(β)>0.
For example,(12)a3(β)=-fpβΠμ+ν+3μ+μT+c+r2+2w,β3=(ν+3μ+μT+c+r2+2w)μfpΠ.
By straightforward calculations and reminding that the coefficientsai are decreasing functions of β, we found that
(13)a2(β3)<0⟹β2<β3,a1(β2)<0⟹β1<β2,a0(β1)<0⟹β0<β1.
This way,β0<β1<β2<β3 and if we take β<β0, all the coefficients ai are positive. But, for R0<1, we can see that the condition β<β0 is fulfilled. Indeed, the constant term a0 of the polynomial P(λ) can be written as
(14)a0(β)=-A0β+B0=B0(1-A0βB0)=B0(1-R0).
Using this form for the coefficienta0 we can see that if R0<1, then a0(β)>0 so β<β0.Remark 2.
ForR0>1 we have a0<0, and the disease-free steady-state solution is unstable. Indeed, if λ1,λ2,λ3, and λ4 are the roots of polynomial P(λ)=0, we have λ1λ2λ3λ4=a0<0, and it is impossible that the four roots have negative real part.


### 3.1.2. Endemic Equilibria

From (7) we get that nontrivial solutions are possible if
(15)AII3+BII2+CII+D=0.The explicit expressions for the coefficientsA and D are
(16)A=β3δσ(μ+μT)(μ+μT+c+(1-q)r1+qr2)>0,(17)D=h((r1+a)[ab+wr2+r1(w+μ)](μ+μt)+μm(r1+a))(1-R0),
where parameters a, b, h, and m are defined as in (5).The coefficientsB and C can be written in the following general form:
(18)B=β2fB(β),C=βfC(β),
where
(19)fB(β)=-B1β+B2,fC(β)=-C1β+C2.The coefficients{B}i=1,2 and {C}i=1,2 are all positive and depend on the parameters given in Table 1. See Appendix A for the explicit form of these coefficients.Changes in the signs of the coefficientB and C as function of transmission rate β can be explained using the above defined functions fB(β) and fC(β), respectively. The functions fB(β) and fC(β) both are linear and decreasing functions of β.Consider the polynomial function(20)P(x)=Ax3+Bx2+Cx+D.From (17) we can see that for R0>1 the coefficient D is negative, so we have P(0)=D<0. On the other hand, because the coefficient A is always positive, there must be a value x* such that, for x>x*, it holds that P(x)>0. Since function P(x) is continuous, this implies the existence of solution for the equation P(x)=0.To determine how many possible endemic states arise, we consider the derivativeP′(x)=3Ax2+2Bx+C, and then we analyse the following cases.(1)
IfΔ=B2-3AC≤0, P′(x)≥0 for all x, then P(x) is monotonically increasing function and we have a unique solution, that is, a unique endemic equilibrium.(2)
IfΔ≥0, we have solutions of the equation P′(x)=0 given by
(21)x2,1=-B±B2-3AC3Aand P′(x)≥0 for all x≥x2 and x≤x1. So, we need to consider the positions of the roots x1 and x2 in the real line. We have the following possible cases.(i)
IfC≤0, then for both cases B≥0 and B<0, we have x1<0, x2>0 and P′(x)>0 for all x≥x2≥0. Given that P(0)=D<0, this implies the existence of a unique endemic equilibrium.(ii)
IfB≥0 and C≥0, then both roots x1 and x2 are negative and P′(x)>0 for all x≥0.(iii)
IfB<0 and C>0, then both roots x1 and x2 are positive and we have the possibility of multiple endemic equilibria. This is a necessary condition, but not sufficient. It must be fulfilled also that P(x1)≥0.LetβB be the value of β such that fB(βB)=0 and βC the value of β such that fC(βC)=0. Moreover, let βR0 be the value for which the basic reproduction number R0 is equal to one (the value of β such that coefficient D becomes zero).Lemma 3.
If the conditionβR0<βC<βB is met, then system (1) has a unique endemic equilibrium for all β>βR0 (Table 3).Proof.
Using similar arguments to those used in the proof of Lemma 1,  we have, given the conditionβR0<βC<βB, that for all values of β such that β<βR0, all polynomial coefficients are positive; therefore, all solutions of the polynomial are negative and there is no endemic equilibrium (positive epidemiologically meaningful solution).
ForβR0<β<βC the coefficients B and C are both positive, while the coefficient D is negative; therefore, appears only one positive solution of the polynomial (the greatest one), so we have a unique endemic equilibrium.
ForβC<β<βB, the coefficient C is negative and B is positive. According to the cases studied above we have in this situation a unique endemic equilibrium.
Finally, forβ>βB the coefficients B and C are both negative, and according to the study of cases given above we also have a unique positive solution or endemic equilibrium.Let us first consider biologically plausible values for the reinfection parametersη and δ, that is, values within the intervals 0≤δ≤1, 0≤η≤1. This means that the likelihood of both variants of reinfections is no greater than the likelihood of primary TB. So, we are considering here partial immunity after a primary TB infection.Lemma 4.
For biologically plausible values(δ,η)∈[0,1]×[0,1] system (1) fulfils the condition βR0<βC<βB.Proof.
Using straightforward but cumbersome calculations (we use a symbolic software for this task), we were able to prove that if we consider all parameters positive (as it is the case) and taking into account biologically plausible values(δ,η)∈[0,1]×[0,1], then fB(βC)>0 and D(βB)>0 and it is easy to see that these inequalities are equivalent to βR0<βC<βB.We have proven that the conditionβR0<βC<βB implies that the system can only realize two epidemiologically meaningful (nonnegative) equilibrium states. Indeed, if we consider the disease transmission rate β as a bifurcation parameter for (1), then we can see that the system experiences a transcritical bifurcation at β=βR0, that is, when R0=1 (see Figure 2). If the condition βR0<βC<βB is met, the system has a single steady-state solution, corresponding to zero prevalence and elimination of the TB epidemic for β<βR0, that is, R0<1, and two equilibrium states corresponding to endemic TB and zero prevalence when β>βR0, that is, R0>1. Moreover, according to Lemma 4 this condition is fulfilled in the biologically plausible domain for exogenous reinfection parameters (δ,η)∈[0,1]×[0,1]. This case is summarized in Table 2.Table 2
Qualitative behaviour for system (1) as a function of the disease transmission rate β, when the condition βR0<βC<βB is fulfilled.
Interval
Coefficients
Type of equilibrium
β
<
β
R
0
A
>
0, B>0, C>0, D>0
Disease-free equilibrium
β
R
0
<
β
<
β
C
A
>
0, B>0, C>0, D<0
Unique endemic equilibrium
β
C
<
β
<
β
B
A
>
0, B>0, C<0, D<0
Unique endemic equilibrium
β
>
β
B
A
>
0, B<0, C<0, D<0
Unique endemic equilibriumTable 3
Qualitative behaviour for system (1) as function of the disease transmission rate β, when the condition βB<βC<βR0 is fulfilled. Here, Δ1 is the discriminant of the cubic polynomial (20).
Interval
Coefficients
Type of equilibrium
β
<
β
B
A
>
0, B>0, C>0, D>0
Disease-free equilibrium
β
B
<
β
<
β
C
A
>
0, B<0, C>0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
C
<
β
<
β
R
0
A
>
0, B<0, C<0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
R
0
<
β
A
>
0, B<0, C<0, D<0
Unique endemic equilibriumFigure 2
Bifurcation diagram (solutionx of polynomial (20) versus β) for the condition βR0<βC<βB. βR0 is the bifurcation value. The blue branch in the graph is a stable endemic equilibrium which appears for R0>1.From Table2 we can see that although the signs of the polynomial coefficients may change, other new biologically meaningful solutions (nonnegative solutions) do not arise in this case. The system can only display the presence of two equilibrium states: disease-free or a unique endemic equilibrium.The basic reproduction numberR0 in this case explains well the appearance of the transcritical bifurcation, that is, when a unique endemic state arises and the disease-free equilibrium becomes unstable (see blue line in Figure 2).However, the change in signs of the polynomial coefficients modifies the qualitative type of the equilibria. This fact is shown in Figures5 and 7 illustrating the existence of focus or node type steady-sate solutions. These different types of equilibria as we will see in the next section cannot be explained using solely the reproduction number R0.In the next section we will explore numerically the parametric space of system (1), looking for different qualitative dynamics of TB epidemics. We will discuss in more detail how dynamics depends on the parameters given in Table 1, especially on the transmission rate β, which will be used as bifurcation parameter for the model.Let us consider here briefly two examples of parametric regimes for the model in order to illustrate the possibility to encounter a more complex dynamics, which cannot be solely explained by changes in the value of the basic reproduction numberR0.Example I. Suppose β=βR0, this implies that R0=1 and D=0; therefore, we have the equation:
(22)P(x)=Ax3+Bx2+Cx=x(Ax2+Bx2+C)=0.It is easy to see that besides zero solution, ifB<0, C>0 and B2-4AC>0,  (22) has two positive solutions x1 and x2. So, we have in this case three nonnegative equilibria for the system.The conditionB<0 for β=βR0 means fB(βR0)<0, and this in turn implies that βB<βR0. On the other hand, the condition C>0 implies fC(βR0)>0 and therefore βR0<βC.Gathering both inequalities we can conclude that ifβB<βR0<βC, then the system has the possibility of multiple equilibria.Since the coefficientsA and B are both continuous functions of β, we can always find a ϵ neighbourhood of βR0, |β-βR0|<ϵ such that the signs of these coefficients are preserved. Although in this case we do not have the solution x=0, we eventually could still have two positive solutions and consequently, multiple equilibrium states; see the green line in Figure 3.Figure 3
PolynomialP(x) for different values of β with the condition βB<βR0<βC. The graphs were obtained for values of δ=3.0 and η=2.2. The dashed black line indicates the case β=βR0. The figure shows the existence of multiple equilibria.Example II. Suppose we take numerical values for the parameters in Table 1 such that the condition βB<βC<βR0 is fulfilled.Ifβ<βB, then all coefficients of the polynomial (20) are positive and there is not nonnegative solutions. In this case, the system has only a disease-free equilibrium.ForβB<β<βC and βC<β<βR0 the signs of the coefficients of the polynomial are A>0, B<0, C>0, and D>0,A>0, B<0, C<0, D>0, respectively.In both cases the polynomial has two possibilities:(a)
three real solutions: one negative and two positive solutions forΔ1<0,(b)
one negative and two complex conjugate solutions forΔ1>0.HereΔ1 is the discriminant for the polynomial (20).In the (a) case we have the possibility of multiple endemic states for system (1). This case is illustrated in numerical simulations in the next section by Figures 8 and 9.We should note that the valueβ=βB is not a bifurcation value for the parameter β.Ifβ=βB, then A>0, B=0, C>0, and D>0. In this case we have
(23)Δ1=14D2A2+127C3A3>0.The discriminantΔ1 is a continuous function of β, for this reason this sign will be preserved in a ϵ neighbourhood of βB.We should be able to find a bifurcation value solving numerically the equation(24)Δ1(β*)=0,
where β* can be bounded by the interval βB<β*<βR0 (see Figure 4).Figure 4
Bifurcation diagram for the conditionβB<βC<βR0. β* is the bifurcation value. The blue branch in the graph is a stable endemic equilibrium which appears even for R0 < 1.Figure 5
Numerical simulation forR0=3.585422172, δ=0.9, η=0.01, and β=0.00052. The system goes toward a focus type stable stationary equilibrium.


## 3.1.1. Disease-Free Equilibrium

ForII=0 we get the disease-free steady-state solution:
(8)P0=(S0,E0,II0,IN0,R0)=(Πμ,0,0,0,0).Lemma 1.
The disease-free steady-state solutionP0 is locally asymptotic stable for R0<1.Proof.
The Jacobian matrix of system (6) evaluated in P0=((Π/μ),0,0,0,0) is
(9)J=[-μ0-βΠμ000-ν-μ(1-p)βΠμ000qνfpβΠμ-μ-μT-c-r10w0(1-q)ν(1-f)pβΠμ-μ-μT-c-r2w00r1+cr2+c-2w-μ].
The characteristic equation for J have the form
(10)(λ+μ)(λ4+a3λ3+a2λ2+a1λ+a0)=0.
Given the polynomial
(11)P(λ)=λ4+a3λ3+a2λ2+a1λ+a0=0,
in the special case when a1,a2,a3>0, 3 roots of the polynomial P(λ) have negative real part and if(i)
a0=0, the 4th root, or largest eigenvalue, is zero,(ii)
a0>0, all eigenvalues are negative,(iii)
a0<0, the largest eigenvalue has positive real part.
Thus, the stability of disease-free steady-state solution is determined solely by the sign of the constant terma0 of the polynomial P(λ) [36].
The coefficientsa0,a1,a2,a3 are all decreasing functions of β and they are linear functions with respect to the parameter β, so they all take the general form ai(β)=-Aiβ+Bi with i=0,1,2,3 and Ai,Bi>0. We can define βi=Bi/Ai. It is easy to see that for β<βi we have ai(β)>0.
For example,(12)a3(β)=-fpβΠμ+ν+3μ+μT+c+r2+2w,β3=(ν+3μ+μT+c+r2+2w)μfpΠ.
By straightforward calculations and reminding that the coefficientsai are decreasing functions of β, we found that
(13)a2(β3)<0⟹β2<β3,a1(β2)<0⟹β1<β2,a0(β1)<0⟹β0<β1.
This way,β0<β1<β2<β3 and if we take β<β0, all the coefficients ai are positive. But, for R0<1, we can see that the condition β<β0 is fulfilled. Indeed, the constant term a0 of the polynomial P(λ) can be written as
(14)a0(β)=-A0β+B0=B0(1-A0βB0)=B0(1-R0).
Using this form for the coefficienta0 we can see that if R0<1, then a0(β)>0 so β<β0.Remark 2.
ForR0>1 we have a0<0, and the disease-free steady-state solution is unstable. Indeed, if λ1,λ2,λ3, and λ4 are the roots of polynomial P(λ)=0, we have λ1λ2λ3λ4=a0<0, and it is impossible that the four roots have negative real part.


## 3.1.2. Endemic Equilibria

From (7) we get that nontrivial solutions are possible if
(15)AII3+BII2+CII+D=0.The explicit expressions for the coefficientsA and D are
(16)A=β3δσ(μ+μT)(μ+μT+c+(1-q)r1+qr2)>0,(17)D=h((r1+a)[ab+wr2+r1(w+μ)](μ+μt)+μm(r1+a))(1-R0),
where parameters a, b, h, and m are defined as in (5).The coefficientsB and C can be written in the following general form:
(18)B=β2fB(β),C=βfC(β),
where
(19)fB(β)=-B1β+B2,fC(β)=-C1β+C2.The coefficients{B}i=1,2 and {C}i=1,2 are all positive and depend on the parameters given in Table 1. See Appendix A for the explicit form of these coefficients.Changes in the signs of the coefficientB and C as function of transmission rate β can be explained using the above defined functions fB(β) and fC(β), respectively. The functions fB(β) and fC(β) both are linear and decreasing functions of β.Consider the polynomial function(20)P(x)=Ax3+Bx2+Cx+D.From (17) we can see that for R0>1 the coefficient D is negative, so we have P(0)=D<0. On the other hand, because the coefficient A is always positive, there must be a value x* such that, for x>x*, it holds that P(x)>0. Since function P(x) is continuous, this implies the existence of solution for the equation P(x)=0.To determine how many possible endemic states arise, we consider the derivativeP′(x)=3Ax2+2Bx+C, and then we analyse the following cases.(1)
IfΔ=B2-3AC≤0, P′(x)≥0 for all x, then P(x) is monotonically increasing function and we have a unique solution, that is, a unique endemic equilibrium.(2)
IfΔ≥0, we have solutions of the equation P′(x)=0 given by
(21)x2,1=-B±B2-3AC3Aand P′(x)≥0 for all x≥x2 and x≤x1. So, we need to consider the positions of the roots x1 and x2 in the real line. We have the following possible cases.(i)
IfC≤0, then for both cases B≥0 and B<0, we have x1<0, x2>0 and P′(x)>0 for all x≥x2≥0. Given that P(0)=D<0, this implies the existence of a unique endemic equilibrium.(ii)
IfB≥0 and C≥0, then both roots x1 and x2 are negative and P′(x)>0 for all x≥0.(iii)
IfB<0 and C>0, then both roots x1 and x2 are positive and we have the possibility of multiple endemic equilibria. This is a necessary condition, but not sufficient. It must be fulfilled also that P(x1)≥0.LetβB be the value of β such that fB(βB)=0 and βC the value of β such that fC(βC)=0. Moreover, let βR0 be the value for which the basic reproduction number R0 is equal to one (the value of β such that coefficient D becomes zero).Lemma 3.
If the conditionβR0<βC<βB is met, then system (1) has a unique endemic equilibrium for all β>βR0 (Table 3).Proof.
Using similar arguments to those used in the proof of Lemma 1,  we have, given the conditionβR0<βC<βB, that for all values of β such that β<βR0, all polynomial coefficients are positive; therefore, all solutions of the polynomial are negative and there is no endemic equilibrium (positive epidemiologically meaningful solution).
ForβR0<β<βC the coefficients B and C are both positive, while the coefficient D is negative; therefore, appears only one positive solution of the polynomial (the greatest one), so we have a unique endemic equilibrium.
ForβC<β<βB, the coefficient C is negative and B is positive. According to the cases studied above we have in this situation a unique endemic equilibrium.
Finally, forβ>βB the coefficients B and C are both negative, and according to the study of cases given above we also have a unique positive solution or endemic equilibrium.Let us first consider biologically plausible values for the reinfection parametersη and δ, that is, values within the intervals 0≤δ≤1, 0≤η≤1. This means that the likelihood of both variants of reinfections is no greater than the likelihood of primary TB. So, we are considering here partial immunity after a primary TB infection.Lemma 4.
For biologically plausible values(δ,η)∈[0,1]×[0,1] system (1) fulfils the condition βR0<βC<βB.Proof.
Using straightforward but cumbersome calculations (we use a symbolic software for this task), we were able to prove that if we consider all parameters positive (as it is the case) and taking into account biologically plausible values(δ,η)∈[0,1]×[0,1], then fB(βC)>0 and D(βB)>0 and it is easy to see that these inequalities are equivalent to βR0<βC<βB.We have proven that the conditionβR0<βC<βB implies that the system can only realize two epidemiologically meaningful (nonnegative) equilibrium states. Indeed, if we consider the disease transmission rate β as a bifurcation parameter for (1), then we can see that the system experiences a transcritical bifurcation at β=βR0, that is, when R0=1 (see Figure 2). If the condition βR0<βC<βB is met, the system has a single steady-state solution, corresponding to zero prevalence and elimination of the TB epidemic for β<βR0, that is, R0<1, and two equilibrium states corresponding to endemic TB and zero prevalence when β>βR0, that is, R0>1. Moreover, according to Lemma 4 this condition is fulfilled in the biologically plausible domain for exogenous reinfection parameters (δ,η)∈[0,1]×[0,1]. This case is summarized in Table 2.Table 2
Qualitative behaviour for system (1) as a function of the disease transmission rate β, when the condition βR0<βC<βB is fulfilled.
Interval
Coefficients
Type of equilibrium
β
<
β
R
0
A
>
0, B>0, C>0, D>0
Disease-free equilibrium
β
R
0
<
β
<
β
C
A
>
0, B>0, C>0, D<0
Unique endemic equilibrium
β
C
<
β
<
β
B
A
>
0, B>0, C<0, D<0
Unique endemic equilibrium
β
>
β
B
A
>
0, B<0, C<0, D<0
Unique endemic equilibriumTable 3
Qualitative behaviour for system (1) as function of the disease transmission rate β, when the condition βB<βC<βR0 is fulfilled. Here, Δ1 is the discriminant of the cubic polynomial (20).
Interval
Coefficients
Type of equilibrium
β
<
β
B
A
>
0, B>0, C>0, D>0
Disease-free equilibrium
β
B
<
β
<
β
C
A
>
0, B<0, C>0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
C
<
β
<
β
R
0
A
>
0, B<0, C<0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
R
0
<
β
A
>
0, B<0, C<0, D<0
Unique endemic equilibriumFigure 2
Bifurcation diagram (solutionx of polynomial (20) versus β) for the condition βR0<βC<βB. βR0 is the bifurcation value. The blue branch in the graph is a stable endemic equilibrium which appears for R0>1.From Table2 we can see that although the signs of the polynomial coefficients may change, other new biologically meaningful solutions (nonnegative solutions) do not arise in this case. The system can only display the presence of two equilibrium states: disease-free or a unique endemic equilibrium.The basic reproduction numberR0 in this case explains well the appearance of the transcritical bifurcation, that is, when a unique endemic state arises and the disease-free equilibrium becomes unstable (see blue line in Figure 2).However, the change in signs of the polynomial coefficients modifies the qualitative type of the equilibria. This fact is shown in Figures5 and 7 illustrating the existence of focus or node type steady-sate solutions. These different types of equilibria as we will see in the next section cannot be explained using solely the reproduction number R0.In the next section we will explore numerically the parametric space of system (1), looking for different qualitative dynamics of TB epidemics. We will discuss in more detail how dynamics depends on the parameters given in Table 1, especially on the transmission rate β, which will be used as bifurcation parameter for the model.Let us consider here briefly two examples of parametric regimes for the model in order to illustrate the possibility to encounter a more complex dynamics, which cannot be solely explained by changes in the value of the basic reproduction numberR0.Example I. Suppose β=βR0, this implies that R0=1 and D=0; therefore, we have the equation:
(22)P(x)=Ax3+Bx2+Cx=x(Ax2+Bx2+C)=0.It is easy to see that besides zero solution, ifB<0, C>0 and B2-4AC>0,  (22) has two positive solutions x1 and x2. So, we have in this case three nonnegative equilibria for the system.The conditionB<0 for β=βR0 means fB(βR0)<0, and this in turn implies that βB<βR0. On the other hand, the condition C>0 implies fC(βR0)>0 and therefore βR0<βC.Gathering both inequalities we can conclude that ifβB<βR0<βC, then the system has the possibility of multiple equilibria.Since the coefficientsA and B are both continuous functions of β, we can always find a ϵ neighbourhood of βR0, |β-βR0|<ϵ such that the signs of these coefficients are preserved. Although in this case we do not have the solution x=0, we eventually could still have two positive solutions and consequently, multiple equilibrium states; see the green line in Figure 3.Figure 3
PolynomialP(x) for different values of β with the condition βB<βR0<βC. The graphs were obtained for values of δ=3.0 and η=2.2. The dashed black line indicates the case β=βR0. The figure shows the existence of multiple equilibria.Example II. Suppose we take numerical values for the parameters in Table 1 such that the condition βB<βC<βR0 is fulfilled.Ifβ<βB, then all coefficients of the polynomial (20) are positive and there is not nonnegative solutions. In this case, the system has only a disease-free equilibrium.ForβB<β<βC and βC<β<βR0 the signs of the coefficients of the polynomial are A>0, B<0, C>0, and D>0,A>0, B<0, C<0, D>0, respectively.In both cases the polynomial has two possibilities:(a)
three real solutions: one negative and two positive solutions forΔ1<0,(b)
one negative and two complex conjugate solutions forΔ1>0.HereΔ1 is the discriminant for the polynomial (20).In the (a) case we have the possibility of multiple endemic states for system (1). This case is illustrated in numerical simulations in the next section by Figures 8 and 9.We should note that the valueβ=βB is not a bifurcation value for the parameter β.Ifβ=βB, then A>0, B=0, C>0, and D>0. In this case we have
(23)Δ1=14D2A2+127C3A3>0.The discriminantΔ1 is a continuous function of β, for this reason this sign will be preserved in a ϵ neighbourhood of βB.We should be able to find a bifurcation value solving numerically the equation(24)Δ1(β*)=0,
where β* can be bounded by the interval βB<β*<βR0 (see Figure 4).Figure 4
Bifurcation diagram for the conditionβB<βC<βR0. β* is the bifurcation value. The blue branch in the graph is a stable endemic equilibrium which appears even for R0 < 1.Figure 5
Numerical simulation forR0=3.585422172, δ=0.9, η=0.01, and β=0.00052. The system goes toward a focus type stable stationary equilibrium.


## 4. Numerical Simulations

In this section we will show some numerical simulations with the compartmental model (1). This model has fourteen parameters that have been gathered in Table 1. In order to make the numerical exploration of the model more manageable, we will adopt the following strategy.(i)
First, instead of fourteen parameters we will reduce the parametric space using four independent parametersβR0, β, δ, and η. The parameters β, δ, and η are the transmission rate of primary infection, exogenous reinfection rate of latently infected, and exogenous reinfection rate of recovered individuals, respectively. βR0 is the value of β such that basic reproduction number R0 is equal to one (or the value of β such that coefficient D in the polynomial (20) becomes zero). On the other hand, βR0 depends on parameters given in the list Λ={Π,c,f,μ,ν,p,q,w,μt,r1,r2}. This means that if we keep all the parameters fixed in the list Λ, then βR0 is also fixed. In simulations we will use βR0 instead of using basic reproduction number R0.(ii)
Second, we will fix parameters in the listΛ according to the values reported in the literature. In Table 4 are shown numerical values that will be used in some of the simulations, besides the corresponding references from where these values were taken. Mostly, these numerical values are related to data obtained from the population at large, and in the next simulations we will change some of them for considering the conditions of extremely high incidence/prevalence of TB in semiclosed communities. In any case, these changes will be clearly indicated in the text.(iii)
Third, for any pairs of valuesδ and η we can compute βB and βC, that is, the values of β such that B=0 and C=0, respectively, in the polynomial (20). So, we have that the exploration of parametric space is reduced at this point to the study of the parameters βR0, βB, βC, and β. According to the chosen values for δ, η, and βR0, we have six possible orderings for the parameters βR0, βB, and βC (see Appendix B).The dynamic behavior of system (1) will depend of these orderings. In particular, from Table 5, it is easy to see that if β≤min(βR0,βB,βC) then the system has a unique equilibrium point, which represents a disease-free state, and if β≥max(βR0,βB,βC), then the system has a unique endemic equilibrium, besides an unstable disease-free equilibrium.(iv)
Fourth and finally, we will change the value ofβ, which is considered a bifurcation parameter for system (1), taking into account the previous mentioned ordering to find different qualitative dynamics.Table 4
Numerical values for the parameters in the listΛ. Some of the given numerical values for the model parameters are mainly related to the spread of TB in the population at large and are basically taken as reference. Other values assuming for the parameters, different than those given in this table will be clearly indicated in the text.
Parameter
Description
Value
Π
Recruitment rate
200 (assumed)
c
Natural cure rate
0.058 [23, 33, 34]
ν
Progression rate from latent TB to active TB
0.0256 [33, 34]
μ
Natural mortality rate
0.0222 [2]
μ
T
Mortality rate due to TB
0.139 [2, 33]
w
Relapse rate
0.005 [2, 33, 34]
q
Probability to develop TB (slow case)
0.85 [2, 33]
f
Probability to develop TB (fast case)
0.70 [2, 33]
p
Proportion of new infections that produce active TB
0.05 [2, 33, 34]
r
1
Treatment rates forII
0.50 (assumed)
r
2
Treatment rates forIN
0.20 (assumed)Table 5
Different possible orderings forβR0, βB, and βC. In every case A>0, Δ1 is the cubic discriminant of the equation P(x)=0, Δ is the discriminant of the quadratic equation P′(x)=0, where P(x) is the polynomial (20).
Interval
Coefficients
Equilibria
β
R
0
<
β
C
<
β
B
β
<
β
R
0
B
>
0, C>0, D>0
Disease-free equilibrium
β
R
0
<
β
<
β
C
B
>
0, C>0, D<0
Unique endemic equilibrium
β
C
<
β
<
β
B
B
>
0, C<0, D<0
Unique endemic equilibrium
β
B
<
β
B
<
0, C<0, D<0
Unique endemic equilibrium
β
C
<
β
R
0
<
β
B
β
<
β
C
B
>
0, C>0, D>0
Disease-free equilibrium
β
C
<
β
<
β
R
0
B
>
0, C<0, D>0
Two equilibria ifP(x2)≤0 or Δ1<0; none if P(x2)≥0 or Δ1>0
β
R
0
<
β
<
β
B
B
>
0, C<0, D<0
One equilibrium forΔ1<0 or Δ1>0
β
B
<
β
B
<
0, C<0, D<0
Unique endemic equilibrium
β
C
<
β
B
<
β
R
0
β
<
β
C
B
>
0, C>0, D>0
Disease-free equilibrium
β
C
<
β
<
β
B
B
>
0, C<0, D>0
Two equilibria ifP(x2)≤0 or Δ1<0; none if P(x2)≥0 or Δ1>0
β
B
<
β
<
β
R
0
B
<
0, C<0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
R
0
<
β
B
<
0, C<0, D<0
Unique endemic equilibrium
β
R
0
<
β
B
<
β
C
β
<
β
R
0
B
>
0, C>0, D>0
Disease-free equilibrium
β
R
0
<
β
<
β
B
B
>
0, C>0, D<0
Unique endemic equilibrium
β
B
<
β
<
β
C
B
<
0, C>0, D<0
One equilibrium (Δ1>0), three equilibria (Δ1<0)
β
C
<
β
B
<
0, C<0, D<0
Unique endemic equilibrium
β
B
<
β
R
0
<
β
C
β
<
β
B
B
>
0, C>0, D>0
Disease-free equilibrium
β
B
<
β
<
β
R
0
B
<
0, C>0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
R
0
<
β
<
β
C
B
<
0, C>0, D<0
One equilibrium (Δ1>0), three equilibria (Δ1<0)
β
C
<
β
B
<
0, C<0, D<0
Unique endemic equilibrium
β
B
<
β
C
<
β
R
0
β
<
β
B
B
>
0, C>0, D>0
Disease-free equilibrium
β
B
<
β
<
β
C
B
<
0, C>0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
C
<
β
<
β
R
0
B
<
0, C<0, D>0
Two equilibria (Δ1<0) or none (Δ1>0)
β
R
0
<
β
B
<
0, C<0, D<0
Unique endemic equilibriumIt is especially interesting to explore the consequences of modifications in the values of the reinfection parameters without changing the values in the listΛ, because in this case the threshold βR0 remains unchanged. Thus, we can study in a better way the influence of the reinfection in the dynamics of the TB spread.The values given for the reinfection parametersδ and η in the next simulations could be extreme, trying to capture this way the special conditions of high burden semiclosed communities.Example I (Case
β
R
0
<
β
C
<
β
B
,  δ=0.9,η=0.01). Let us consider here the case when the condition βR0<βC<βB is met. We know from the previous section that this condition is met under biologically plausible values (δ,η)∈[0,1]×[0,1].According to Lemmas3 and 4, in this case the behaviour of the system is characterized by the evolution towards disease-free equilibrium if β<βR0 and the existence of a unique endemic equilibrium for β>βR0. Changes in the parameters of the list Λ alter the numerical value of the threshold βR0 but do not change this behaviour.First, we consider the following numerical values for these parameters:δ=0.9, η=0.01, and β=0.00052. We also fix the list of parameters Λ according to the numerical values given in Table 4.The basic reproduction number for these numerical values givesR0=3.585422172. The initial conditions considered were
(25)S(0)=4980,E(0)=0,II(0)=20,IN(0)=0,R(0)=0.We also have the following values:(26)βR0=0.0001450317354,βB=0.01087387065,βC=0.0002715343808.These values clearly meet the conditionβR0<βC<βB, and according to Lemma 3 the system must have in this case a unique endemic equilibrium for all β>βR0.Figure5 shows that under the above described situation, the system will converge to an endemic equilibrium given by the focus type stationary steady solution:
(27)S∞=1616,R∞=4080,IN∞=103,II∞=195,E∞=1150.By straightforward calculations we can show that this focus is stable, and no matter what initial conditions are taken for the system, the solutions always evolve to this endemic state.Figure6 shows the trajectories of the system for multiple initial conditions in a three-dimensional phase space in which the horizontal axes are susceptible S and recovered R individuals, while the vertical axis is the prevalence II+IN+E.Figure 6
Phase space representation of the evolution of the system toward a stable focus type equilibrium. In this representation were used multiple initial conditions and the following values:R0=3.585422172, δ=0.9, η=0.01, and β=0.00052.Figure 7
Numerical simulation forR0=3.585422172, δ=0.01, η=0.9, and β=0.00052. In this case the system converges to a stable node type equilibrium.Figure 8
Numerical simulation forR0=0.9653059690, δ=3.0, and η=2.5. The system can evolve to two different equilibria II∞=0 (red lines) or II∞=285 (dark green lines) according to different initial conditions.Figure 9
Numerical simulation forR0=0.9653059690, δ=3.0, and η=2.5. Phase space representation of the system with multiple equilibrium points.Example II (Case
β
R
0
<
β
C
<
β
B
,  δ=0.0,η=0.9).  For our next numerical simulation we consider the following values for the used parameters: δ=0.01, η=0.9, β=0.00052, and as before the list of parameters Λ is fixed according to Table 4.The basic reproduction number for these parameters as before gives the same valueR0=3.585422172. The used initial conditions were
(28)S(0)=4980,E(0)=0,II(0)=20,IN(0)=0,R(0)=0.We also have the following values:(29)βR0=0.0001450317354,βB=0.01226355348,βC=0.0003132229272.These values meet the conditionβR0<βB<βC, and as in the previous simulation the system evolves toward a unique endemic equilibrium, but this time the dynamical properties of the equilibrium have changed.In fact, Figure7 shows the evolution of the system toward a stable node type endemic equilibrium:
(30)S∞=1938,R∞=974,IN∞=60,II∞=156,E∞=4530.In our model, considering biologically plausible domain for exogenous reinfection parameters(δ,η)∈[0,1]×[0,1], the condition βR0<βC<βB is fulfilled. Under this condition we have a unique endemic equilibrium for β>βR0. The emergence by a transcritical bifurcation of this endemic state is properly explained by the basic reproduction number R0. However, changes in the reinfection parameters δ, η can modify the qualitative nature of the dynamics of the disease, in addition to changing the numbers of individuals in the different compartments of the model in the endemic equilibrium state, without having any change in the value of the basic reproduction number R0, which in this case fails to describe these variations in the dynamics of the disease.Example III (Case
β
B
<
β
C
<
β
R
0
,
δ
=
3
.
0
,
η
=
2
.
5
).  There is now evidence that rates of secondary tuberculosis in high endemic communities (for example semiclosed communities), in patients with LTB or/and already treated for primary disease, are actually higher than in people presenting with primary infection [21, 22]. Taking this into consideration we consider now the following numerical values for the parameters: β=0.00014, δ=3.0, η=2.5. In this case the basic reproduction number takes the value R0=0.9653059690. Additionally we have
(31)βR0=0.0001450317354,βB=0.0001066568066,βC=0.0001225687204.For these parameter we have that the conditionβB<βC<βR0 is fulfilled and the system has the possibility of multiple equilibria. In fact, we have in this case the following stationary points P=(S,R,Ii,In,E):
(32)P1=(9009,0,0,0,0),P2=(8507,182,9,5,2166),P3=(3221,1406,285,103,1566).P1 is a stable disease-free equilibrium point (stable node), P3 is a stable endemic equilibrium (stable focus), and P2 is an unstable equilibrium point (saddle point).Figure8 shows the convergence to II∞=0 or to II∞=285  according to  with different initial conditions.In Figure9 is shown another representation (phase space) of the evolution of the system toward P1 or to P3 according to different initial conditions. The representation is a three-dimensional phase space in which the horizontal axes are susceptible S and recovered R individuals, while the vertical axis is the prevalence II+IN+E.For the previously numerical values, the system experiences abackward bifurcation [37] at the value β*=0.0001261648723 with βB<β*<βR0. For β>β*, the system possesses two stable equilibrium points and one unstable (see Figure 4).Example IV (Case
β
B
<
β
R
0
<
β
C
,
δ
=
3
.
0
,
η
=
2
.
5
). Consider now a more extreme situation with η=2.5, δ=3.0, and p=0.7 (the other parameters kept the same values given in Table 4). In this case the condition βB<βR0<βC is fulfilled.This example is shown in order to illustrate more complex and rich dynamics that might admit system (1), which is mathematically possible and could in principle be a model case for an extreme hypothetical situation in a semiclosed high burden community. For these parameters we have
(33)βR0=0.0001679568390,βC=0.0001729256777,βB=0.0001489092005,
which clearly satisfy the condition βB<βR0<βC. Therefore, as was explained in the previous section, the system has the possibility of multiple equilibria.In fact, for the bifurcation valueβ1=0.0001673533706 of the disease transmission rate, which satisfies the condition βB<β1<βR0, the system acquires two positive equilibria, apart from the disease-free equilibrium.Whenβ=βR0 appear three positive equilibrium points and the disease-free equillibrium becomes unstable. For β2=0.0001688612368 with βR0<β2<βC the system admits a unique and stable endemic equilibrium (see Figure 10).Figure 10
Bifurcation diagram (solutionx of polynomial (20) versus β) for the condition βB<βR0<βC. The system experiences multiple bifurcations at β1, βR0, and β2.We take now the valueβ=0.0001675, which satisfies the condition β1<β<βR0.With these numerical values the basic reproduction number isR0=0.9972800211<1, and therefore, the disease-free equilibrium is stable.We have in this case the following stationary pointsP=(S,R,Ii,In,E):
(34)P0=(5148,0,0,0,0),P1=(3372,1041,122,60,482),P2=(2828,1283,190,88,651).P0 is the stable disease-free equillibrium point (stable node), P1 is an unstable equilibrium point (saddle point), and P2 is a stable endemic equilibrium (stable focus). Figure 11 shows the convergence to II∞=0 or to II∞=190 according to the initial condition.Figure 11
Numerical simulation forR0=0.9972800211, δ=3.0, and η=2.5. The system can evolve to two different equilibria II∞=0 or II∞=190 according to the initial condition.In Figure12 is shown another representation (phase space) of the evolution of the system toward P0 or to P2 according to the initial conditions.Figure 12
Numerical simulation forR0=0.9972800211, δ=3.0, and η=2.5. Phase space representation of the system with multiple equilibrium points.Let us take now the valueβ=0.0001683, which satisfies the condition βR0<β<β2. In this case, the basic reproduction number has the value R0=1.002043150. We still have that the condition βB<βR0<βC is fulfilled and the system in this case has four equilibrium points P=(S,R,Ii,In,E):
(35)P0=(5148,0,0,0,0),P1=(5042,76,5,3,20),P2=(3971,734,69,36,298),P3=(2491,1413,246,109,750).P0 is the unstable disease-free equillibrium point (saddle point ), P1 is a stable endemic equilibrium point (node), P2 is an unstable equilibrium (saddle point), and P3 is a stable endemic equilibrium point (focus).Figure13 shows the phase space representation of this case.Figure 13
Numerical simulation forR0=1.002043150, δ=3.0, and η=2.5. The system can evolve to two different equilibria P1 (stable node) or P3 (stable focus) according to the initial condition. P0 and P2 are unstable equilibria.For further numerical analysis, we set all the parameters in the listΛ according to the numerical values given in Table 4, leaving free the parameters β, η, and δ related to the primary transmission rate and reinfection rates of the disease.We will explore the parametric space of system (1) and relate it to the signs of the coefficients of the polynomial (20).In Figure14, we consider values of β such that R0>1. We can observe from this figure that as the primary transmission rate of the disease β increases, and with it the basic reproduction number R0, the system under biological plausible condition, represented in the figure by the square (δ,η)∈[0,1]×[0,1], evolves such that initially (for lower values of β) coefficients B and C are both positive, then B remains positive and C becomes negative and finally both coefficients become negative.Signs of coefficientsB and C as functions of exogenous reinfection rate of latent δ and exogenous reinfection rate of recovered η for R0⩾1. The parameter β has the values: (a) β=βR0=0.0001450317354, (b) β=0.0002450317354, (c) β=0.0003450317354, and (d) β=0.001145031735.
(a)
R
0
=
1
(b)
R
0
=
1
.
689504264
(c)
R
0
=
2
.
379008527
(d)
R
0
=
7
.
895042636This change in the coefficients signs as the transmission rateβ increases agrees with the results summarized in Table 2 when the condition βR0<βC<βB is fulfilled.Next, in order to explore another mathematical possibilities we will modify some numerical values for the parameters in the listΛ in a more extreme manner, taking a hypothetical regime with Λ*={μ=0.03885, μt=0.01520, p=0.8, ν=0.0266, f=0.8, q=0.85, w=0.005, c=0.4, r1=0.5, r2=0.2}.In Figure15 besides signs of B and C we consider also the signs of the discriminant Δ of the quadratic equation P′(x)=0, where P(x) is the polynomial (20).Signs of coefficientsB, C, and discriminant Δ=B2-3AC as functions of exogenous reinfection rate of latent δ and exogenous reinfection rate of recovered η for R0⩾1. The parameter β has the values: (a) β=βR0=0.0002277727471, (b) β=0.0002287727471, (c) β=0.0002477727471, and (d) β=0.0005277727471.
(a)
R
0
=
1
(b)
R
0
=
1
.
004390340
(c)
R
0
=
1
.
087806817
(d)
R
0
=
2
.
317102260From Figure15, in particular we can see that the domain with B<0, C>0, Δ>0, and P(x1)>0 represented in red allows the possibility of multiple endemic equilibria for system (1). However, despite the extreme numerical values for the parameters taking in Λ*, this domain is still far from the domain of biologically plausible values for δ and η, represented in the figures by the square (δ,η)∈[0,1]×[0,1].As the transmission rate of the diseaseβ increases, and with it the basic reproduction number R0, this red domain of the parametric space becomes increasingly smaller until finally it disappears. In fact, the red domain is only significant when the basic reproduction number R0 is near one.In Figure16, we show some numerical simulation where basic reproduction number R0 is less than one. The red domain indicates the possibility of multiple endemic equilibria for the system even for R0<1. We can see that this domain is far from the domain of biologically plausible values for δ and η represented in the figure by the square (δ,η)∈[0,1]×[0,1]. As the transmission rate of the disease β decreases, and with it the number R0, this parameter domain moves away from the square. In all represented cases B>0 and C>0 inside the square (δ,η)∈[0,1]×[0,1].Signs of coefficientsB, C, and discriminant Δ=B2-3AC as functions of exogenous reinfection rate of latent δ and exogenous reinfection rate of recovered η for R0⪕1. The parameter β has the values: (a) β=0.0002177727471, (b) β=0.0002027727471, (c) β=0.0001777727471, and (d) β=0.0001277727471.
(a)
R
0
=
0
.
9560965911
(b)
R
0
=
0
.
8902414781
(c)
R
0
=
0
.
7804829561
(d)
R
0
=
0
.
5609659131


## 5. Discussion and Conclusions

In order to consider high incidence and prevalence of TB and LTB in semiclosed communities, we have used in this work a compartmental SEIR model with five possible pathways to TB disease. The extra nonlinear terms considered in the model lead to a more complex mathematical treatment in comparison with previously used models (see e.g., [23, 25–29, 32, 35, 36]). But the special form of some coefficients obtained from the analysis of standard SEIR models with constant transmission rate allowed us to move forward with some analytical results that were confirmed later by numerical simulations.In this paper we follow a new methodology for numerical exploration of this kind of models, which allowed us to handle the high dimensionality of its parameter spaces and thus study its different dynamic behaviors. We found that the transmission rateβ can be used as a bifurcation parameter and that the system undergoes qualitative changes in the dynamics when this parameter varies. In this context the analysis of the parametric space is reduced to the study of variations of the parameters βR0, βB, βC, and β. We divided the parametric space into six possible arrangements for these parameters, which in turn determine all possible different qualitative behaviours of the system dynamics.From model (1) we can see that reinfection requires latently infected, recovered, and actively infectious individuals. The basic reproduction number R0 has to do solely with infections produced when an infectious individual is introduced into an uninfected population, but reinfection does not alter the dynamics of an uninfected population. So, the number R0 does not completely describe the dynamics of the model when the reinfection is incorporated, as was noted before by Feng et al. in [26]. Unlike the model published by these authors, which uses a single parameter for exogenous reinfection, in our model we use two parameters related to two possible pathways of reinfection (reinfection of latently infected and reinfection of recovered individuals).This is a reason why our model shows a more complex and richer dynamics.We have showed through theoretical analysis in Section3 and numerical simulations in Section 4 that if we accept as valid the plausible assumption that exposure to mycobacterium induces an immune response, which is partially protective against reinfection, then the system for semiclosed communities (1) reproduces well, common observed trends in TB epidemiology that are similar to what happens in population at large, which is basically that, for R0<1, there is only one disease-free status, while for R0>1, there exists a unique endemic state with nonzero prevalence. For R0=1 occurs a transcritical bifurcation from which emerges an endemic stable state.Moreover, according to Lemmas3 and 4, any values of reinfection parameters in this parametric regime: (δ,η)∈[0,1]×[0,1] would lead to the same qualitative dynamics and will not affect this already classical behavior in SEIR models. In this case only one of the aforementioned arrangements (βR0<βC<βB) emerges as valid under this biologically plausible condition.Since the two parameters related to exogenous reinfection of latently infected and recovered individuals do not affect the value of the numberR0, even under the plausible assumption of partial immunity, variation of reinfection parameters can make that for the same value of the number R0, the quality of dynamics and the number of affected by disease individuals (incidence and prevalence) drastically change. For example, Figures 5 and 7 show two types of dynamics, that is, convergences to different stationary points, a focus and a node for the same basic reproduction number R0. Some evidence of this variability in tuberculosis epidemiology due to dynamic balance between primary infection and reinfection has been presented in several works (see e.g., [26, 38]).Taking less plausible assumption, but already evidenced in several works [5, 21, 22, 26], of an increased susceptibility to reinfection over primary infection in some cases leads us to a further study of model (1). For δ>1 and η>1 system (1) experiences a rich and complex dynamics with successive and different kind of bifurcations as the transmission rate β changes. These cases incorporate possible multiple endemic states, regardless of whether the values for the basic reproduction number R0 were less than or greater than 1. So, these behaviors cannot be explained using only this number. It is in this context that the use of the disease transmission rate β as bifurcation parameter instead of R0 acquires real usefulness.Some important implications of the simulations withδ>1 and η>1 lie in the fact that many of the measures taken to stop and control an epidemics are designed to reduce the value of the basic reproduction number R0 such that disease-free status for R0<1   is achieved. However, in this parametric regime, reinfection might cause the system to fall into a state unable to eliminate endemic disease, although it fulfills that R0<1. Thus, semiclosed communities with this kind of regime will become in genuine high transmission pockets of TB inserted in the general population [4]. Indeed, semiclosed communities such as prisons might become in a reservoir for disease transmission to the population at large and should be a source of public concern [4, 6, 7].The theoretical approach and numerical simulations presented in this paper for the study of the impact of reinfection on TB dynamics in semiclosed communities could have important implications at multiple levels, including vaccine design, control program design, epidemiology of tuberculosis in regions where the risk of reexposure is high, and for systems-based computer models which to date assume that primary infection will confer at least some degree of (stable) memory immunity to a secondary infection, but that in fact also have to consider less plausible assumptions about an increased susceptibility to reinfection.



---

*Source: 648291-2013-05-09.xml*