Backward bifurcation and hysteresis in models of recurrent tuberculosis

An epidemiological model is presented that provides a comprehensive description of the transmission pathways involved for recurrent tuberculosis (TB), whereby cured individuals can become reinfected. Our main goal is to determine conditions that lead to the appearance of a backward bifurcation. This occurs when an asymptotically stable infection free equilibrium concurrently exists with a stable non-trivial equilibria even though the basic reproduction number R0 is less than unity. Although, some 10-30% cases of TB are recurrent, the role of recurrent TB as far as the formation of backward bifurcation is concerned, is rarely if ever studied. The model used here incorporates progressive primary infection, exogenous reinfection, endogenous reactivation and recurrent TB as transmission mechanisms that contribute to TB progression. Unlike other studies of TB dynamics that make use of frequency dependent transmission rates, our analysis provides exact backward bifurcation threshold conditions without resorting to commonly applied approximations and simplifying assumptions. Exploration of the model through analytical and numerical analysis reveal that recurrent TB is sometimes capable of triggering hysteresis effects which allow TB to persist when R0 < 1 even though there is no backward bifurcation. Furthermore, recurrent TB can independently induce backward bifurcation phenomena if it exceeds a certain threshold.


Introduction
In 1882 Mycobacterium tuberculosis was identified by Robert Koch [1] as the aetiological agent responsible for tuberculosis (TB). Despite all attempts to control its spread by modern medical science, TB has become one of the most widespread and serious of all infectious diseases today [2]. The disease is transmitted from one person to another in tiny microscopic droplets when a person with pulmonary TB expels bacteria into the air by either through coughing, sneezing, singing, laughing or other related activities that involve airborne pathways. Amongst all infectious diseases, TB is one of the leading causes of death worldwide, and second only to human immunodeficiency virus (HIV) [3,4]. Approximately a quarter of the global population harbours the TB bacteria and another eight to nine million new cases of tuberculosis emerge every year [5]. Extraordinarily, TB is a treatable disease and can be prevented and cured through the use of prophylaxis and therapeutics for individuals with latent and clinically active TB respectively. Such treatment should in theory be an effective strategy for controlling the spread of TB. However, it has failed in practice due to the inability to distribute sufficient drug treatment, usually in the form of antibiotics, to the world's population combined with the difficulties of ensuring compliance to the required lengthy treatment program. Moreover, erratic treatment has led to the evolution of multi-drug resistant tuberculosis, giving rise to the fear that TB may become an untreatable disease in the not too distant future. These problems, taken together, have led the World Health Organization (WHO) to formulate a post-2015 global "Stop TB Strategy" [2] to "end the global TB epidemic." A TB episode may have an exogenous or endogenous origin. Exogenous refers to a disease episode that results from recent exposure to some external infectious source (typically, contact with an infectious person). Endogenous designates situations where the individual is already harbouring the causing agent, which is under some healthy control by the immune system, but which destabilizes and leads to disease. TB is unusual in that it can arise via both exogenous and endogenous infection. The pathogenesis of TB is characterized by the infection either remaining dormant, often for a long period that may last years, or progress directly to active TB where clinical symptoms immediately manifest. The latter process is referred to as fast primary progression. The particular course of the disease depends on the host's immune response towards the tubercle bacilli. Thus, the exposure to tubercle bacilli does not necessarily result in the manifestation of clinical forms of TB. Studies suggest that only 5-10% of individuals progress directly to the active stage after exposure to bacilli [6,7]. The other component of the population of exposed individuals develop dormant TB and may remain latently infected, possibly for the rest of their lifetime. However, destabilization of the immune system by the pathogen within the latently infected host can trigger endogenous reactivation, in which latent bacilli are reactivated and cause clinical Mycobacterium tuberculosis. The lifetime risk of a latently infected individual to progress to the infectious stage is approximately 5-10% [8].
We are particularly interested in recurrent TB, which is defined as the emergence of a second episode of TB after the first episode has been successfully cured [9]. This often arises following treatment, because an individual that recovers from a first episode of the disease does not necessarily gain permanent immunity to a second. Approximately 10-30% of all cases of TB are due to recurrent tuberculosis [10,11], and multiple episodes are largely attributable to ineffective or poorly implemented tuberculosis control programs. Although recurrent TB is recognized as a serious problem it receives little attention [10]. Through the use of advanced molecular fingerprinting techniques, TB recurrence has been classified into two fundamental forms of infection: i) relapse of the original infecting strain, and ii) reinfection with a new strain of Mycobacterium tuberculosis. The role of reinfection and relapse to the overall burden of tuberculosis recurrence is not well understood and this has potential public health implications [12]. It is important to note that the system adopted by the World Health Organization in recording and reporting TB cases does not differentiate between true relapse (i.e., reactivation of latent TB) and reinfection (exogenous acquisition of TB) and classifies as relapse any recurrence of TB [11]. This greatly affects the collection and analysis of data making it even more difficult to assess the specific role of reinfection.
In [13] it was found that persons who had TB have a greatly increased risk of developing TB when reinfected. Moreover, the study suggests that the incidence rate of TB attributable to reinfection after successful treatment is four times higher than that of new TB [13]. A key goal of this paper is to investigate the effect of recurrent TB on the formation of backward bifurcations through detailed mathematical modelling. To our knowledge, no other study has examined this possibility. Similar to other infectious diseases, the ability for TB to persist in a population or go extinct is to a large degree governed by the basic reproduction number, R 0 , which is defined as the number of secondary cases an infected individual will generate when introduced into a completely susceptible population [14][15][16][17]. This fundamental epidemiological quantity determines whether an infection will be able, on average, to reproduce itself (R 0 > 1) in the population or not (R 0 < 1). Characteristically, when R 0 is below unity, the introduction of a few infected individuals in a susceptible population will only lead to disease die-out, as the disease is unable to reproduce itself or transmit through the population effectively. Conversely, when R 0 is above unity an epidemic may trigger and long-term disease persistence is feasible. Classical epidemic models therefore often have two intrinsic equilibria: a disease free equilibrium (DFE) and a non-trivial endemic equilibrium. By endemic is meant an equilibrium in which the number of infectives is greater than zero. The stability of these equilibria switch at the (transcritical) bifurcation point when R 0 = 1. Fig 1 illustrates the more typical forward bifurcation by plotting the equilibrium number of infectives (I Ã ) in a population as a function of R 0 . This shows that a stable DFE in which I Ã = 0 exists when R 0 < 1. However, when crossing the threshold to a regime where R 0 > 1, there is a change of stability where the DFE becomes unstable while the previously unstable endemic equilibrium stabilizes.

Backward bifurcations
For decades, it has been widely accepted that the condition R 0 < 1 is an essential requirement for the elimination of a disease. However, this viewpoint has been recently challenged. Instead, Recurrent TB induces a rare hysteresis effect the phenomenon of backward bifurcation offers a different interpretation since it shows that although R 0 < 1 and the DFE is stable, there might still be another stable endemic equilibrium coexisting simultaneously. Thus even though R 0 < 1, a population may still reside at an endemic equilibrium at which the disease persists indefinitely. When there are multiple stable equilibria coexisting simultaneously, the final equilibrium a population will reach is dependent on the initial conditions of its sub-populations. Fig 2 provides a typical bifurcation diagram that shows the key features of a backward bifurcation. Note that three equilibria coexist when R 0 is in the range 0 < R c < R 0 < 1, where R c is a critical value, which in Fig 2 is R c = 0.53. In this range, the "middle" equilibrium is unstable, while the other two outer equilibria (the DFE and the endemic equilibrium) are both stable. When R 0 < R c , only the DFE equilibrium exists and is stable.
Multiple coexisting equilibria can lead to interesting dynamics when, for example R 0 varies slowly as shown in Fig 2. Suppose, for example, that the infective population is close to extinction I Ã = 0 (the DFE) and R 0 increases slowly. As soon as R 0 passes through the threshold point R 0 = 1 the number of infectives will suddenly jump from close to I Ã = 0 to the large endemic equilibrium (I Ã > 0), (see Fig 2). Similarly, when the infective population sits close to the endemic equilibrium and R 0 reduces slowly through the threshold point R 0 = 1, rather than switching to the DFE, the infectives remain close and converge to the stable endemic equilibrium. However, as soon as R 0 reduces below the threshold R c , the infective population immediately jumps towards zero (the DFE), while the endemic equilibrium disappears.
Backward bifurcations have major implications for infectious diseases such as TB, since control programs based on reducing R 0 below unity may be ineffective given the disease might be able to easily persist indefinitely under such conditions. Certainly this non-intuitive possibility will have to be taken into account if the new WHO "Stop TB Strategy" is to succeed. Lipsitch and Murray [18] considered the phenomenon of backward bifurcation caused by exogenous reinfection in TB epidemic models with skepticism and suggested that it is unlikely to occur for realistic parameter values. However, TB epidemiology in the last two decades has greatly changed and challenges Lipsitch and Murray [18] line of argument. More recent research supports the notion that individuals experience multiple infections throughout the course of their life as a result of reinfection, especially in high-incidence TB areas [19][20][21], all conditions which could promote the existence of backward bifurcations. More on this controversy will be elaborated in the text.

Model of recurrent tuberculosis
Although, there are numerous TB models that have attempted to include the recurrent TB pathway (for instance see [22][23][24][25][26][27]), none have explored its role in the formation of backward bifurcation. Yang et al. [24] investigated the impact of multiple infections and long latency on the dynamics of recurrent tuberculosis. Their results suggest that a backward bifurcation is expected to occur when a critical value of the disease incubation period is exceeded. However, Yang et al. [24] assume that the reinfection of recovered individuals is negligible, arguing that such a pathway increases non-linearity and makes the model mathematically intractable. The models of Feng et al. [23] and Kar et al. [22] were all based on the assumption that individuals pass through a long latency period before TB reactivated to clinically active TB. However, based on the natural history of TB, fast primary progression of TB forms an important process through which symptomatic TB emerges [28] and needs to be included. Indeed, Feng et al. [23] incorporated a particular recurrent TB pathway in their model; however, their simplifying assumptions hindered further exploration with regard to its role in causing a backward bifurcation. The models developed by Gomes et al. [25] and Herrera et al. [26] attempted to incorporate exogenous re-infection, partial immunity to reinfection and primary progression. However, neither group examined how recurrent TB reinfection pathways could lead to a backward bifurcation. Instead their main objective was to study the reinfection threshold. Hence, our main aim is to study reinfection amongst recovered individuals and deduce the epidemiological implication it has especially with regards to the possible formation of a backward bifurcation. For this purpose we develop a TB model that includes fast primary progression and possible reinfection pathways.
The model is represented graphically in Fig 3 and is based on the following processes. We assume that every individual in the population belongs to one of four broad classes: susceptible individuals (S), exposed individuals (E) (these are infected individuals who are not able to transmit infection), individuals with active TB (I) (who manifests symptoms and are able to pass on the infection), and recovered/treated individuals (R). Individuals have the potential to move through these four classes, as for example in the loop Susceptible individuals are generated by recruitment through births and immigration at a rate Λ. Upon contact with an infective, a small proportion q of infected susceptible individuals follow the fast primary progression route (i.e they move directly to the infective class) while the rest (1 − q) move to the exposed class where they pass through a long latency period before reactivation and becoming infectious. Infected individuals who recover from the disease then move to the recovered class.
The model includes three different pathways for exogenous reinfection as shown in Fig 3: (i) Path A = pλE, where exposed individuals, who have partial immunity against exogenous reinfection (See [23]), become reinfected and progress to active TB, where recovered individuals become reinfected and progress to the exposed sub-population and (iii) Path C = σθλR, where recovered individuals become reinfected and progress to active TB.
In our model, the parameter p measures the degree of partial protection against TB among latently infected individuals and pλE is the exogenous reinfection incidence rate. σ measures the probability of fast progression to the infectious class after reinfection. Following the studies of [22,23,25,26,29], p = 1 implies that the body does not render protection against exogenous reinfection while 0 < p < 1 implies the body is partially immune against exogenous reinfection (i.e latent infection provides immunity against new infections) [27,30]. Note that p > 1 would imply that an individual with latent TB infection has increased susceptibility to become newly infected, as compared to the susceptibility of the general population [27]. This relates to studies which have found that recovered individuals are more likely to be susceptible to future TB infection than TB-naive individuals [13].
The parameter θ (0 < θ < 1) quantifies the amount of exogenous reinfection among TB recovered individuals via paths B and C while σ represent the probability of fast progression after reinfection amongst recovered individuals. θ < 1 indicates that recovered individuals have acquired some degree of partial protective immunity to TB, while θ > 1 indicates increased susceptibility.
The full model equation is given in terms of the rates of change of each of the sub-populations S, E, I, and R namely: Note that the total population at time t is given by N(t), The model assumes a frequency dependent incidence rate. Here we have followed the convention of working with the so-called 'force of infection' λ defined as where c represents the host-host contact. The parameter β is the probability of a contact being infectious [31] and I/N denotes the likelihood that the encounter is with an individual with active TB (assuming random mixing) [32]. This form of incidence is considered to be more appropriate for infections in human populations [33]. That is individuals become infected when they come into contact with infected individuals regardless of the size of the human population [33,34]. It is important to note that although a large number of previous TB models utilised frequency dependence incidence rates, most assume that the disease-induced death (μ d ) is negligible for the purpose of mathematical simplification. Indeed the analysis becomes mathematically difficult or even intractable without making such assumptions, since otherwise the total population of the model N(t) might never remain constant. Both Feng et al. [23] and Porco et al. [7] resorted to this assumption to ensure that the total population N(t) remained constant, even though this was not true of their more general model formulation. In the analysis here, we make no such assumptions or simplifications and consider the total population to be varying in time.
In more detail, as Eq (1) shows, susceptible individuals move to the latently infected class E upon an effective contact with an infected individuals (I) at a rate λS. The exposed sub-population (E) increases with the infection of susceptible individuals at a rate (1 − q)λS and reinfection (recurrent TB) of recovered individuals at a rate (1 − σ)θλR. It decreases by exogenous re-infection (pλEI), and endogenous reactivation (kE) upon which exposed individuals move to the infectious class. The infected sub-population is generated by fast primary progression of TB susceptibles (qλS), exogenous reinfection amongst exposed sub-population (at rate pλE), exogenous reinfection of recovered individuals (σθλR) and endogenous reactivation (k). The sub-population is decreased by per-capita recovery due to treatment r, and by disease induced death at rate μ d . Finally the recovered sub-population (R) is generated by the recovery of infected individuals and exogenous reinfection (at rate θλR). Note that we model natural death which affects all classes of the population at the same rate via the mortality parameter μ.
Our model assumes all immigrants are susceptible. While this excludes the realistic possibility of immigration of infectives (either latent individuals or individuals with active TB), it helps untangle the conditions that result in backward bifurcation. Previous modelling studies have already demonstrated that the immigration of individuals with latent TB or active TB are pathways that can trigger backward bifurcations [35,36]. Progression to active TB is not uniform, as some infected individuals are more likely to progress to active TB than others. A number of models incorporating long and variable rate of progression have been constructed and analysed (for instance see [23,37,38]). Feng et al. [39] investigated the impact of variability in latency using arbitrary continuous distributions and found that such generalization did not result in qualitative differences in terms of the model dynamics. Before the Feng et al. [23] study, Blower et al. [37] formulated a differential equation model with two latent cohorts: one cohort consisted of those who rapidly develop TB after primary infection, while the second cohort involved individuals who develop the infection slowly through endogenous reactivation. Feng et al. [39] showed that the artificial divisions (as in Blower et al. [37]) play no role in the qualitative dynamics.
We also chose a single latent compartment over two latent compartments (fast and slow TB progression) for the purpose of mathematical tractability. The increased number of reinfection pathways needed with two compartments would add a considerable degree of model complexity, and become very difficult to analyse. Moreover, the original study of [23], where the role of reinfection in inducing backward bifurcation was first identified, consisted of a single latent compartment. The paper by [18] that criticised Feng et al. [23] also had a single latent compartment. Since our goal is to investigate how recurrent TB (reinfection of recovered individuals) can impact the backward bifurcation phenomenon found by Feng et al. [23], to keep the same single latent compartment structure. Table (1) provides a detailed list and description of model parameters as well as typical parameter values used here, as obtained from the literature.

Basic properties
Following the methods in Garba et al. [46], it is not difficult to prove that when all model parameters are nonnegative the state variables S(t), I(t), E(t) and R(t) are all positive for all time t. The equilibria of the model are found by setting the rates of all variables in the left-hand side of Eq (1) to zero. Clearly the equations have an intrinsic disease free equilibrium (DFE) given by ð " S; " E; " I; " The stability of the DFE is controlled by the basic reproduction number R 0 which represents the average number of new infections generated by an infected individual when introduced into an entirely susceptible population. R 0 may be determined using the next generation operator method (see [47]) as shown in S1 Appendix, where it is found that: It is possible to decouple the expression for R 0 to account for slow TB progression and fast primary progression The slow TB component of R 0 can be obtained by observing that the average infectious period is given as 1 m þ r þ m d and the probability of progressing from latent compartment to infective class is given as k m þ k : The average time an individual who starts in the latent compartment is expected to spend in the infectious compartment is, Multiplying this average time by (1 − q)cβ yields the slow TB R 0 . Moreover, multiplying mean An important result is the following Theorem: This general result has been reviewed in [48], and hence we do not prove Theorem 1 here. The theorem implies that it is possible to eradicate the disease from the community when R 0 < 1 if the initial sizes of the sub-populations of model (1) are in the basin of attraction of the disease free equilibrium.
Interestingly, the formula for the basic reproduction number Eq (3) does not include the reinfection parameters p and θ despite the fact that these terms should contribute significantly to the emergence of new cases of TB infection. Hence, this already suggests that R 0 alone is unable to completely quantify some key dynamical features of the TB epidemic, and is in fact the first sign that a backward bifurcation might be involved. It will emerge that the reinfection parameters p and θ play an important role and are responsible for the presence of the backward bifurcation intrinsic to this model.

Backward bifurcation analysis
We begin by identifying the model's endemic equilibrium points being mindful that there may in fact be several such points coexisting simultaneously. As before, to find the endemic equilibria (S Ã , E Ã , I Ã , R Ã ) we set the rate Eq (1) to zero and solve for the equilibrium quantities S Ã , E Ã , I Ã and R Ã in terms of the force of infection λ (See S3 Appendix). Substituting these equilibrium quantities into the force of infection Eq (2) yields where a 3 ¼ yp; a 2 ¼ cypðb 0 À bÞ; a 1 ¼ c½qym þ yk þ mpðb 1 À bÞ; One notes from Eq 4 that the reproductive number , and also that the root λ = 0 corresponds to the DFE, where I Ã = 0. Now the roots of the cubic equation substituted in S Ã , E Ã , I Ã , R Ã yield the endemic equilibrium for any specific set of model parameters (see S3 Appendix).

No recurrent TB (i.e θ = 0): Quadratic P 2 (λ)
We first examine the important particular parameter subset in which σ = θ = 0, that is, in the absence of recovered individuals becoming reinfected. This model has the same basic reproduction number as for the case σ, θ > 0, as σ and θ do not appear in the R 0 expression (3). For σ = θ = 0 the third degree polynomial Eq (9) collapses to the quadratic: where In the above, we see that c 0 < 0 corresponds to R 0 > 1 and vice versa. Thus, in the absence of recurrent TB we deduce c 1 < 0, R 0 < 1 and c 2 1 À 4c 2 c 0 > 0 (see Theorem 1 in S1 Appendix) which indicates conditions for the existence of a backward bifurcation, based on the roots of the quadratic Eq (10).
Note that Case (iii) of Theorem 1 in S1 Appendix stipulates the condition that 4 ¼ c 2 1 À 4c 2 c 0 > 0 which means there are two real positive endemic equilibria as required for a backward bifurcation to appear. In fact 4 = 0 provides the critical point for the backward bifurcation where the two positive endemic equilibria collide and annihilate each other leaving the DFE as the only equilibria.
By setting 4 = 0 we can determine the critical value of the transmission coefficient denoted by β c . For mathematical convenience let Now the discriminant 4(β) may be expressed in terms of β. Let β c be the critical value of β for which the discriminant equals zero i.e., The critical value of basic reproduction number denoted by R c is obtained by replacing parameter β in R 0 with β c which yields where the right-hand term in large brackets is just β c . In fact R c defines a sub-threshold domain of bistable equilibria of the model system (1) in the sense that within the region R c < R 0 < 1 the model Eq (1) has two positive endemic equilibria simultaneously existing with a stable disease free equilibrium. Thus, the backward bifurcation for Eq (1) occurs for values of the basic reproduction number R 0 that lie between R c < R 0 < 1. The associated backward bifurcation for the model without the reinfection pathways A and B (i.e, σ = θ = 0) shown in Fig 4 is obtained by plotting λ as a function of β. Fig 4 shows that model (1) has a disease free equilibrium which corresponds to λ = 0 and two non-trivial endemic equilibria which, according to numerical simulations, one is locally asymptotically stable (LAS) and the other is unstable (saddle). Shortly we will apply center manifold to examine stability and confirm coexistence of these three equilibria (see S2 Appendix). We thus find coexistence of two positive equilibria when R 0 < 1, hence confirming that the model exhibits the phenomenon of backward bifurcation for R c < R 0 < 1.
It can be seen from Eqs (11) and (12), that when exogenous reinfection parameter p increases, R c decreases, and vice-versa. Hence, we can deduce from expression (12) that R c is inversely proportional to the level of exogenous reinfection p. This observation is confirmed by Fig 5 which illustrates the effect of increasing exogenous reinfection p on R c . That is, with low values of exogenous reinfection the critical value R c is high implying that the extent of the backward bifurcation regime becomes smaller as R c becomes closer to unity (R 0 = 1). The threshold implies that TB can be eliminated from the community if the basic reproduction number is maintained below R c (i.e R 0 < R c ). More formally we have the following lemma: Lemma 1 For model Eq (1), when σ = θ = 0,

Recurrent TB: Model with all reinfection pathways (A, B and C) θ > 0: Cubic P 1 (λ)
We now return to the fully general model with all parameters p, σ, θ positive. Recall that the sign of the roots of the cubic polynomial (Eq (4)) P 1 (λ), tell us the signs of the equilibrium populations for the number of infected individuals (via Eq (2)). Observe that the coefficients in the cubic polynomial a 3 , a 2 , a 1 and a 0 (see Eq (5)) are all real numbers. For any non-negative model parameters, a 3 is always positive while a 2 , a 1 and a 0 can be either positive, zero or negative depending on β 0 , β 1 and β R , respectively (see Eqs (6)-(8)). A comprehensive analysis of the roots of the cubic (9) may be carried out using Descartes rule of signs as described in S4 Appendix. A simpler more intuitive approach is to examine the roots at the transcritical bifurcation point R 0 = 1. Such an analysis gives an understanding of the type of bifurcation that is likely to occur in the vicinity of R 0 = 1, and provides conclusions that coincide with the more detailed analysis based on Descarte's rule of signs (S4 Appendix). Conveniently when R 0 = 1 the cubic polynomial (9) reduces to the quadratic equation f ðlÞ ¼ a 3 l 2 þ a 2 l þ a 1 ¼ 0: Keeping in mind Eq (5), a simple study of the roots shows that if either a 2 < 0 and a 1 < 0 (i.e., if β > β 0 and β > β 1 ) or if β 1 < β < β 0 , the quadratic equation has one positive endemic equilibria. As seen in Fig 6, this is the signature of a backward bifurcation. Namely, when R 0 is slightly below unity, the model Eq (1) has two positive endemic equilibria but only one when R 0 ! 1. Furthermore, if a 2 > 0 and a 1 > 0 (i.e β < β 0 and β < β 1 ) then at the point R 0 = 1, the model has no positive endemic equilibria. This characteristic is indicative of forward bifurcation as seen in Fig 7. However, if R 0 is increased slightly above unity then model Eq (1) has one positive endemic equilibrium point. It can be shown that if β 0 < β < β 1 , then the reduced equation has two positive real roots, indicating that the model Eq (1) exhibits hysteresis (see Figs 8 and 9), as will be discussed in detail shortly.
This discussion based on the point R 0 = 1 is summarized in Lemma 2.

Existence of backward bifurcation reinfection threshold at R 0 = 1
The concept of existence of a critical point at R 0 = 1 was fully developed by [49]. For the general model (1), we are interested in determining the critical value of exogenous reinfection p c required to allow the formation of a backward bifurcation at R 0 = 1.  Define where The expression F r (Eq 14) is associated with recurrent TB since it is the only term containing the reinfection parameter θ. In the absence of recurrent TB due to reinfection (i.e θ = 0) the expression F r reduces to zero. We now state the following Lemma 3: To see this, first note that referring to Fig 6, it is apparent that for our model a backward bifurcation requires that at R 0 = 1 there is a single positive endemic equilibrium. Note again that when R 0 = 1 (where a 0 = 0), the equilibria relate to the roots of the quadratic f(λ) = a 3 λ 2 + a 2 λ + a 1 = 0. The roots are: Recurrent TB induces a rare hysteresis effect We are interested in the point where the bifurcation changes from forward to backward. But when R 0 = 1 (β = β R ) this is just the point where the positive endemic equilibrium vanishes to zero, and from Eq (15) must occur when a 1 = 0 (β = β 1 ); or equivalently β 1 = β R (see Eq (5)). Also note that β 1 is a function of p (i.e. β 1 = β 1 (p)) while β R is not. We therefore equate Eqs (7) and (8) for β 1 (p) and β R and solve to find the critical value p c for which β 1 (p c ) = β R . After some algebraic manipulation we find the required backward bifurcation reinfection threshold is given by Eq 13. A more detailed derivation of p c can be found in S2 Appendix, which takes into account the stability of the equilibria through the use of center manifold theory.
Numerically, Lemma 3 is illustrated by Figs 6 and 7 which respectively show that model Eq 1 has a backward bifurcation when p > p c and forward bifurcation when p < p c . Recurrent TB induces a rare hysteresis effect

Relation with models of Lipsitch & Murray [18] and Feng et al. [23]
It is interesting to compare the critical point p c for the reinfection threshold given in Eq (13) with that found by Feng et al. [23] in their much simpler but still important model. First note that Feng et al. [23] do not include fast primary progression which is equivalent to setting q = 0 (see Fig 3). They also assume that the exogenous reinfection rates are θ = 1, σ = 0 and that the disease induced death rate is negligible with μ d = 0. Under these conditions Feng et al. [23] find that backward bifurcations occur only if p > P feng , where and the latter approximation assumes that k ( r. Recurrent TB induces a rare hysteresis effect In their controversial paper, Lipsitch and Murray [18] argued that in the real world recovered individuals gain immunity to reinfection, and thus reinfection amongst exposed individuals must be less than the probability of progressing to the infectious stage of TB. They showed that this implies given that k ( μ. As such, Lipsitch and Murray [18] argued that backward bifurcations should not be expected in the real world. Our extended model has some interesting insights with regard to these studies. First note that after inclusion of the assumptions made by Feng et al. [23] in our model (eg., setting q = 0), the threshold p c for backward bifurcation (Eq (13)) simplifies to That is, backward bifurcations are possible only when p > p c % p feng , and thus the result of Feng et al. [23] is retrieved.
Consider now the extended model with more realistic reinfection pathways (0 < q ( 1, σ > 0), but still assuming that μ d = 0, k ( r and μ ( r. For this approximation F r % σθ, and the backward bifurcation threshold Eq 13 simplifies to: It is interesting to compare this result to the Lipsitch and Murray [18] threshold P L discussed above. If recovered individuals gain high immunity from having been infected, then the reinfection pathway R ! I is relatively small (see Fig 3), with σθ ( 1. In this case Eq (17) shows that p c % P L + q. This is similar to the [18] criterion, and suggests that backward bifurcation is unlikely to occur in the real world, if as Lipsitch and Murray [18] claim that in reality p < P L . However, if recovered individuals gain only mild immunity against reinfection, then the reinfection pathway R ! I and σθ can be relatively large. Note that θ can be greater than unity as suggested by Verver et al. [13]. Gomes et al. [41] estimate θ to be in the range [1.61, 7.79]. In this situation it is quite possible that p c < P L . Hence, with relatively low immunity amongst recovered individuals [13] backward bifurcation can occur despite the fact that the Lipsitch and Murray [18] prediction would literally predict otherwise. This does not mean that Lipsitch and Murray [18] have erred, but that their result may need modifications when discussing the presence of more complex reinfection pathways. In fact even just for the simplified reinfection pathways of the original Feng et al. [23] model, the validity of the Lipsitch and Murray [18] argument has recently been called into question given the difficulties of comparisons with real world processes (Wang et al. [32]).
In the more recent literature, numerous studies have pointed out that initial infection may not confer protection against exogenous reinfection especially in high-risk populations [50,51]. Some studies have demonstrated that it is possible for exogenous reinfection to outweigh endogenous reinfection [20,21]. This is supported by the fact that the majority of new TB cases (about 90%) occur as a result of reinfection rather than endogenous reactivation [19][20][21]43]. In this situation, p c > k/(μ + k), and backward bifurcations can occur even according to the Lipsitch and Murray [18] criteria. Other medical research shows that reinfection rates after successful treatment are much higher than rates of new TB; sometimes approximately four times higher [13,41]. Furthermore, similar to vaccine conferred immunity, the protection rendered by latent TB infection wanes with time and it is uncertain whether latent infections would provide a similar immunity decades after the first episode of TB [43].

Hysteresis
For the usual forward bifurcation (eg. , Fig 4), a model has two locally stable branches at the transcritical point R 0 = 1: i) an infection free equilibrium that is locally asymptotically stable when R 0 < 1 and ii) an endemic equilibrium which is stable for R 0 > 1. However, the scenario where there is only one endemic equilibria when R 0 > 1 may not always be the case. For example, Reluga et al. [52] noted in their study of epidemic models with structured immunity that it is possible that more than one endemic equilibria may coexist even though the basic reproduction number R 0 > 1. This leads to an unusual phenomenon of forward bifurcation with hysteresis, which can be triggered in our TB model when reinfection is taken into account. Thus Eq (1) exhibits a hysteresis effect where multiple endemic equilibria coexist when R 0 > 1, as shown in Fig 8. The two 'outer' equilibria are stable while the interior equilibrium (red) is unstable. Table 2 clarifies that three endemic equilibria coexist for R 0 > 1 if β 0 < β < β 1 . For some parameter regimes with hysteresis, the endemic equilibria may also be found in the region where R 0 < 1 and where disease is not expected, as shown in Fig 9. In this scenario, (similar to a backward bifurcation) TB persists for R 0 < 1 even though the bifurcation at R 0 = 1 is a forward bifurcation. We know of no other epidemic modelling study reporting this feature. Hysteresis loops to the left of the epidemic threshold R 0 = 1 have epidemiological implication in that, although there is no backward bifurcation (as ascertained by the fact that the hysteresis effect occurs when (p < p c )) policy makers need to reduce R 0 below another threshold R c to eradicate TB. That is, reducing R 0 below unity will be necessary but not sufficient in eradicating TB within the community.

No reinfection path A, (i.e p = 0, θ > 0)
The case when there is no reinfection among recovered individuals but reinfection of exposed individuals (i.e p > 0 and θ = 0) was studied by Kar et al. [22], where an exogenous reinfection threshold was established. In this section we focus on the scenario when there is no exogenous reinfection among exposed individuals (i.e p = 0, reinfection path A is omitted). Our goal is to demonstrate that in such a situation recurrent TB due to reinfection of recovered individuals only (θ > 0) can induce the phenomenon of backward bifurcation. Table 2. Generalization of the model equilibria of Eq (1).

Range of R 0 Conditions Equilibria of model system (1) Type of bifurcation
One positive endemic equilibrium Forward bifurcation One positive endemic equilibrium Forward bifurcation β 0 < β < β 1 Two positive endemic equilibria Associated to hysteresis One positive endemic equilibrium Forward bifurcation β β 0 , β β 1 One positive endemic equilibrium Forward bifurcation β 0 < β < β 1 Three positive endemic equilibria Hysteresis In model system (1) setting p = 0 we obtain the following subsystem: The model system (18) equilibrium points (S Ã , E Ã , I Ã , R Ã ) can be expressed in terms of force of infection λ Ã , obtained by solving the following quadratic equation; Note that subsystem (18) has the same basic reproduction number as model (1).
It is easy to see that R 0 = 1 implies b 0 = 0. Thus, the following equality is satisfied; This combined with the condition b 1 < 0, which is necessary for backward bifurcation to occur, and with some algebraic manipulation leads to the criterion for backward bifurcation: Thus, we arrive at the following Theorem; Theorem 2 The model subsystem (18) at R 0 = 1 has: Furthermore, a similar condition to (21) can be obtained by setting p = 0 in the center manifold results of S2 Appendix, hence corroborating Theorem 2. Thus, if θ > θ c = 5 model (18) will exhibit backward bifurcation. However, if θ < θ c = 5 model (18) does not exhibit backward bifurcation (i.e has only forward bifurcation) for the parameters given in Fig 10 caption.
With the existing evidence that recovered individuals have increased susceptibility to reinfection, four times higher than that of new TB [13,27], Theorem 2 suggests that the contribution of recurrent TB in the general TB burden can significantly alter TB dynamics especially in a scenario where recurrent TB independently triggers the phenomenon of backward bifurcation.
It is important to note that although previous TB models have attempted to incorporate recurrent TB pathways they do not investigate whether recurrent TB alone can trigger bi-stability, but rather concentrate on backward bifurcation caused by exogenous reinfection of exposed individuals (i.e p > 0). Selecting value of parameter θ from the estimated interval, i.e θ 2 3.87 [1.61, 7.79] [41], we verify the threshold given in Theorem 2.

Impact of incorporating recurrent TB parameters
Recall that recurrent TB due to reinfection is denoted by reinfection pathways B and C in Fig  3. Given that reinfection path A does induce backward bifurcation when p > p c it is of interest to investigate how the additional recurrent reinfection paths B and C can impact the backward bifurcation. In Fig 13 the endemic equilibrium I Ã is plotted as a function of the basic reproduction number R 0 for scenarios with different recurrent TB contributions. The figure shows that recurrent TB due to reinfection among treated/recovered individuals shifts the backward bifurcation to the left as well as widens the gap between bifurcation curves. In summary we Recurrent TB induces a rare hysteresis effect report the following important implications that arise from our analysis. Raising the intensity of recurrent TB: (i) Widens the gap between the bifurcation curves, (ii) Reduces the critical value R c , (iii) Increases the number of infected individuals (TB burden) (see Fig 14), (iv) Reduces the reinfection threshold p c that induces the phenomenon of backward bifurcation.

Discussion
Until just before the turn of the century the role of exogenous reinfection in the transmission of TB was usually believed to be minimal. van Rie et al. [50] wrote: "For decades it has been assumed that postprimary tuberculosis is usually caused by reactivation of endogenous infection rather than by a new, exogenous infection." However, these views are no longer considered accurate and our understanding of the role of exogenous reinfection has been completely revised. Warren et al. [53] refuted the unitary concept of pathogenesis of tuberculosis proposed in 1960s: that is, tuberculosis results from a single infection with a single Mycobacterium tuberculosis strain, and such infections were thought to confer protective immunity against exogenous reinfection. Thus, exogenous reinfection was thought to be uncommon. Murray et al. [19] found that their data for exogenous reinfection among US immigrants strongly suggested that reinfection likely plays a major role in high-incidence TB areas. Lipsitch and Murray [18] argued that backward bifurcation should not occur when taking into account biologically realistic parameter values. Their argument was built on the premise that exogenous reinfection among exposed individuals should be less than the probability of Recurrent TB induces a rare hysteresis effect progressing to clinically active TB due to endogenous reactivation, which translates in mathematics to p < P L = k/(μ + k). But Feng et al. [23] found that backward bifurcation can only take place if p > P feng > P L . For this reason, Lipsitch and Murray [18] concluded that backward bifurcations are unlikely to be relevant in the context of TB despite non-existence of a data to support their claim. The argument of Lipsitch and Murray [18] was based on the Feng et al. [23] TB model which as mentioned, failed to incorporate key TB pathways such as primary progression and recurrent TB due to reinfection where some recovered/treated individuals revert directly to the infective stage. Yet, these pathways are critical to TB epidemiology, an aspect that was even pointed out by Lipsitch and Murray [18] as a weakness of the Feng et al. [23] model. Lipsitch and Murray [18] called for further research that would account for these omitted pathways. However, modern research no longer seems to support Lipsitch and Murray's [18] argument since exogenous reinfection among individuals with latent TB can outweigh endogenous reactivation [20,21]. And this implies that the protection provided by latent TB infection is not strong enough to prevent individuals becoming reinfected. This is supported by the fact that majority of new TB cases (about ninety percent) are as a result of reinfection rather than endogenous reactivation [13,[19][20][21]. Moreover, studies of TB reinfection in high-HIV burden countries have demonstrated the strong negative impact of HIV on immunity, which is likely to outweigh the possible protection conferred by latent infection [43]. Backward bifurcation is likely to occur in this scenario, since reinfection will be greater than endogenous reactivation (i.e, p > k/(μ + k)).
The analysis conducted here has shown that when exogenous reinfection is significant (resulting in relatively small p c -see S5 Appendix), as is currently understood to be not atypical, backward bifurcation can indeed occur for relatively low values of the reinfection parameter p. Furthermore, we observed that if we omit the exogenous reinfection path A required to cause backward bifurcation in the study of Feng eta al. [23], then recurrent TB alone (i.e., paths B and C in our model) may still yield backward bifurcations. For instance, if p = 0 (thereby omitting pathway A), and the recurrent TB rate parameter θ exceeds a certain threshold θ > θ c then backward bifurcation can indeed occur (see Figs 11 and 12). To the best of our knowledge, this result has not been observed in previous TB modelling studies. In addition, recurrent TB can induce forward bifurcation with hysteresis, rather than just the usual forward bifurcation. Interestingly, the hysteresis loop depends on reinfection parameters, and this may lead to an unusual scenario where the hysteresis crosses the threshold R 0 = 1 thus entering the region where only backward bifurcation is expected (see Fig 9). Epidemiologically this implies that TB will persist when R 0 < 1 even though we have only a forward bifurcation.
In the literature it has been observed that individuals who previously have had active TB and who were successfully treated are more likely to gain active TB another time [13]. However there is little understanding as to why this occurs. Gomes et al. [41] suggested two alternative mechanisms that might explain this: a) previous infection increases susceptibility of individuals to reinfection; b) population heterogeneity, in which some individuals are more atrisk than others, might lead to this conclusion. In this paper, we have almost exclusively explored the former possibility. However the latter possibility may also be at play as examined by Gomes et al. [41]. In the latter case, heterogeneity would be less likely to create changes in reinfection parameters such as θ, and thus might not lead to the same bifurcation phenomena found in this study. Nevertheless, a wealth of theoretical studies do suggest that heterogeneous infection processes are often involved in creating backward bifurcations, and thus we similarly expect complex dynamical phenomena [54][55][56].
In future work, it would be interesting to additionally consider the possibility that some infected individuals who are treated do not become completely cured. This situation would lead to another cohort of individuals characterised by the fact that treatment has failed, and thus require adding an extra compartment which distinguishes complete recovery and incomplete recovery (of individuals who are infectious). The extra infectives from the latter compartment will tend to increase TB prevalence and thus widen the bifurcation curves. While this possibility is of importance, it falls outside the main scope of the present article, and would lead to a model that is very difficult to analyse mathematically.