DYNAMICAL MODELS OF TUBERCULOSIS AND THEIR APPLICATIONS

The reemergence of tuberculosis (TB) from the 1980s to the early 1990s instigated extensive researches on the mechanisms behind the transmission dynamics of TB epidemics. This article provides a detailed review of the work on the dynamics and control of TB. The earliest mathematical models describing the TB dynamics appeared in the 1960s and focused on the prediction and control strategies using simulation approaches. Most recently developed models not only pay attention to simulations but also take care of dynamical analysis using modern knowledge of dynamical systems. Questions addressed by these models mainly concentrate on TB control strategies, optimal vaccination policies, approaches toward the elimination of TB in the U.S.A., TB co-infection with HIV/AIDS, drug-resistant TB, responses of the immune system, impacts of demography, the role of public transportation systems, and the impact of contact patterns. Model formulations involve a variety of mathematical areas, such as ODEs (Ordinary Differential Equations) (both autonomous and non-autonomous systems), PDEs (Partial Differential Equations), system of difference equations, system of integro-differential equations, Markov chain model, and simulation models.

1. Introduction.Tuberculosis (TB) is a disease that affects human and animal populations.Ancient Egyptian mummies show deformities consistent with tubercular decay [20,23].TB was probably transmitted from animals to humans in areas where agriculture became dominant and animals were domesticated.The growth of human communities probably increased the recurrence of TB epidemics leading to its currently overwhelmingly high levels of endemicity in some developing nations.McGrath estimates that a social network of 180 to 440 persons is required to achieve the stable host pathogen relationship necessary for TB infection to become endemic in a community [58].Historically the terms phthisis, consumption and white plague were used as synonym for TB.
TB was a "fatal" disease.In earlier times, some physicians refused to visit the late-stage TB to keep their reputation.TB was responsible for at least one billion deaths during the nineteenth and early twentieth century and the leading cause of 362 C. CASTILLO-CHAVEZ AND B. SONG human death for centuries.Today, "only" 3 million deaths worldwide are attributed to TB every year.World Health Organization's (WHO) data shows that most cases of TB are in developing countries.Twenty three counties in East Asia and Africa account for over 80% of all cases around the world [84].
It was not clear how TB was transmitted until Robert Koch's brilliant discovery of the tubercle bacillus in 1882 (Koch also identified the cause of anthrax).He identified Mycobacterium tuberculosis as the causative agent of TB.The tubercle bacilli live in the lungs of infected hosts.They spread in the air when infectious individuals sneeze, cough, speak or sing.A susceptible individual may become infected with TB if he or she inhales bacilli from the air.The particles containing Mycobacterium tuberculosis are so small that normal air currents keep them airborne and transport them throughout rooms or buildings [83].Hence, individuals who regularly share space with those with active TB (the infectious stage of the disease) have a much higher risk of becoming infected.These bacilli become established in the alveoli of the lungs from where they spread throughout the body if not suppressed by the immune system.
The hosts' immune responses usually limit bacilli multiplication and, consequently, the spread that follows initial infections.About 10% of infected individuals eventually develop active TB.Most infected individuals remain as latently infected carriers for their entire lives.The average length of the latent period (noninfectious stage) ranges from months to decades.However, the risk of progression toward active TB increases markedly in the presence of co-infections that debilitate the immune system.Persons with HIV co-infections progress faster towards the active TB state than those without them [68].
Most forms of TB can be treated.Effective and widespread treatment for active and latently infected individuals has been available for about five decades.Streptomycin is still used today to treat TB but in combination with pyrazinamide. .Data taken from [22].
Isoniazid and rifampin are thought to be the most effective in the fight against M. tuberculosis.The widespread introduction of antibiotics reduced mortality by 70% from 1945 to 1955 in the U.S.A. albeit most major reductions in TB mortality rates had already been achieved before their introduction [4,29,54].Latent TB can be handled with isoniazid but treatment is effective only if applied for at least six months.Active cases must be treated for nine months with multiple drugs (isoniazid, rifampin, pyrazinamide) and complex regimens.Treatment covers over 95% of the cases in the U.S.A. despite its high cost [85].Antibiotic-resistant strains are easily generated when treatment is not completed.The consequences of incomplete treatment may be serious [15].Lack of treatment compliance has serious consequences due to its dramatic impact on the evolution of antibiotic resistant strains [49].The expenses associated with treatment programs for those with active TB are so high that their effective implementation is out of the reach for most developing nations.
As shown in Fig. 1, the mortality associated with TB in the U.S.A. continues to exhibit a downward trend.The annual case rate of TB had been declining steadily but raised slightly in the 1980s and early 1990s in the U.S.A..The change in this trend had been labeled as a period of TB reemergence.TB reemergence over the past decade and a half has challenged existing prevention and control TB programs in developing nations.
In this paper, we review some of the literature associated with TB models and their theoretical impact-particularly those aspects where the authors or their collaborators have contributed.Some results appearing in this paper by Song and Castillo-Chavez have not been published.
The paper is organized as follows.Section 2 introduces the notation that we try to use throughout the manuscript.Section 3 reviews some of the earliest known TB models.Section 4 deals with the exploration of the impact of various epidemi- Data taken from [78,79].
ological factors as well as the role of close and casual contacts on TB dynamics.section 5 looks at the impact of demography on TB dynamics.In section 6, we review some cell-based models for TB transmission at the immune system level.A Markov chain model on TB projections is described in section 7. Models dealing with TB control strategies are discussed in section 8.A model dealing with the role of public mass transportation on TB evolution and control is reviewed in section 9.
Finally, a list of challenges associated with modeling TB dynamics is outlined.
2. Notation.The population of interest is divided into several compartments (classes, categories, or subpopulations) dictated by the epidemiological stages (host statues).For the most part, in the context of TB, four or five epidemiological stages are identified (see Table 1).We shall do our best to denote these subclasses using uniform symbols as we discuss a multitude of models.Tables 1 and 2 list the definitions and symbols (subpopulations and parameters) that we try to use.
If it is necessary to subdivide a population into subpopulations, subscripts will be V Vaccinated possibly reduced susceptibility to TB used to distinguish them.For instance, I s and I r represent the drug-sensitive and drug-resistant infectious TB classes, respectively.Active TB, case TB, index TB, mature TB, open case, and lesion case all mean active-TB infectious case here.We use β to measure of the likelihood of transmission or the "force" of infection.However, the meaning of β or its interpretation often changes from model to model.

3.
Early dynamical models.The first model for the transmission dynamics of TB was built in 1962 by Waaler [81], the chief statistician of the Norwegian Tuberculosis Control Services.Waaler divided the population into three epidemiological classes: noninfected (susceptible), infected non-cases (latent TB), and infected cases (infectious).He formulated the infection rate as an unknown function of the number of infectious individuals.He used a particular linear function to model infection rates in the implementation of his model.The incidence (new cases of infections per unit time) was assumed to depend only on the number of infectious.Furthermore, the equations for the latent and infectious classes were assumed to be uncoupled from the equation for the susceptible class.The central part of this model is given by the following linear system of difference equations: where the incidence rate aI t is proportional to the number of infectious; e is the per-capita progression rate from latent-TB to infectious-TB cases; g is the percapita treatment rate (treated individuals will become members of latent-TB class again.);d 2 is the per-capita death rate of the latent-TB class; and d 3 is the percapita death rate of the infectious-TB class.Using data from a rural area in south India for the period of 1950 to 1955, Waaler [35] estimated the parameters of this linear model to be a = 1, e = 0.1, d 2 = 0.014, g = 0.10085, d 3 = 0.07.Because the eigenvalues all have norms close to 1 (1.04),Waaler predicted that the time trend of TB is unlikely to increase (it may decrease, albeit slowly).This linear model did not model the mechanics of transmission.However, the parameters, estimated from a specific area in India, set useful ranges for the estimation of parameters in developing nations.
Brogger developed a model [10] that improved on Waaler's.Brogger not only introduced heterogeneity (age) but also changed the method used for calculating infection rates.The infection rate in Brogger's model was a combination of linear and nonlinear infection terms.In fact, it was given by the term βS(1 − Z + Z I N ), where Z was an adjusting parameter used to differentiate between normal infection, superinfection, and direct leaps (within a very short period, an uninfected individual becomes a lesion case or an active-TB case).Two extreme cases were covered in the model: Z = 1 making the incidence be βS I N , the familiar version of today, and, Z = 0 giving an infection rate proportional to the number of susceptibles.The prevalence I N was used to adjust all flow rates including those from infected to open cases.This was not surprising as Brogger wanted to use prevalence as an indicator of the effectiveness of control policies.His aim was to compare different control strategies that included finding and treating more cases, the utilization of vaccination, and mass roentgenograph.The data of two WHO/UNICEF projects in Thailand from 1960 to 1963 were used to estimate the parameters incorporated into his model.Brogger chose those parameters that "best" fits available data.Control strategies (additional new parameters) were "squeezed" into the model.Simulations were run and comparisons made to evaluate the value of different strategies of TB control programs.This model did not formulate clearly the relationship between infection rate and prevalence.ReVelle classified this important relationship in 1967.
Using Brogger and Waaler's model as a template, ReVelle introduced the first nonlinear system of ordinary differential equations that models TB dynamics [63,64].In modeling the infection rate, he did not follow the typical mass action law, given by the bilinear function βSI of Kermack and McKendrick [47].It was ReVelle who first, at least in the context of TB dynamics, rigorously explained why the infection rate depends linearly on the prevalence using the probabilistic approach that is common today (homogeneous mixing).The form βS I N for the infection rate is found in most epidemic models used today.Mathematically it is well known that if the total population size N remains constant over time or if it asymptotically approaches a constant then the use of an infection rate proportional to SI does not change the qualitative properties of the model.However, when modeling epidemics for developing countries, as Revell did with his model, βS I N seems a more appropriate form of modeling the infection rate.ReVelle modeled TB dynamics via a system of non-linear differential equations but he ignored population structure.Nine compartments were introduced in Revell's nonlinear model.The total population was governed by the Malthus model because he wanted to apply it to developing nations.Making projects (in Waaler's words "time trend of tuberculosis") was not Revell's main theme.In fact, his main objective seemed to be associated with the evaluation and implementation of control polices and their cost.He developed an optimization model and used it to select control strategies that could be carried out at a minimal cost.It is worth to mention that Waaler also developed a model in 1970 that would minimize the cost of alternative tuberculosis control measures [82].
All dynamical models prior to Ferebee's work were motivated by the study of TB in developing nations.The two data sets used were from Thailand and India.No specific model seems to have been developed for the U.S.A.. Ferebee, associate chief of the research section of the Public Health Service Tuberculosis Program of U.S.A., changed this trend.She set up a discrete model, based upon a set of simple assumptions, to model the dynamics of TB in the U.S.A. [34].She used the same compartments as Waaler did, that is, susceptible, infected and infectious.The basic time unit was a year.Hence, within one year, new infected people would become infectious and contribute to the pool of new infected, namely, some individuals would move from the susceptible to the infected class and from the infected to the infectious class within a year.Ferebee described her algorithm, methods of estimation of relevant parameters, and the number of infected people in the U.S.A..The results showed that the number of new cases would decrease, but slowly, if vaccination were not applied to the US population.It is worth mentioning that the estimation of demographic parameters was solely based on the 1963 US data.Despite its shortcomings, this work indeed gave the first rough estimate and forecast of TB cases in the U.S.A..She stated that her assumptions were checked for consistency with bits and pieces of information obtained from a variety of sources.These assumptions have since played an important role theoretically and practically in the context of TB.In other words, the first US study "defined" an appropriate parameter range.We shall outline underlying assumptions: i.There are 25 million infected individuals and 125 million susceptible; ii.The per-capita progression rate from infected to infectious is k 1 = 1/625 per year; iii.Primary infected people exhibit a higher per-capita progression rate in the first year (k 2 = 1/12); that is, one out of every 12 new infections will progress to the infectious stage during the first year; iv.Each new infectious individuals will infect 3 people; v.No significant additional death is ascribed to tuberculosis.
Her assumption k 2 >> k 1 has directed model simulations and constructions.This assumption has been recently incorporated via the inclusion of additional compartments for fast TB and slow TB (see section 4.1).
Earlier mathematical models for TB transmission were developed mostly by statisticians.Their approach followed a pattern: build a mathematical model for TB transmission; with help from a data set, estimate parameters; find numerical solutions; and predict or make inferences about the relative value of alternative control strategies.There was no qualitative analysis of the models, and the long-time behavior (asymptotic properties) of the models was also not studied.
The continuous decline of TB incidence in developed nations and the introduction of effective antibiotics suggested that elimination of active TB in developed nations was possible.This view may have been the main reason why there was almost no theoretical work on TB dynamics from the 1970s to the early 1990s.The story has changed over the last decade because of the reemergence of TB (new outbreaks in the U.S.A. and in many developed nations).In the following sections we shall review some of the most recent models and the theoretical results.

Intrinsic mechanics of transmission.
4.1.Slow and fast routes.The initially exposed individuals (infected individuals) have a higher risk of developing active TB.With time passing, those individuals still face the possibility of progressing to infectious TB, but the rate of progression slows down.In other words, the likelihood of becoming an active infectious case decreases with the age of the infection.Bearing this in mind, several researchers constructed a series of dynamical models for TB progression and transmission in scenarios that took these factors into consideration [6,7,8,19,32,60].We shall review some of this work.In the simplest model (that we know), the population of interest is partitioned into three epidemiological classes: susceptible, latent, and infectious.The infection rate given by βSI (using the mass action law) is divided.A portion pβSI gives rise to immediate active cases (fast progression), while the rest (1 − p)βSI gives rise to latent-TB cases with a low risk of progressing to active TB (slow progression).The progression rate from latent TB to active TB is assumed to be proportional to the number of latent-TB cases, that is, it is given by kE, where k ranges from 0.00256 to 0.00527 (slow progression).The total incidence rate is pβSI + kE.The version in [6] is given by following system: where the parameters are defined in Table 2.The qualitative dynamics of model (3)(4)(5) are governed by the basic reproductive number This dimensionless quantity measures the average number of secondary infectious cases produced by a "typical" infectious individual in a population of susceptibles at a demographic steady state.The first term in (6) gives the new cases resulting from fast progression while the second those resulting from slow progression.Sensitivity and uncertainty analysis were carried out.Simulation results showed that TB dynamics were quite slow for acceptable parameter ranges.Waaler's model also supported slow TB dynamics [81].Model (3)(4)(5) requires that p be known a priori (it is not allowed to change) and, the R 0 derived from the model depends linearly on population size.A model that removes these restrictions is reviewed next.e −(µ+r1)(τ −s) − ṗ(τ − s)e −(µ+r2+d)(t−τ ) dsdτ.
The following system of integro-differential equations is used to model TB dynamics with a variable latent period: where σ = βc is the force of infection per infective;   (10)(11)(12)(13)(14)(15) when q = 0.In region I, the disease-free equilibrium is globally asymptotically stable; in regions II and IV, one of strains disappears; and III represents the coexistence region.(10)(11)(12)(13)(14)(15) when q > 0. The coexistence of the two strains is impossible.
The introduction of an arbitrary distribution of latency period did not change the qualitative dynamics of TB; that is, a forward bifurcation diagram characterizes the dynamics of the last model (7-9).4.3.Multiple strains.Incomplete treatment, wrong therapy, and co-infection with other diseases, for instance, HIV, may give rise to new resistant strains of TB (multiple drug resistant, or MDR strains).Form First Lady Eleanor Roosevelt was one of the victims of MDR TB [24,62].Models that include multiple strains of TB have been developed [7,15,17].A recently published two-strain TB model include drug-sensitive and drug-resistant strains [15,17].Hence, two subclasses of latent and infectious individuals are required.The subscripts s and r stand for drug-sensitive and drug-resistant types.The model is given by the following set of equations: It can be seen from Equation (15) that the treatment rate for the I r class is equal to zero, meaning that TB due to this strain is not treatable by current antibiotics.The proportion of treated infectious individuals who did not complete treatment is p+q.The proportion p modifies the rate at which they depart from the latent class; qr 2s I s gives the rate at which individuals develop resistant-TB due to their lack of compliance with TB treatment.The proportion of successfully treated individuals is 1 − p − q.The dimensionless quantities give the basic reproductive number of strains j, where j = 1, 2. The asymptotic behavior of model (10)(11)(12)(13)(14)(15) is determined by R j .Fig. s 3 and 4 show the bifurcation diagram for this two-strain model.These diagrams shows that naturally resistant and natural types can co-exist, albeit the region of coexistence is small (region IV in Fig. 3).Furthermore, in [15] it was shown that antibiotic-induced resistance results in the substantial expansion of the region of coexistence.In fact, regions IV and III become a single large region of coexistence.In other words, antibioticinduced resistance guarantees the survival of resistant strains.A two-strain model that incorporates the effects of multiple drug resistance can also be found in [7].4.4.Multiple strains and variable latent period.A model that considers both multiple strains and variable latent period is proposed by Feng [33].Drug-resistant and drug-sensitive strains are modeled, but only the age of the infection with drugsensitive strain is considered.They introduce a function p(θ) as the proportion of the sensitive-strain that are active at infection-age θ to distinguish active TB and inactive TB.The model does not make a difference between active and inactive TB for the drug-resistant strain because after acquiring drug-resistent TB, an individual dies quickly.Consequently for the drug-resistant strain they only count active TB.The total population is divided into three classes: susceptible (S(t)), infections with drug-sensitive strain (I s (t)), and infections with drug-resistant strain (I r (t)).Letting i s (θ, t) be the infection-age density of infected individuals with the drugsensitive strain at time t, the model framework takes the following system: where  (20)(21)(22)(23).When R 0 < R p , the diseasefree equilibrium is globally asymptotically stable.However, when R p < R 0 < 1, there are two endemic equilibria.The upper ones are stable, and the lower ones are unstable.
birth rate into the population; N (t) = S(t) + I s (t) + I r (t) the total population size; The disease-induced mortalities for drug-sensitive TB and drug-resistant TB are denoted by d i respectively.It is assumed that a fraction r of the treated individuals with the drug-sensitive strain does not recover due to incomplete treatment, and the remaining fraction 1 − r is successfully treated and become susceptible again.
It is also assumed that a fraction q of those who do not finish their treatment will generate drug-resistant TB and the remaining fraction 1 − q of them will keep as infectious.
Model [16][17][18][19] shows "that nonantibiotic-induced coexistence is possible but rare for naturally resistant strains, while coexistence is almost the rule for strains that result from the lack of compliance with antibiotic treatment by TB-infected individuals" [33].4.5.Exogenous reinfection.Immediately after primary infection, an infected individual, on average, has a higher risk of progression (becoming an infectious case).Infected individuals who do not become infectious within a short time period may still develop active TB via exogenous or endogenous reinfection or both.To see the impact of exogenous reinfection on the dynamics of TB, we discuss the model (see [32]) given by the following set of equations: The term pβcE I N models exogenous reinfection, that is, the potential reactivation of TB by continuous exposure of latently-infected individuals to those who have active infections.
The dynamics generated from model (20)(21)(22)(23) are "surprising" as they show that R 0 = 1 is not always the key threshold.We introduces a method that deals with bifurcation problem arising from this model as well as more general epidemic models in the following subsection.4.6.An approach to determine the direction of the bifurcation at R 0 = 1.For the most part, in epidemic models, there are two distinct bifurcations at R 0 = 1 forward (supercritical) and backward (subcritical).A forward bifurcation happens when R 0 crosses unity from below; a small positive asymptotically stable equilibrium appears and the disease-free equilibrium losses its stability.On the other hand, a backward bifurcation happens when R 0 is less than unity; a small positive unstable equilibrium appears while the disease-free equilibrium and a larger positive equilibrium are locally asymptotically stable.Epidemiologically, a backward bifurcation "says" that it is not enough to only reduce the basic reproductive number to less than one to eliminate a disease and that when R 0 crosses unity, hysteresis takes place.This phenomenon that probably was first found in epidemiological models by Huang et al. [45] in 1992 in a study of an HIV/AIDS model has now become common.(Huang et al. identified it as a subcritical bifurcation.)Recent studies with models supporting backward bifurcations include those of Hadeler and Castillo-Chavez in a variable core group model [41]; Dushoff et al. [30] in models for fatal diseases; Feng et al. [32] in a TB model with exogenous reinfection; and Kribs-Zaleta and Velasco-Hernandez [50] in a model with vaccination.Most recently, van den Driessche and Watmough [27,28] have shown the existence of this behavior in epidemic models with delay while Castillo-Chavez and Huang have done the same for models with age-structure [18].
Center manifold theory has been used to decide the local stability of a nonhyperbolic equilibrium (linearization matrix has at least one eigenvalue with zero real part) [14,37,86].We shall describe a theory that not only can determine the local stability of the nonhyperbolic equilibrium but also t settles the question of the existence of another equilibrium (bifurcated from the nonhyperbolic equilibrium).This theory is based on the general center manifold theory.To describe it, consider a general system of ODEs with a parameter φ: Without loss of generality, it is assumed that 0 is an equilibrium for System (24) for all values of the parameter φ, that is 24) around the equilibrium 0 with φ evaluated at 0. Zero is a simple eigenvalue of A and all other eigenvalues of A have negative real parts; A2: Matrix A has a nonnegative right eigenvector w and a left eigenvector v corresponding to the zero eigenvalue.Let f k be the kth component of f and The local dynamics of (24) around 0 are totally determined by a and b. i. a > 0, b > 0. When φ < 0 with |φ| ≪ 1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when 0 < φ ≪ 1, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium; ii. a < 0, b < 0. When φ < 0 with |φ| ≪ 1, 0 is unstable; when 0 < φ ≪ 1, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; iii. a > 0, b < 0. When φ < 0 with |φ| ≪ 1, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when 0 < φ ≪ 1, 0 is stable, and a positive unstable equilibrium appears; iv. a < 0, b > 0. When φ changes from negative to positive, 0 changes its stability from stable to unstable.Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.
The results of Theorem 4.1 are summarized in Table 3. Huang et al. [45] and Dushoff et al. [30] initially decomposed the center manifold and found out that the sign of a can be used to determine the direction of the bifurcation.In doing so, they computed the Taylor expansion of f (x, φ) only on its state variables.Van den Driessche and Watmough [28] used a Taylor expansion on the state variables and the parameter φ in the analysis.Let us sketch our proof below.
Proof.Let E c and E s be the generalized eigenspaces of A for the zero eigenvalue and all other eigenvalues, respectively.It follows from the center manifold theory that the center manifold W c is one dimensional, and R n = E c ⊕ E s .We parameterize the center manifold by c(t) and decompose it into E c and E s , that is, where c(t) ∈ E c and h(c, φ) ∈ E s .Because the center manifold is tangent to E c at the origin, h(c, φ) is a higher order term (h(c, φ) has at least order 2).It also follows by the invariance of the center manifold under the flow that Applying Taylor expansion to the right-hand side of ( 29) at (0, 0) and noticing that h(c, φ) is of higher order, we obtain that where D 2 xx is the Hessian matrix; I n is the identity matrix of order n; ⊗ is the Kronecker product.Using φφ f (0, 0) = 0 and the fact that ch(c, φ) is of higher order, we simplify the above expansion for f as (higher order terms are dropped) Multiplying both sides of (29) by v and using the fact that v•h = 0 and vD x f (0, 0) = 0, we finally obtain the following equation for c(t): Obviously, at φ = 0 a transcritical bifurcation takes place in Equation (31).The bifurcation diagram is shown in Table 3.
Table 3. Summary of Theorem 4.1.In the bifurcation diagrams, the vertical axis represents equilibrium points x * , and the horizontal axis is the parameter φ.Solid lines and dashed lines symbolize stable (S) and unstable (U), respectively.
a and b stability of 0 stability and sign of The requirement that w is nonnegative in Theorem 4.1 is not necessary.When some components in w are negative, we still can apply this theorem, but one has to compare w with actual the equilibrium because the general parameterization of the center manifold before the coordinate change is provided that x 0 is a nonnegative equilibrium of interest (Usually x 0 is the diseasefree equilibrium).Hence, x 0 − 2bφ a > 0 requires that w(j) > 0 whenever x 0 (j) = 0.If x 0 (j) > 0, then w(j) need not be positive.Applying Corollary 4.1 to model (20)(21)(22)(23), we give a rigorous proof that model (20)(21)(22)(23) undergoes a backward bifurcation.We do it in a rather simple case σ = 1 though when σ = 1 the arguments are identical.Let φ = βc be the bifurcation parameter.Introducing with . The linearization matrix of system (32)(33)(34) around the disease-free-equilibrium when φ = φ * is It is clear that 0 is a simple eigenvalue of D x f .A right eigenvector associated with The rest of the second derivatives appearing in the formula for a in (26) and b in ( 27) are all zero.Hence, We collect these results in the theorem below: Theorem 4.2.If p > p 0 = k µ 1 + k µ+r+d , the direction of the bifurcation of system (32-34) (or system (20)(21)(22)(23)) at R 0 = 1 is backward.
Model (20-23) exhibits totally different behavior as it supports multiple stablesteady states via a backward bifurcation (subcritical bifurcation).The results are counter-intuitive.that is, when the basic reproductive number R 0 = βc µ+r+d k µ+k < 1, the disease-free equilibrium is only locally stable.The system supports multiple endemic equilibria (see Fig. 5).The results of this model reveal how hard it may be to eliminate TB if exogenous reinfection is important.Vynnycky and Fine have also incorporated reinfection in their non-autonomous and linear PDE model (see section 8).4.7.Generalized households.Both close and casual contacts may give rise to infections.Individuals having lengthy exposure may have a substantially higher risk of infection than those having only casual contacts with active-TB cases.For thousands of years, no one knew that TB was an infectious disease.In the nineteenth century, doctors in western Europe and the United States thought that TB was hereditary because it often ran in families.After Doctor Koch identified the pathogen responsible for TB, it became clear that prolonged contacts with an active-TB case might generate new infections.In fact, now it is known that most new infections are due to close contacts.The case of a teacher-librarian with active TB who infected the children in her classroom but not the children who visited the library [61] supports the incorporation of differences between casual and close contacts in the investigation of TB dynamics at the population level.The report of the Asian tourist who infected at least 6 airline passengers while travelling to and fromand within-the United States further supports this view [49,46].In fact, in 1998, WHO declared that flights more than 8 hours may expose passengers and crew to TB infections [59].Hence, distinguishing the type of contacts in the transmission of TB is important.Aparicio and colleagues proposed a new dynamical model that incorporated close and casual contacts [2].They focused on the active-TB cases within their social networks (family members, officemates, classmates, any persons who have prolong contacts with an index case).This "new" epidemiological unit, called the generalized household or the epidemiologically active cluster, was used to study the transmission dynamics of TB.We describe the approach next.The population is partitioned into TB-active clusters and TB-inactive clusters.Once one member of an inactive cluster develops active TB, the whole cluster becomes an epidemiologically active-TB cluster.Specific contacts within each epidemiologically active cluster are not modeled, it is assumed that the risk of infection within an epidemiological active cluster just depends on the "life" of the cluster, that is, on the average length of TB's infection period distribution.The flow chart for the basic cluster model is shown in Fig. 6.Letting N 2 = S 2 + E 2 , the model equations are The basic reproductive number for this model is It can be seen that R 0 depends nonlinearly on the parameter β (risk of infection on an epidemiologically active cluster of size n) and linearly on the average generalized household size, n.Song et al. [71] further developed the basic cluster model (35)(36)(37)(38)(39) via the incorporation of casual contacts.An extended version of the cluster model for the TB transmission that includes the close and casual contacts is given by the following nonlinear system: Here β is assumed to depend on the average cluster size n; p denotes the average fraction of time spent by the source case within his/her generalized household; and 1 − p, the average fraction of time spent by this source-case outside the cluster.The rate of infection within clusters becomes pβ(n)S 1 , while the rate of infection outside is , where N is the total population size, and N −n represents the average total number of individuals outside the cluster.Hence, (1−p)β * I N −n S 1 gives the number of new infections per unit time in the N 1 population, that is, the incidence from S 1 to E 1 .The term (1 − p)β * I N −n S 2 gives the incidence from S 2 to E 2 .There are no new cases of active TB within each epidemiologically active cluster, and consequently, the infection rate is pβ(n)S 1 .The generation of new cases from casual contacts is modeled in the typical way.When p = 1 and β(n) is a constant, the extended cluster model  becomes the basic cluster model (35)(36)(37)(38)(39).The basic reproductive number associated with the generalized cluster model is where K = Λ µ is the asymptotic carrying capacity of the total population.Two specific forms of β(n) are discussed below: Case 1.In the first case, we consider that β(n) is a piecewise-defined continuous function that is, where n L is the critical cluster size after which β(n) begins to decrease, and n M is an upper bound for cluster size.Consequently, for n ≤ n L , for n ≥ n L .
Since the maximum cluster size n M is significantly smaller than the carrying capacity K, then K K−n ≈ 1, whenever n < n L R 0 increases linearly with cluster size, and Hence, an increase in n translates an increase in TB transmission but the increase is non-linear.Initially, this increase is almost linear but as n becomes larger the rate of increase decreases, because the time spent by an infectious individual per contact cancels out increases in cluster size.Hence, R 0 (n) is bounded by a constant value, and this bound limits the size of TB prevalence.
Case 2. The second case assumes β(n) to be inversely proportional to n (that is, β(n) = β1 n ).It can be seen from 45 that R 0 (n) is the sum of contributions from both the within and out of cluster new secondary cases of infection.The ratio E(n) of within to between cluster contributions is given by Function E(n) increases with n and reaches its maximum value at Hence,n * here defines the optimal cluster size; that is, the value maximizes within cluster transmission.
Singular perturbation theory and multiple time scales techniques are used to study the global dynamics of the cluster models [71].Since the average infectious period (about 3-4 months) is much shorter than the average latent period, which has the same order as the host population, two different time scales can be identified.Time is measured using the average latency period 1/k as the unit of time (that is, τ = kt).Variables S 2 and E 2 are re-scaled by Ω, the total asymptotic population size (Ω = Λ µ ); N 1 could be re-scaled by Ω instead, it is re-scaled by the balance factor The re-scaled model equations are given by where ǫ = k β+γ , m = β β+γ < 1 and B = µ k .The terms in right-hand side of system (46)(47) all have the same order of magnitude whenever ǫ ≪ 1.Therefore, y 1 , y 2 , and y 3 are fast variables and x 1 and x 2 are slow variables.Hoppensteadt's early theorem [44] help us to show the global stability of the endemic equilibrium when ǫ is small.A Liapunov function and a Dulac function are selected to establish the global stability of the disease-free equilibrium.Hence, a global forward bifurcation (see Fig. 7) characterizes the dynamics of the cluster models.5. Models with density dependent demography.Demography plays an important role in the transmission dynamics of TB since the average rate of progression from infected (non-infectious) to active (infectious) TB is very slow.In fact, the average latent period has the same order of magnitude as the life-span of the host population.Two distinct demographic scenarios are studied.In the first, exponential growth is observed over a long time scale, and in the second exponential growth is observed over a short time scale (quasi-exponential growth).
Consequently, three different recruitment rates are used in the study of dynamical models for TB transmission with demography: constant recruitment rate (Λ) [15], linear recruitment rate (rN ) [74], and logistic recruitment rate (rN (1 − N/K)) [74].The general model is given by where the recruitment rate B(N ) includes the three demographies described above, namely, rN , rN (1 − N/K), and Λ.The basic epidemiological reproductive number is given by However, this non-dimensional number is not enough to characterize the dynamics of model (51-54).
We shall consistently use the following compressed notation m r = r + r 2 + d, n r = r + r 1 + k, m µ = µ + r 2 + d, n µ = µ + r 1 + k, and σ = βc to simplify the discussions.First, we study the dynamics of model (56-58) with B(N ) = rN .That is, it is assumed that the total population exhibits exponential growth in the absence of TB (the net growth rate of the population, in the absence of the disease, is r − µ).Total population size increases exponentially if r > µ, and remains constant if r = µ.The case where r < µ is trivial.Hence we assume that r ≥ µ.In the presence of TB, the total population may (theoretically) decrease exponentially even when r > µ if d is large enough; that is, technically, a fatal disease can impact population growth (see also May and Anderson,[57]; Busenberg and Hadeler [11]).Realistic examples of situations where a disease has or is likely to to have an impact on the demographic growth can be found in the work on myxomatosis [51] and on HIV [1,56].Three non-dimensional threshold parameters R 0 , R 1 , and R 2 provide a full characterization of the possible dynamical regimes of system (56-58) under the various demographic regimes.The basic reproductive number gives the average number of secondary cases produced by a typical infectious individual during his/her entire life in a population of mostly susceptibles.It is implied from R 0 < 1 that the infected populations goes to zero, while R 0 > 1 implies that the infected populations grows (initially) exponentially (together with the total population N ).In this last case, there are two possibilities: N grows faster than I, or N does not grow faster than I.In the first case, the fraction u = I N approaches zero as time increases, and the additional threshold parameter plays a role; that is, R 1 discriminates between the last two possibilities.If R 1 < 1, then lim t→∞ u = 0, while R 1 > 1 implies that lim t→∞ u > u * > 0. Our assumption r > µ implies that R 0 > R 1 is true.If the infectious (I) population changes faster than the total population (N ) then (a fatal) disease can drive the population to extinction (even when R 1 > 1).The threshold parameter that decides this last situation is given by where u * is a positive constant (independent of µ (see ( 67)); that is, R 2 determines whether or not the total population size grows exponentially.In fact, population size would decrease exponentially (from TB) only if R 2 < 1. System (56-58) is homogeneous of degree one, and hence, it can support exponential solutions [39,40].However, we are interested in the global dynamics rather than the linear (local) stability of the homogeneous systems.Global analysis requires the rewriting of system (56-58) using the projections u = I N , v = E N .The equations for u, v are given by the following quadratic system: Note that both u and v are independent of N and µ.The subset To further simplify the quadratic system (62-63), we introduce the new variables x and y and re-scale time t.Specifically, we let The re-scaled system is where In the new system, Ω becomes which is positively invariant under the flow of system (65-66).Transformation 64 not only reduces the number of parameters but, more importantly, it fixes the horizontal isocline and decomposes the vertical isocline into a degenerate quadratic curve.Under the standard classification of Ye et al. [87], system (65-66) is a quadratic system of the second type.The following two theorems characterize the dynamics of system (65-66), and hence of system (62-63).
The standard classification of planar quadratic differential systems rules out the existence of closed orbits or limit cycles.(Other approaches can be used to draw the same conclusion, for example, see Busenberg and van den Driessche [12]; Lin and Hethcote [53].The full structure of system (56-58) is provided in Theorem 5.2 below: Theorem 5.2.Consider system (56)(57)(58) and assume that r > µ: i.If R 0 < 1, then (∞, 0, 0) is globally asymptotically stable. ii.
) is globally asymptotically stable and Hence, whenever R 0 < 1 the disease dies out while the total population increases exponentially.Although the disease spreads (that is, the number of infected grows in total number) when R 1 < 1 < R 0 , the proportions I N and E N approach zero.One can observe from the case 1 < R 1 < R 0 that disease-induced mortality can lead to the extinction of a population that would otherwise increase exponentially (a fatal disease can indeed regulate a population).Note that R 2 is positive since u * is positive and independent of µ.We also have established that when r < µ, (0, 0, 0) is globally asymptotically stable even though lim = 0. Note that R 1 < R 0 whenever r > µ.Theorem 5.2 provides a complete characterization of the dynamic structure of model (56)(57)(58).The global dynamics are shown in Fig. 8. Here, we provide the proof that the disease-free equilibrium is a global attractor.Proof.The disease-free equilibrium is ( r−µ r , 0, 0).It is straightforward to show that the endemic equilibrium is unique whenever R 0 > 1 and R * 2 > 1 and the diseasefree equilibrium is locally stable whenever R 0 ≤ 1.Here, we only need to establish the global stability of the disease-free equilibrium under the assumption R 0 This actually produces a differential inequality on the function f (t); that is, It follows that lim > 0 and R 0 < 1.If R 0 = 1, f (t) no longer decays exponentially, but it still vanishes as time goes to infinite.We estimate the derivative of f (t) again as follows: Hence, f (t) is a decreasing function and If liminf t→∞ I(t) > 0, then lim .
For system (51)(52)(53)(54), if R 0 ≤ 1, the disease-free equilibrium is a global attractor; if R 0 > 1 and R * 2 > 1, there exists a unique endemic equilibrium that is globally asymptotically stable under some assumptions.Hence, a global forward bifurcation describes the dynamics if R * 2 > 1 (see Fig. 7).To show the global stability of the endemic equilibrium, an equivalent monotone system to the original one was identified and a strong version of the Poincaré-Bendixon theorem applied.Further when the disease dies out, the decay is of exponential form with a rate proportional to R 0 − 1. Particularly, interesting dynamics are observed when 1 < R 0 and R * 2 < 1.Under these conditions, only the disease-free equilibrium is supported.One cannot study the stability of the disease-free equilibrium by the regular linearization approach since the vector field at the origin is not analytic.Simulations, however, strongly suggest that the disease-free equilibrium is a global attractor.A summary of the results obtained for this model can be found in Table 4. Feng et al. [31] obtained the global forward bifurcation in the case where the recruitment rate is a constant Λ.

6.
Cell-based models: Battle within a host.The study of TB dynamics within an individual host is also within the realm of epidemics on populations of cells.We review some of these models in a rather superficial way as our work may find applications in immunology.
Co-infection with HIV.Co-infection of TB and HIV has been considered to be the main cause for the reemergence of TB in the U.S.A..Here we review a model for the dynamics of TB and HIV at the cell-level.Four populations are considered in a host who has HIV and TB.The model assumes that the compartments are the lymph tissues [48].Letting T (t) be the armed CD4 + and CD8 + T cells; M (t) the macrophage population; H(t) the HIV population; and T b (t) the M. tuberculosis population, the model equations of Kirschner [48] are where e 1 is the source rate of T cells in the absence of infection; e 2 is the rate of change in the T-cell supply; M 0 represents the equilibrium value for the macrophage population; N 1 and N 2 are the numbers of new viral particles produced by an infected T cell and macrophage, respectively; γ 3 is the rate at which CD8 + T cells kill the virus; γ 4 is the rate at which macrophages kill the virus; T cells clear M. tuberculosis at rate γ 5 ; macrophages clear M. tuberculosis at rate γ 6 .In constructing this model, the growth rate of T cells is of Michaelis-Menten type, where µ m is the natural death rate of the macrophage population and g the bursting rate that results from the proliferation of bacilli inside the cytoplasm of a macrophage.If an inactive macrophage dies of natural causes, it releases N bacilli (on average); if it dies of bursting, then g bacilli are released.Usually, we assume g ≥ N .Encounters of T cells with inactive macrophages stimulate the replication of T cells, which is modeled by the term rM 1 T .An exponential growth rate for bacilli is assumed in the absence of macrophages.The total dynamical analysis for this model has not been completed.

Geographic distribution.
A stochastic Markov chain model was proposed to study the impact of various geographic factors.Debanne et al. [25] constructed a multivariate Markovian model to project the distribution of TB cases across the U.S.A. by state, race and ethnic group.In their simulations, the transition probabilities from susceptible to infected are estimated via the annual risk for infection.Similar to the basic reproductive number, a contagious parameter is used to define the average number of susceptible infected by an infectious case [75,76].The effect of HIV is considered as an exogenous input.Because there is a high degree of uncertainty for the value of most parameters (only lower and upper bounds of parameters are available), a sensitivity analysis is conducted.The selection of parameters are used mostly to calibrate the model to data.The initial state (1980) was computed using inverse iteration (difference equations) back to the time when the actual data were available.The same model framework was applied to different racial and ethnic groups and different states by choosing appropriate parameters.Simulations results showed that TB cases change from state to state as well as among racial and ethnic groups.These simulations were put together and gave a picture of TB trends for the whole nation.Results showed that from 1980 to 2010 TB cases in the U.S.A. experienced an intermediate increase followed by a continuing decline.For instance, by 2010 the TB rate would be 4.6 per 100,000.This figure is higher than the expected case rate required to meet the CDC's goals.
8. Control strategies.A variety of dynamical models are proposed to study TB vaccination and elimination strategies.Our new results are included in subsections 8.3 8.4, and 8.5.

8.1.
Vaccination.A TB vaccine called BCG (Bacillus of Calmette and Guerin) has been available for many decades, but its effectiveness in preventing TB is controversial [67].Results of field trials of the vaccine have differed widely.Some indicates protection rates as high as 70% to 80%, while others indicates the vaccine is completely ineffective in the prevention of TB.We review three models with vaccinations, including an age-dependent model.Two vaccine models, preexposure and postexposure, are proposed by Lietman and Blower [52].In constructing their models, the susceptible and latent classes are divided into unvaccinated and vaccinated subclasses, specified by subscripts u and v, respectively.Preexposure vaccinations only vaccinate newborns at a proportion C with probability of success q.The immunity provided by vaccination is not permanent, leading to an expected effective vaccination time ω that should be less than the average life-span.The formula for fast progression is the same as in model (3-5) and slow progression rates in this model are assumed to be dependent on the number of infectious via two general functions k u (I) and k v (I).The model equations are Postexposure vaccinations only apply to the individuals with latent TB.A third subclass of the latent-TB class is added to represent the waned state of vaccination for those with latent TB.The model equations for postexposure vaccinations are An interesting property of models (74)(75)(76)(77)(78) and (79)(80)(81)(82)(83) lies in their ability to identify distinct vaccination strategies for different nations that wish to eliminate TB.For developed countries, where the prevalence of latent TB is low, only a preexposure vaccine with treatment of active TB would be necessary.For developing nations where the prevalence of latent TB is high, a combination of preexposure and postexposure vaccine with treating active TB would be the most effective strategy.One must with care deal with these results, particularly for developing countries with rapidly expanding populations, as the models have no real demography.BCG vaccination is age-dependent.The impact of BCG vaccination may be better assessed with the use of age-structured models.In fact, age has been included in most simulation models [10,25,34].A linear partial differential equation model for TB transmission is proposed by Vynnycky and Fine [80].They include both endogenous and exogenous reinfection in the model, plus the effects of BCG vaccination.There is no contact structure for the general population in this model.The infection and reinfection rates are modeled as explicit time-dependent linear functions (piecewise linear function).Numerically, the Euler method is applied to solve the corresponding linear non-autonomous PDE model.Parameters are estimated using data from across the world.Some parameters are estimated from the best fit required by the use of least squares methods.Their target populations are those of England and Wales, where HIV is not a significant factor.Hence, HIV is excluded from the model.Results from their model showed that ARI for age groups 0-15 years, 16-19 years and 20 years above are 4%, 9% and 14%, respectively.
To search for an optimal strategy for TB vaccination, another dynamical model including age structure, contact structure, and vaccination is proposed by Castillo-Chavez and Feng [16].The model includes five subpopulations: susceptible (s), vaccinated (v), latent (l), infectious (i), and treated (j).Their model is given by the following system of PDEs: Parameter Λ is the birth rate; µ(a) is the age specific per-capita natural death rate; β(a) is the age-specific (average) probability of becoming infected through contacts with infectious individuals; c(a) is the age-specific per-capita contact/activity rate; σ and δ represent the reduction in risk of infection due to treatment and vaccination, respectively (0 ≤ σ, δ ≤ 1); k is the per-capita rate of progression to active TB; r is the per-capita treatment rate; and ψ(a) is the per-capita vaccination rate.In addition, p(t, a, a ′ ) gives the "probability" that an individual of age a has a contact with an individual of age a ′ given that the individual of age a had a contact with a member of the population.The infection force is u)du describes the contact structure corresponding to proportionate mixing [13]; and n(t, a) = s(t, a) + v(t, a) + l(t, a) + i(t, a) + j(t, a) is the total population density.
The basic reproductive number R 0 can be calculated following the next generation operator approach [26].Let S(ξ) denote demographic steady state of the population in the absence of disease; A(τ, ξ, η) denote the expected infectivity of an individual who was infected τ units of time ago while at age η; (that is, A(τ, ξ, η) denotes the average infectivity that can be exercised on an uninfected individual at age ξ provided the uninfected population finds itself at the steady demographic state S(ξ))).Under the special assumption of proportionate-mixing [13]; that is, The demographic steady state (the state where the infection is absent) of the system is given by the following nonuniform age-distribution: The probability that an individual of age α + τ , who was infected τ units of time ago, is still in class i, is given by γ(τ, α) = τ 0 ke −(µ+k)u e −(r+µ)(τ −u) du =: K(τ )e −µτ where Thus the expected infectivity is given by 89) leads to the computation of the reproductive number associated with the vaccination strategy ψ.Namely, and K(τ ) is given by (90).Hence, in the absence of a vaccine; that is, whenever ψ(a) = 0 the above formula reduces to Whenever the basic reproductive number R(ψ) < 1, the disease-free steady state is globally stable provided that the total susceptible population has reached the demographic steady state n(a).
Two optimization problems are proposed: if the goal is to bring R 0 (ψ) to less than a pre-assigned value then find the vaccination strategy ψ(a) that minimizes the total cost associated with this goal (reduced prevalence to a target level); if the budget is a fixed cost then find the vaccination strategy ψ(a) that minimizes R 0 (ψ); that is, that minimizes the prevalence.They found two possible optimal strategies, one-age strategy vaccinating the susceptible population at exactly age A; and two-age strategy vaccinating part of the susceptible population at exactly age A 1 and the remaining susceptible population at a later age A 2 .The optimal strategies (one or two ages) also depends on the data, particularly on the type of cost function used.8.2.Treatment of latent TB.Earlier models (prior to the 1970s) targeted the evaluation of control strategies of TB [10,63,81,82], such as vaccination strategies.However, these "optimal" strategies have not worked well toward the elimination of TB globally or even regionally.The reasons behind the lack of success of these policies are debatable.Either these strategies have not been applied by the policymakers or they are not truly "optimal".For instance, although 12% of GNP of the U.S.A. is spent on health care [5], the amount spent on prevention is very small.WHO estimates that if the amount of aid spent on TB treatment programs could be increased from 15 to 100 million dollars yearly, then 1.2 million deaths could be avoided every year [65].That is, the death toll would be reduced by over 30%.This leaves the far-reaching question of what are the best strategies for the complete elimination of TB.Part of this work shows that the focus should include control measures in the latent-TB class.The reason is simple.The huge pool of latent-TB patients is a timebomb or reservoir of infection [62].Forces or new diseases that compromise the immune system may lead to new TB outbreaks, as it has occurred with HIV.Mathematical models that stress the importance of treating individuals with latent TB are introduced by Blower [88] and Castillo-Chavez and Song [73].Adding an early latent class and long-term latent class into the model (3-5) generates the following system: Early latent-TB individuals progress to active TB at the rate pω and to long-term latent TB at the rate (1 − p)ω; long-term latent-TB individuals develop active TB at the rate k; and treatment rates are r 0 , r 1 , and r 2 for early latent TB, long-term latent TB, and active TB, respectively.This two-stage latent TB model distinguishes the risk of progression to active TB for primary infections and longterm infections, which is consistent with the fact that p = 0.05 >> k = 0.00256.(One has seen this in Ferebee's assumptions in section 3.) Model (91-94) allows one to explore the role of treating early latent-TB cases.Results from this two-stage latent TB model show that treatment of 25% of early latent-TB cases together with treatment of 80% of active-TB cases may result in the elimination of TB.
where A(t) has the form: Here α i is constant to be determined via simulations.N (t), now assumed to be independent of the disease, is a known "external" input to the epidemiological model.The values of N (t) are in fact input from extrapolated published U.S.A. demographic data.The transmission rate β is assumed to be a constant; k, TB's activation rate, is also assumed to be constant; r i (i = 1, 2, 3), the treatment rates defined before, are also assumed to be constant; p is the rate at which primary latent-TB cases become permanent latent-TB cases; and µ(t) and d(t) are functions of time.Estimates for some of these parameters are listed in Table 5. Fig. 9 shows ).Fig. 12 illustrates the effect of HIV on the control of TB.It is clear that HIV delays the achievement of CDC's goal but has no permanent impact on the long-term persistence of TB.However, prolonging TB's "survival" enhances the likelihood of its evolution, a situation that is not explored here.We have introduced the impact of HIV/AIDS on TB progression during the past two decades via a temporary perturbation on the distribution of TB progression times.Our selection of this perturbation is driven by our desire to fit active-TB data since our primary goal is to look at not the coevolution of co-infections but at HIV/AIDS co-infections on the ability of the U.S.A. to meet CDC's target by 2010.Our model suggests that if emphasis is placed on treating at least 20% of the latently-infected individuals then CDC's target may be met by 2020.Our model also shows that reemergence of diseases that compromise the immune system (or recurrent outbreaks) would make it very difficult to control TB unless treatment emphasis is put on the earlier (non-detectable) stages of TB disease.
Theoretically sufficient conditions for TB extinction and persistence are derived in terms of upper limits and lower limits of the mortality functions for the nonautonomous model.We place this mathematical result in subsection 8.5 8.4.Regression approach.To partially back up our conclusions from a statistical viewpoint, we use regression to study the trend of new TB cases each year.We let the number of new cases be the response variable, denoted by Y , and time be the predictor, denoted by X.As can be seen from Fig. 1, the scatter plot of annual new cases Y versus year X appears to be exponential.Intuitively, a logarithmic transformation is taken on the response variable.The quadratic regression is the best fit.Fig. 13 shows the fitted curve, 90% confidence bands, and 90% prediction bands.The quadratic regression model turns out to be appropriate after the regression assumptions are verified.Consequently, we can use Equation (98) to predict the number of new cases in the near future.The results shows that the predicted case rate in 2010 is 7.1659/1, 000, 000.The same rate from our deterministic model (95-97) is 5.2952/1, 000, 000.Both figures draw the same conclusion; that is, CDC's goals for TB elimination are unrealistic within the proposeed time horizon.8.5.Asymptotic behavior of the non-autonomous model.An asymptotic analysis of model (95-97) is carried out and the results are discussed [72] in this section since they play a role in the parameterization of the model.The analysis helps establish a criterion for disease persistence (that is, a threshold condition) which must be met by the parameterized model.The long-term behavior of our system is determined by the asymptotic property of the functions N (t), µ(t), and d(t).The following theorem characterizes such behavior.Proof.We will use the following equalities, which are straightforward in real analysis.Whenever lim A exists, the following limit equalities hold: The proof is based on the fluctuation lemma [43] and its extension by Thieme (see Theorem 2.3 in [77]).Applying Theorem 2.3 from [77] to Equations (95) and (97) directly, one obtains that 0 and I(t) is bounded, it follows that I ∞ = 0.A similar argument results in L ∞ = 0.The first part of the theorem is proved.
The theorem provides conditions for differentiation of the two important biological states: disease elimination or persistence for this non-autonomous system.These thresholds are helpful not only in verifying the reasonableness of published parameters but also in the selection of reasonable ranges of unknown parameters.Our model generalizes the results established for related autonomous systems by Feng et al. [31] and Song et al. [74].9. TB transmitted by public transportation.In Argentina, the TB incidence rate in the 1990s was 42/100,000 of the population, but this value was misleading since in the inner city of Buenos Aires it was 160/100,000 (four times higher than the national average).Buenos Aires has 12 million people, and 9.2 million passengers are carried by the bus system, accounting for 82% of the movement in public transportation [19].A specific model targeting the population (in Buenos Aires) that incorporates the impact of public transportation (buses) has been developed by Castillo-Chavez et al. [19].
The city is divided into N neighborhoods.Each neighborhood is further subdivided according to whether or not an individual frequently takes a bus.Type I individuals are those who seldom take buses or do not take them at all while type II individuals are those who frequently take buses.An individual of any type falls into one of four epidemiological groups at any time susceptible (S), infected but not infectious (E), infectious (I), and treated (T ).
Type I and type II individuals have different levels of activity, which are modeled by the contact rates C 1 i and C 2 i , where i indexes the neighborhood.The mixing structure of the population is driven by the bus system, which depends on the average time spent on the bus by the type I members of each neighborhood.To describe the model, let Q j i = S j i + E j i + I j i + T j i be the total number of individuals of type-j in the ith neighborhood, j = 1, 2. The following parameters are required: a i = per-capita average contact rate of type I individuals in neighborhood i; b i = per-capita average contact rate of type II individuals in neighborhood i; σ i = per-capita rate of getting off the bus by type II individuals in neighborhood i; ρ i = per-capita rate of boarding a bus by type II individuals in neighborhood i.
Then, 1 ρi = average time on a bus by a type II person; ρi ρi+σi = probability of staying in the bus (type II person); σi ρi+σi = probability of staying off the bus (type II person).Proportionate mixing is assumed [13].The mixing probabilities are calculated using the above definitions.
is the mixing probability of type I individuals within the same neighborhood i; is the mixing probability between type I and type II individuals within the same neighborhood i; is the mixing probability between type II and type I individuals within the same neighborhood i; ρi+σi is the mixing probability between type II individuals in the ith neighborhood; the mixing probability between type II individuals in the ith neighborhood and type II individuals in the jth neighborhood.
These mixing probabilities satisfy P 11 i + P 12 i = 1 and P 21 i + P 22 i + N j=1 The recruitment rate Λ j i vary across neighborhoods and types.The natural mortality rate µ is assumed to be identical for all neighborhoods and all epidemiological groups.Treatment rates r i , progression rate k i , and the loss of immunity rate (treated person becomes susceptible again) α i depend on the neighborhood but not on the types.The model equations are Q j i = S j i + E j i + I j i + T j i , where i = 1, 2, 3, . . ., N , and j = 1, 2. The superscripts refer to the type and the subscripts index neighborhoods.The equations of type I and type II populations and the population across neighborhoods are coupled by the incidence rates B 1 i (t) and B 2 i (t), where This research found that the larger the difference of prevalence between neighborhoods, the larger the basic reproductive number.After estimating the relevant parameters, it was found that on average 100 people enter and leave the bus hourly, and that one TB infection per 1,000 travelers was generated per hour of travel.Using another model they found bus travel could be responsible for about 30% of new cases of TB [9,19].They also found that variations in TB transmission were most sensitive to transmission within the transportation system.10.Questions and conclusions.
10.1.An old prediction.In the context of TB control, the importance of R 0 was established in 1937, more than two decades before the introduction of first dynamical model for TB [36].W. H. Frost, an epidemiologist at John Hopkins University, addressed the fundamental role of the reproductive number in the following way: "However, for the eventual eradication of tuberculosis it is not necessary that transmission be immediately and completely prevented.It is necessary only that the rate of transmission be held permanently below the level at which a given number of infection spreading (i.e., open) cases succeed in establishing an equivalent number to carry on the succession.If in successive periods of time, the number of infectious hosts is continuously reduced, the end-result of this diminishing ratio, if continued long enough, must be extermination of tubercle bacillus.".... "This means that under present conditions of human resistance and environment the tubercle bacillus is losing ground, and that the eventual eradication of tuberculosis requires only that the present balance against it be maintained."[36] TB is a slow disease, and therefore, the basic reproductive numberR 0 plays a fundamental role in the study its dynamics and control.Its role (transcritical bifurcation) goes well with the observed downward trend of TB mortality and incidence rates.In fact, it suggests that R 0 is "moving" downward as parameters (naturally) change.Mathematically, the results are not surprising, because the models used are modifications of the frameworks developed by Kermack and McKendrick [47] and Ross [66].10.2.Challenging questions.We have collected a number of dynamical models, the results, and insights that have generated in the study of TB dynamics.The historical evolution of dynamical models of TB follows a common pattern in biology from linear to nonlinear, from one strain to multiple strains, from homogenous to heterogeneous, from deterministic to stochastic, from empirical to theoretical (and vice versa).The models are given by system of difference equations, differential equations (ODEs and PDEs), integro-differential equations, and Markov chains.Although we have seen a remarkable progress in the development of a theoretical framework for the study of the dynamics of TB and other epidemiological processes, there many interesting and challenging topics and questions remain.A partial list includes i. Immigration All models "essentially" assume closed populations, ignoring the effects of immigration.Immigration is probably the critical factor in the generation of new TB cases.In the U.S.A. alone, over 40% of total new cases have been among immigrants in the past few years.

ii. Race and ethnicity
There is evidence showing that case rates of TB are different among different groups of people.These changes may be related to variations in susceptibility to the tubercle bacilli.We have seen a complex Markov chain model that takes this into account.However, more work is required if we are to better understand the role of race and ethnicity in disease dynamics.iii.Genetics The major reduction on TB mortality rates was achieved long before the introduction of antibiotics.Can this reduction be explained, at least partially, by the evolution of human susceptibility?Recent work by Aparicio et al. [4] provide a good start.The work on HIV and genetics by Hsu Schmitz [69,70] suggests valuable approaches.iv.Sanitarium The sanitarium waxed and waned historically.It played a critical role in isolating and curing active-TB cases when antibiotics were not available.Models that incorporate the role of isolation on TB control are rare.Frost concluded that it could delay the number of cases [36].v. Global dynamics Theoretically, we characterized the global dynamics of a few models.For most models, the characterization of their global dynamics remains an open question.Multiple strain models, like (10-14) and models with fast and slow progression, like (3-5) should be further analyzed.vi.Time dependence The case of time-dependent coefficients is not only more realistic but often necessary [73].Time-dependent parameters lead to the study of non-autonomous models.To study the long-term behavior of models with time-dependent coefficients, threshold values (like the basic reproductive number) need to be further developed.Obviously, the classic approach for computing the basic reproductive number is not helpful [26].We found upper and lower limits for the persistence and eradication of TB, but it seems hard to get the sharp threshold explicitly.The methods of averages by Ma et al. [55] may be helpful in the study of this problem.Time-dependent models provide a useful way of connecting parameters to data.The work of Aparicio et al. [4] shows this conclusively.vii.Mean latent period The distribution of the latent period is unknown as well as its mean.Knowledge of the shape of this distribution seems critical for control.The results of Feng et al. [31] have shown that it may not have an important qualitative role, but it certainly plays a critical quantitative role.

Figure 1 .
Figure 1.Annual new cases of TB in the United States from 1953 to 2000.Data taken from [22].

Figure 5 .
Figure 5. Backward bifurcation diagram when exogenous reinfection is included in model(20)(21)(22)(23).When R 0 < R p , the diseasefree equilibrium is globally asymptotically stable.However, when R p < R 0 < 1, there are two endemic equilibria.The upper ones are stable, and the lower ones are unstable.

Figure 6 .
Figure 6.Flow chart of the basic cluster model.

Figure 7 .
Figure 7. Global forward bifurcation diagram for the TB models.

5. 1 .
Linear recruitment rate.Currently, most deaths caused by TB represent a small proportion of the deaths for most populations, in other words, d is often insignificant.Therefore, a linear recruitment rate B(N ) = rN with reasonable r values is likely to support exponential growth on a TB-infected population.The use of a logistic recruitment rate B(N ) = rN (1 − N/K) to model the demography in general is also likely to result in logistic growth for the total population N in the presence of TB.To simplify our analysis, we further assume that the infected and reinfected proportions are equal β ′ = β.The use of the variables, N , E and I, is enough.Hence, model(51)(52)(53)(54) reduces to: dN dt = B(N ) − µN − dI, (a). (0, 0, 0) is globally asymptotically stable and lim t→∞ I N = u * , lim t→∞ E N = v * when R 2 < 1, (b).(∞, ∞, ∞) is globally asymptotically stable and lim t→∞ I N = u * , lim t→∞

Figure 8 .
Figure 8. Bifurcation diagram for linear recruitment rate.

6 . 2 .
H+T bC+H+T b , the growth of TB is of Logistic type, r T b T b (K−T b ), and all interactions are characterized by the mass action law.Model(70)(71)(72)(73) undergoes local transcritical and Hopf bifurcations.Bistable equilibria are also found in some ranges of parameter space.The selection of µ T b and r T b as bifurcation parameters helps establish the existence of a local transcritical bifurcation.This cell-based model supports the hypothesis that the presence of TB worsens the clinical picture of HIV-infected individuals and that HIV helps activate TB infections.It also strengthens the view that TB treatment for HIV/TB infected individuals is essential.Multiple levels of macrophages.Another cell-based model involving macrophages, T cells, and M. tuberculosis was proposed by Mathematical and Theoretical Biology Institute students in 1996[42].The function modeling the engulfing of bacteria by macrophages is characterized by three distinct stages: M 0 , normal stage (before encountering bacilli), M 1 inactive stage (containing bacilli but inactive), and M 2 active state (activated by the help of T cells).M. tuberculosis is measured by the bacilli density in the blood of a host.The equations are , a) = β(a)c(a)B(t)(s(t, a) + σj(t, a) + δv(t, a))

Figure 9 .Figure 10 .
Figure 9. New cases of TB and data.

Figure 12 .Figure 13 .
Figure 12.Impact of HIV.The lower curve represents no HIV effect; the upper curve represents the case rate when HIV is included; both are the same before 1983.Dots represent real data.

Table 1 .
Symbols and definitions of subpopulations EExposed infected but unable to infect others (latent or carrier)IInfectious active-TB infections, i.e., he/she can infect others T Treated treated (from latent-TB or active-TB infection)

Table 2 .
Symbols and definitions of parameters.

Table 4 .
Summary of the model with logistic recruitment rate.
[72,73]mination of TB in the U.S.A..A comprehensive and executable model that leads to TB elimination in the U.S.A. is proposed by Song[72,73].Collecting public TB and demographic data for the U.S.A. for the past half century, a model of non-autonomous systems of ordinary differential equations is used to fit U.S.A. tuberculosis incidence over the past five decades.It is shown that tuberculosis in the U.S.A. may be controlled to the point of meeting CDC's criterion of one case per million, but only by the year 2020.This goal may be accomplished only if at least 20% of latently-infected individuals are treated.The effect of HIV/AIDS after 1983 is included in the analysis via a variation of our model that incorporates a function that accelerates TB progression over a window of time.It is shown that TB's case rate may be controlled despite increases in the rate of TB progression due to HIV.From the census and projection data, the total population size N (t) is included in the model as an external input.Two latently-infected classes are introduced, primary latent/exposed class (L 1 ) and a permanently latent class (L 2 ).A non-autonomous ODE model with two latent classes and the incorporation of HIV stands dL 1

Table 5 .
Estimated parameters of TB transmission for the U.S.A..However, treatment of more latently-infected individuals, for instance, raising r 1 and r 2 to 20% per year, would help reach CDC's target of 1/1,000,000 in a more reasonable period of time (see Fig.s 10 and 11