Trade-off between BCG vaccination and the ability to detect and treat latent tuberculosis

ARTICLE IN PRESS Journal of Theoretical Biology 261 (2009) 548–560

Contents lists available at ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Trade-off between BCG vaccination and the ability to detect and treat latent tuberculosis David J. Gerberry  Department of Mathematics, Purdue University, West Lafayette, IN 47907-2067, USA

a r t i c l e in f o

a b s t r a c t

Article history: Received 26 February 2009 Received in revised form 16 July 2009 Accepted 27 August 2009 Available online 4 September 2009

 Policies regarding the use of the Bacille Calmette-Guerin (BCG) vaccine for tuberculosis vary greatly throughout the international community. In several countries, consideration of discontinuing universal vaccination programs is currently under way. The arguments against mass vaccination are that the effectiveness of BCG in preventing tuberculosis is uncertain and that BCG vaccination can interfere with the detection and treatment of latent tuberculosis. In this work, we pose a dynamical systems model for the population-level dynamics of tuberculosis in order to study the trade-off which occurs between vaccination and detection/treatment of latent tuberculosis. We assume that latent infection in vaccinated individuals is completely undetectable. For the case of a country with very low levels of tuberculosis, we establish analytic thresholds, via stability analysis and the basic reproductive number, which determine the optimal vaccination policy, given the effectiveness of the vaccine and the detection/treatment rate of latent tuberculosis. The results of this work suggest that it is unlikely that a country detects and treats latent tuberculosis at a high enough rate to justify the discontinuation of mass vaccination from this perspective. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Tuberculosis BCG Vaccination Mathematical modeling Threshold analysis LTBI treatment

1. Introduction In the fight against tuberculosis, the Bacille Calmette-Gue rin (BCG) vaccine is one of the first lines of defense. While BCG is safe, relatively inexpensive and the most widely used vaccine in the world (Araujo et al., 2008), its use has always been surrounded by uncertainty. The main uncertainty stems from a lack of consistency in clinical trials of the vaccine’s efficacy (i.e. the reduction in the incidence of disease among people who have received the vaccine compared to the incidence in unvaccinated people). Over several decades, hundreds of studies of BCG have reported efficacies that range from 0% to 84% (Centers for Disease Control and Prevention, 1996; Hart and Sutherland, 1977). To complicate matters further, evidence exists that rates of protection depend on geographical location, age at vaccination, form of tuberculosis, particular strain of BCG used and other factors. Another dubious aspect of BCG vaccination is its possible interference with the detection of latent tuberculosis. To detect M. tuberculosis infection, the standard technique is the tuberculin skin test (TST) in which a small amount of tuberculin purified protein derivative is injected into the forearm. After a few days the injection site is observed with the presence and size of the

 Tel.: + 1 310 794 8650.

E-mail addresses: [email protected], [email protected]. 0022-5193/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2009.08.029

reaction used to diagnose tuberculosis infection. However, vaccination with BCG can cause positive TST reactions in uninfected individuals and hence the misdiagnosis of tuberculosis. With so much uncertainty surrounding BCG, it is not surprising that policies regarding its use vary throughout the international community. Currently, the World Health Organization (WHO) recommends universal vaccination of infants in countries with a high TB burden and states that countries with a low TB burden ‘‘may choose to limit BCG vaccination’’ to high-risk individuals (World Health Organization, 2004). In most high burden countries, the WHO recommendations are followed to the extent that resources allow. In low burden countries, policies vary from universal vaccination, to vaccination of high-risk individuals only, to no vaccination at all. In a 2005 survey of the European Union, nearly half of the reporting nations were considering changes to their BCG policies (Infuso and Falzon, 2006). The current vaccination policy of the United States as stated by the Centers for Disease Control and Prevention (CDC) (Centers for Disease Control and Prevention, 1996) is: ‘‘The use of BCG vaccine has been limited because (a) its effectiveness in preventing infectious forms of TB has been uncertain and (b) the reactivity to tuberculin that occurs after vaccination interferes with the management of persons who are possibly infected with M. tuberculosis.’’ In this work, our primary research question is the following: How well must a nation be able to detect and treat latent tuberculosis in order to justify the discontinuation of universal

ARTICLE IN PRESS D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

vaccination? As definitive estimates for the efficacy of the vaccine as well as detection and treatment rates for latent tuberculosis do not exist, mathematical modeling provides an effective means to approach the issue. Clearly there is a trade-off between vaccination and detection/treatment of latent TB. If the vaccine has high efficacy and a country has a low detection/treatment rate of latent TB, vaccination will be the best course of action. On the other hand, if a country is capable of detecting/treating latent TB at a high rate and the vaccine has low efficacy, it would be prudent not to vaccinate the population. We aim to find explicit thresholds which can determine the ‘‘best’’ vaccination policy on a countryspecific basis. The structure of our paper is as follows: in Section 2 we formulate a data-based model for the epidemiology of tuberculosis which includes vaccination, treatment and immigration. In Section 3, we characterize the disease-free equilibrium. In Section 4 we establish and examine analytic thresholds that determine when a universal vaccination policy will be beneficial from an epidemiological point of view. In Section 5, we consider a model which represents the most extreme assumptions against vaccination and examine the resulting thresholds.

2. The model 2.1. Definitions and assumptions As with any modeling endeavor, various assumptions about the underlying biology must be made. At this point, we wish to clearly state some assumptions of this work. We define the ‘‘best’’ vaccination policy in only the epidemiological sense. That is, the policy which results in the lowest levels of tuberculosis. We do not consider the opportunity costs or economic costs of vaccination. In an upcoming paper (Gerberry and Milner, 2010), we focus more on the quantitative differences in TB burdens that result from vaccination using data from several countries. We assume that throughout the duration of the vaccine’s efficacy, latent tuberculosis is completely undetectable. This may be an overly negative assumption as it is widely accepted that tuberculin reactivity wanes rapidly after vaccination in infants (Menzies, 2000; Menzies and Vissandjee, 1992) and hence previous BCG vaccination can often be ignored when interpreting tuberculin skin tests (Centers for Disease Control and Prevention, 2000). We do not explicitly assume a particular form of tuberculosis in our model but rather tuberculosis in general. As pulmonary tuberculosis is responsible for the largest public health burden, one may keep it in mind throughout our discussion. In general our philosophy is to make the strictest assumptions against vaccination whenever we are faced with uncertainty. The rationale for doing so is that the prevailing trend is for low burden countries to discontinue mass vaccination programs. Hence, if conditions which contradict the current trends are found, they exist even under the strictest assumptions against vaccination. As the model is somewhat intricate, we describe its formulation in incremental steps beginning with the basic pathology of TB. We then move to formalize the efficacies of the vaccine, treatment of tuberculosis and immigration.

549

(2004). The model incorporates the following major aspects of tuberculosis pathogenesis: (1) fast and slow progression to active disease, (2) exogenous reinfection of latent TB and (3) variable immunity to reinfection after recovery. Individuals enter the population with a recruitment rate of p. The natural mortality rate is denoted by m. A proportion, c, of individuals entering the population are vaccinated at birth and hence enter the vaccinated class, V. The remaining proportion, 1  c, enter the susceptible class, S. The protection of the vaccination wanes over time at a rate of o and results in the movement of individuals from the vaccinated to the susceptible class. We use the frequency-dependent description of disease transmission, with transmission coefficient b. Upon infection with M. tuberculosis, an individual can either progress quickly to active TB (‘‘fast progression’’), or develop a latent infection and progress slowly to active TB (‘‘slow progression’’). We denote the classes of latently infected and unvaccinated, latently infected and vaccinated, and actively infected individuals by E, EV , and I, respectively. The probability of a random individual exhibiting fast progression to active TB is denoted p. Latent infections progress (or reactivate) to active disease at a rate n. In addition, progression from latent to active TB can be induced through a novel external infection known as exogenous reinfection. Again, only a proportion p of such cases will undergo fast progression. We also allow for the possibility that previous latent infection provides protection against such reinfection and thus multiply the transmission term by a factor y1. We assume that active TB disease induces additional death at a rate denoted by mT and clears at a rate g. Once recovered from active tuberculosis (i.e. class R), one is not immune to infection. In fact, some evidence suggests that an individual is made more susceptible to TB infection as a result of previous active disease (Verver et al., 2005). Here, we denote by y2 the factor of susceptibility that results from previous disease (y2 o1 ) protection, y2 4 1 ) increased susceptibility). The total population is denoted by N. A summary of the model, as described to this point, is given in Fig. 1. 2.3. Vaccine efficacies As the model is currently stated and depicted in Fig. 1, vaccinated classes are unnecessary because the protective effects

2.2. Statement of the model Our mathematical model of tuberculosis is composed of a system of seven ordinary differential equations. It is a generalization of the pre-exposure vaccine model proposed by Ziv et al.

Fig. 1. Basic outline of compartmental model for tuberculosis. Natural death rate, m, omitted for clarity of presentation. Effects of vaccination, treatment and immigration not yet depicted.

ARTICLE IN PRESS 550

D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

of the vaccine are not included. To describe these effects, we assume that the vaccine can protect against tuberculosis via one or more methods: (1) by preventing initial infection, (2) by reducing the probability of fast progression to active TB and (3) by reducing the rate of slow progression from latent to active TB. We denote the probability that the vaccine will prevent an initial infection as q1 . The probability of a newly infected individual being protected from fast progression is denoted by q2. The proportion reduction of the slow progression rate is denoted by q3. When thinking of vaccination, the protection described by q1 is usually the first which comes to mind. However, there is evidence that BCG protects against progression to active TB ðq2 ; q3 Þ as well (Smith et al., 2000), and that this may even be the principal mode of protection (Gomes et al., 2004). It is widely accepted that BCG confers no protection against exogenous reinfection (Smith et al., 2000). The efficacies of vaccination are depicted in Fig. 2.

2.4. Treatment The focus of this work is to establish analytic criteria for when vaccination is beneficial to a society from an epidemiological point of view. As discussed, this will depend heavily on the society’s ability to treat tuberculosis. For this reason, we incorporate treatment into the model. Those individuals with active tuberculosis disease, I, are now divided into those that are treated, IT , and those that are not, IT c . We assume that a new active TB infection is detected with probability d. Also, existing active infections are detected at a rate d=2. The probability of successful treatment is denoted s. We let g1 be the natural clearance rate for active TB, g2 be the rate at which individuals complete treatment and assume that disease-induced mortality only affects untreated individuals. The parameters d and s are defined in this manner to take advantage of the extensive WHO Global TB Database (World Health Organization, 2008) which tracks these indicators. To formalize the trade-off between vaccination and detection of latent TB, we make the pessimistic assumption that latent infection in vaccinated individuals is completely undetectable and hence untreatable. We define r as the probability that a case of latent TB in an unvaccinated individual will be detected and

Fig. 2. Outline of compartmental model for tuberculosis with effects of vaccination included.

successfully treated before the individual leaves the latent class (i.e. before progression to active TB or death). Splitting detection and treatment into separate components is not advantageous as real data for these parameters is not readily available. We choose this characterization of treatment as it is consistent with the manner in which detection/treatment of active TB is estimated. While reliable estimates for this proportion do not exist, we believe that this consistency makes our quantification of treatment more natural to the public health community than a simple rate. However, to implement treatment mathematically we utilize a simple rate, denoted t, that directs latently infected individuals to the recovered class. To reconcile the two terms, we take advantage of the relationship between compartmental ODE models and exponential waiting times. Hence, r, defined as the probability of being treated before leaving the latent class, is given by the rate of treatment divided by the total rate of leaving the latent class. If we neglect exogenous reinfection (a reasonable assumption in a low prevalence society), we get that r ðn þ mÞ. Thus, the intuitive r ¼ t=ðt þ m þ nÞ and hence t ¼ 1r proportion r becomes a model parameter rather than rate t. The incorporation of treatment is shown in Fig. 3. 2.5. Immigration The last necessary addition to our model is immigration. In many North American and Western European countries, which are least likely to have a mass vaccination policy, immigration is a key contributor to tuberculosis dynamics. We define a as the proportion of the recruitment rate which is due to immigration. We let f be the percentage of those immigrants that are vaccinated and c be the percentage of which that are latently infected with M. tuberculosis. We assume that actively infected and recovered individuals do not immigrate. As many countries with large immigration rates screen incoming individuals for infectious disease, it is reasonable to assume that very few individuals with active TB immigrate. Theoretically, recovered individuals can immigrate with no TB-related difficulty. However, reliable estimates of the quantity of such immigration do not exist. As recovered individuals are quite similar to latently infected individuals in our model, neglecting immigration of recovered individuals will have a minimal impact. Fig. 4 depicts

Fig. 3. Compartmental model for tuberculosis incorporating treatment. Here, I ¼ IT þ IT c and N ¼ S þ V þ E þ EV þ IT þ IT c þ R. We denote the total inflow of new active tuberculosis infections as O (i.e. O ¼ nE þ y1 pbEI=N þ pbSI=N þ r ðm þ nÞ pð1  q2 Þð1  q1 ÞbVI=N þ ð1  q3 ÞnEV þ y1 pbEV I=N þ py2 bRI=N). Also, t ¼ 1r where r is the proportion of latent infections which are detected and successfully treated before progression to active TB or death.

ARTICLE IN PRESS D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

551

3. Disease-free equilibrium To establish analytic thresholds for when vaccination is a society’s prudent choice, we consider the case of a population at the brink of eradicating tuberculosis. Mathematically, we must therefore consider the disease-free equilibrium (DFE). We recall that the full system is given by S0 ¼ ð1  aÞð1  cÞp þ að1  fÞð1  cÞp þ oV  bSI=N  mS; V 0 ¼ cð1  aÞp þ apð1  cÞf  oV  ð1  q1 ÞbVI=N  mV; E0 ¼ ð1  pÞbSI=N þ apð1  fÞc  ðnE þ y1 pbEI=NÞ  tE  mE;

Fig. 4. Compartmental model of tuberculosis including immigration.

the changes to the model which result from incorporating immigration. With all of the above considerations included, we arrive at our system of ordinary differential equations for the epidemiology of tuberculosis at the population level. The system is given by

ðð1  q3 ÞnEv þ y1 pbEv I=NÞ  mEv ; IT0

¼ d½nEþ y1 pbEI=N þ pbSI=N þ pð1  q2 Þð1  q1 ÞbVI=N þð1  q3 ÞnEv þ y1 pbEv I=N þ py2 bRI=N  g1 IT  mIT þ d=2IT c ;

IT0 c ¼ ð1  dÞ½nE þ y1 pbEI=N þ pbSI=N þ pð1  q2 Þð1  q1 ÞbVI=N þð1  q3 ÞnEv þ y1 pbEv I=N þ py2 bRI=N þ ð1  sÞg1 IT  g2 IT c ðm þ mT ÞIT c  d=2IT c ;

S0 ¼ ð1  aÞð1  cÞp þ að1  fÞð1  cÞp þ oV  bSI=N  mS;

V 0 ¼ cð1  aÞp þ apð1  cÞf  oV  ð1  q1 ÞbVI=N  mV;

R0 ¼ sg1 IT þ g2 IT c þ tE  y2 bRI=N  mR:

ð2Þ

In the transmission terms, we denote the total number of infectious individuals in the population as I ¼ IT þ IT c and the total population as N ¼ S þ V þ E þEV þ IT þ IT c þ R for notational convenience. The system does not possess a DFE due to the inflow of infection resulting from immigration. To study the case where eradicating tuberculosis is a possibility, we must restrict our attention to a closed population. Assuming no immigration, the resulting model is given by

E0 ¼ ð1  pÞbSI=N þ apð1  fÞc  ðnE þ y1 pbEI=NÞ  tE  mE;

Ev0 ¼ ð1  pð1  q2 ÞÞð1  q1 ÞbVI=N þ apcf þ ð1  pÞy2 bRI=N ðð1  q3 ÞnEv þ y1 pbEv I=NÞ  mEv ; IT0 ¼ d½nE þ y1 pbEI=N þpbSI=N þ pð1  q2 Þð1  q1 ÞbVI=N þ ð1  q3 ÞnEv þ y1 pbEv I=N þpy2 bRI=N  g1 IT  mIT þ d=2IT c ;

S0 ¼ ð1  cÞp þ oV  bSI=N  mS;

IT0 c ¼ ð1  dÞ½nEþ y1 pbEI=N þ pbSI=N þpð1  q2 Þð1  q1 ÞbVI=N þ ð1  q3 ÞnEv þ y1 pbEv I=N þpy2 bRI=N þ ð1  sÞg1 IT  g2 IT c ðm þ mT ÞIT c  d=2IT c ; R0 ¼ sg1 IT þ g2 IT c þ tE  y2 bRI=N  mR;

Ev0 ¼ ð1  pð1  q2 ÞÞð1  q1 ÞbVI=N þ apcf þð1  pÞy2 bRI=N

V 0 ¼ cp  oV  ð1  q1 ÞbVI=N  mV; E0 ¼ ð1  pÞbSI=N  ðnEþ y1 pbEI=NÞ  tE  mE; ð1Þ

where I ¼ IT þ IT c and N ¼ S þ V þ E þEV þ IT þ IT c þR. Recent works, conducted independently of that presented here, seem to also use (Ziv et al., 2004) as a basis for examining vaccination and treatment of tuberculosis (Bhunu et al., 2008a, b). However, these studies differ significantly from the current work in that they do not address the trade-off which exists between vaccination and detection/treatment of latent TB. Instead, the authors assume no treatment of latent TB while modeling BCG vaccination and then consider new vaccines which are currently in development that will allow for detection of latent TB. Additionally, all active infections are grouped into one class as opposed to our approach which allows the direct incorporation of WHO-estimated disease indicators d and s. A few simulations are helpful in order to obtain a feeling for the behaviors that the model yields. Here we use biologically reasonable parameters for the epidemiology of tuberculosis in the United States. The derivations of these parameters are discussed in more detail in Gerberry and Milner (2010) and less so in Section 4.3 to follow. In Fig. 5, we show the effect that varying the levels of vaccination and treatment has on the country’s latent TB dynamics. The simulations indicate the sensitivity of the system in regards to treatment and seem to suggest an insensitivity to vaccine coverage. However, this could be linked to the relatively short time period of the simulation (20 years).

Ev0 ¼ ð1  pð1  q2 ÞÞð1  q1 ÞbVI=N þ ð1  pÞy2 bRI=N ðð1  q3 ÞnEv þ y1 pbEv I=NÞ  mEv ; IT0 ¼ d½nEþ y1 pbEI=N þ pbSI=N þ pð1  q2 Þð1  q1 ÞbVI=N þð1  q3 ÞnEv þ y1 pbEv I=N þ py2 bRI=N  g1 IT  mIT þ d=2IT c ; IT0 c ¼ ð1  dÞ½nE þ y1 pbEI=N þ pbSI=N þ pð1  q2 Þð1  q1 ÞbVI=N þð1  q3 ÞnEv þ y1 pbEv I=N þ py2 bRI=N þ ð1  sÞg1 IT  g2 IT c ðm þ mT ÞIT c  d=2IT c ; R0 ¼ sg1 IT þ g2 IT c þ tE  y2 bRI=N  mR: If a population is free of tuberculosis E ¼ Ev ¼ IT ¼ IT c ¼ 0), system (3) reduces to

ð3Þ infection

(i.e.

S0 ¼ ð1  cÞp þ oV  mS; V 0 ¼ cp  oV  mV; R0 ¼  mR: Recalling that S0 ¼ V 0 ¼ R0 ¼ 0 is necessary for an equilibrium, we see that the DFE is given by ðS ; V  ; E ; EV ; IT ; IT c ; R Þ ¼

  ð1  cÞmp þ op cp ; ; 0; 0; 0; 0; 0 : mðm þ oÞ mþo

ð4Þ

ARTICLE IN PRESS 552

D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

US wIMMIG wDEMOG Compare d q1 = 0.24 q2 = 0.24 q3 = 0.24 s = 0.83 r = 0.2 c = 0 ω = 0.0667

proportion of population with Latent TB

proportion of population with Latent TB

US wIMMIG wDEMOG Compare c q1 = 0.24 q2 = 0.24 q3 = 0.24 d = 0.92 s = 0.83 r = 0.2 ω = 0.0667

0.086 0.084

c=0 c = 0.33 c = 0.66 c=1

0.082 0.08 0.078 0.076 0.074 0.072 0.07 0

2

4

6

8

0.086 0.084

d = 0.65 d = 0.75 d = 0.85 d = 0.95

0.082 0.08 0.078 0.076 0.074 0.072 0.07

10 12 14 16 18 20

0

2

4

6

US wIMMIG wDEMOG Compare s q1 = 0.24 q2 = 0.24 q3 = 0.24 d = 0.92 r = 0.2 c = 0 ω = 0.0667

0.086

s s s s

0.084 0.082

= 0.65 = 0.75 = 0.85 = 0.95

0.08 0.078 0.076 0.074 0.072 0.07 0

2

4

6

8

8

10 12 14 16 18 20

time (years)

proportion of population with Latent TB

proportion of population with Latent TB

time (years)

US wIMMIG wDEMOG Compare r q1 = 0.24 q2 = 0.24 q3 = 0.24 d = 0.92 s = 0.83 c = 0 ω = 0.0667

0.086

r=0 r = 0.1 r = 0.2 r = 0.3

0.084 0.082 0.08 0.078 0.076 0.074 0.072 0.07

10 12 14 16 18 20

0

2

4

time (years)

6

8

10 12 14 16 18 20

time (years)

Fig. 5. Simulations for TB in the United States showing the percentage of the population that is latently infected. Each graph depicts four simulations which assume different levels of treatment or prevention while fixing all other parameters. Simulations use b ¼ 11, y1 ¼ 0:60 and y2 ¼ 1:0. The remaining parameter values are specified in either the figure itself or Table 1. (a) Varying vaccination coverage. (b) Varying active TB detection rate. (c) Varying treatment success rate. (d) Varying detection/treatment rate for latent TB.

4. Thresholds for BCG benefit

R0 ¼ gI þ tE  y2 bRI=N  mR:

4.1. Model with no immigration and no treatment of active TB

We now denote the average length of an active infection as 1=g.

We are considering a population in which tuberculosis is nearly extinct. Thus, it is reasonable to study a simplified version of the model which does not include treatment of active TB. While we do consider the more realistic case of including treatment for active TB in Section 4.2, we make this temporary assumption in order to facilitate a stability analysis of the DFE. With this simplification, our model has only one class of infectious individuals and is given by

4.1.1. Stability analysis of DFE Let X ¼ ðS; V; E; EV ; I; RÞT , and rewrite (5) as X 0 ¼ FðXÞ. The DFE equilibrium of (5) is given by  T ð1  cÞmp þ op cp ; ; 0; 0; 0; 0 : X  ¼ ðS ; V  ; E ; Ev ; I ; R ÞT ¼ mðm þ oÞ mþo ð6Þ

S0 ¼ ð1  cÞp þ oV  bSI=N  mS; V 0 ¼ cp  oV  ð1  q1 ÞbVI=N  mV; E0 ¼ ð1  pÞbSI=N  nE  y1 pbEI=N  tE  mE;

Linearizing our system about the DFE, we find that the Jacobian of F about X  is 2 3 m o 0 0 bS 0 6 7 0 0 bV  ð1  q1 Þ 0 7 6 0 m  o 6 6 0 6 DFðX Þ ¼ 6 6 0 6 6 0 4 0 

n  m

0

0

0

m  nð1  q3 Þ

B

0 0

ð1  rÞn rn

ð1  q3 Þn 0

C

0

I0 ¼ nE þ y1 pbEI=N þ pbSI=N þpð1  q2 Þð1  q1 ÞbVI=N þ ð1  q3 ÞnEv þ y1 pbEv I=N  gI  ðm þ mT ÞI þ py2 bRI=N;

A

g

7 0 7 7 7; 0 7 7 0 7 5 m

ð7Þ

Ev0 ¼ ð1  pð1  q2 ÞÞð1  q1 ÞbVI=N  ð1  q3 ÞnEv  y1 pbEv I=N mEv þð1  pÞy2 bRI=N;

ð5Þ

where A ¼ ð1  pÞbS ; B ¼ ð1  pð1  q2 ÞÞð1  q1 ÞbV  ;

ARTICLE IN PRESS D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

C ¼ pbS þ pbV  ð1  q1 Þð1  q2 Þ  mT  m  g:

ð8Þ

Noticing the form of the Jacobian, we immediately have that m  o is an eigenvalue and that m is a double eigenvalue. The remaining eigenvalues are those of the 3  3 submatrix 2 3 n  m 0 A 6 0 m  nð1  q3 Þ B 7 ð9Þ 4 5: ð1  q3 Þn C ð1  rÞn To facilitate the analysis, we temporarily assume that q3 ¼ 0. The eigenvalues of the submatrix are then given by

l1 ¼  n  m; 

l2;3 ¼ 12 C  n  m 7

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðC þ n þ mÞ2 þ4nðAð1  rÞ þ BÞ :

Noting that A and B are nonnegative for biologically feasible parameters and that 0 r r r 1, we see that the eigenvalues are real-valued. It also follows that  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi l3 ¼ 12 C  n  m  ðC þ n þ mÞ2 þ 4nðAð1  rÞ þ BÞ  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi o 12 C  n  m  ðC þ n þ mÞ2 h i ¼ 12 C  n  m  jC þ n þ mj : If C Z0, we have that l3 o 12 ½C  n  m  C  n  m ¼  n  m. l3 o 12 ½C  n  m  jC þ n þ mj If C o0, we have that 1 o 2 ½C  n  m o 0. Hence, the stability of the DFE will be determined by l2. To study the effect of vaccination on the stability of the DFE, we begin by setting l2 ¼ 0 which is equivalent to

nðB  ð1  rÞAÞ : nþm

C¼

ðm þ g þ mT Þmðm þ oÞðn þ mÞ : pðm þ oÞ½ð1  r þ rpÞn þ pmþ ½ðmpq2  pm  nÞq1  mpq2 þ rnð1  pÞmpc ð10Þ

b ¼



We have that for b o b the DFE will be asymptotically stable. Thus, vaccination improves the stability of the DFE whenever it  increases b . Considering (10), we see that vaccination is advantageous in the situation where ðmpq2  pm  nÞq1  mpq2 þ r nð1  pÞ o 0: Solving for r results in ro

pmq1  pmq2 q1 þ q1 n þ mpq2 : nð1  pÞ

ð11Þ

Hence, we have a threshold value for r, which we denote r : r ¼

die out. If R0 4 1, the number of infected individuals rises and an epidemic results. Clearly, the stability of the DFE and R0 are interconnected. Using the reproductive number, we derive a slightly more general threshold value for the treatment of latent tuberculosis. We use the formulation of R0 presented by Diekmann and Heesterbeek (2000). To begin, we note that our model has three disease states: latently infected and unvaccinated ðEÞ, latently infected and vaccinated ðEv Þ and infectious ðIÞ. ~ the probability that a new infection will be of We denote by Y i type i for i ¼ 1; 2; 3, corresponding to E, Ev , I, respectively. To find the probabilities, we recall that the DFE is given by (4). As the total population is p=m, we see that the proportion of the population vaccinated at the DFE is given by pV ¼ cm=ðo þ mÞ and thus the proportion unvaccinated is pS ¼ 1  pV . Note that the proportion of the population vaccinated is not c as we may expect. This is because we assume that the vaccine does not confer lifelong immunity and as its efficacy wanes the individual moves to the susceptible class. The probability that a new infection will occur in an unvaccinated individual is given by the transmission rate to unvaccinated individuals divided by the total rate of transmission. Since vaccination reduces the transmission rate by a factor of q1 , we get that the probability of a new infection occurring in an unvaccinated individual is

bpS =N : bpS =N þ ð1  q1 ÞbpV =N Recalling that a proportion p of new infections will exhibit fast progression to active TB and hence be of type I, we can finalize the term ~1 ¼ Y

Substituting the expressions given by (8) and (6), we find the equivalent equality involving b, the transmission coefficient,

pmq1 ð1  q2 Þ þ q1 n þ mpq2 : nð1  pÞ

ð12Þ

Thus, if a country is able to detect and treat latent tuberculosis infection at a rate which exceeds r  , the discontinuation of BCG vaccination will increase the stability of the DFE. If the detection/ treatment rate of latent TB is below r , BCG vaccination will result in a more stable DFE equilibrium. 4.1.2. Basic reproduction number Along with the stability of the DFE, the basic reproductive number is an important indicator of an epidemic’s growth in a population in which tuberculosis is nearly eradicated. The basic reproductive number, R0 , is defined as the average number of secondary infections caused by a single infectious individual introduced into a susceptible population. If R0 o 1, an average infectious individual is unable to replace itself and the disease will

553

ð1  pÞbpS =N ð1  pÞ½ð1  cÞm þ o ¼ ; bpS =N þ ð1  q1 ÞbpV =N ð1  cq1 Þm þ o

which gives the probability of a new infection being latent and occurring in an unvaccinated individual. In a similar manner, we find that the probability that a new latent infection occurs in a vaccinated individual is ~2 ¼ Y

½1  pð1  q2 Þð1  q1 ÞbpV =N ½1  pð1  q2 Þð1  q1 Þmc ¼ : bpS =N þ ð1  q1 ÞbpV =N ð1  cq1 Þm þ o

Finally, we have that the probability of a new infection being active is ~1  Y ~2 ¼ ~3 ¼ 1  Y Y

p½mð1  cÞ þ o þ ð1  q1 Þð1  q2 Þmc : ð1  cq1 Þm þ o

We next define the transition rate matrix G of transition probability per unit time. In our case, G is a 3  3 matrix because we have three disease states. We let Gij be the probability per unit time of an individual transferring from state j to state i for ia j. We define Gii to be the probability of an individual leaving state i per unit time, whether it be transitioning into another class or death. For our model, we have 2 3 0 0 n  m  t 6 7 0 0 m  ð1  q3 Þn ð13Þ G¼4 5: n ð1  q3 Þn m  mT  g If we consider one infected individual and let ~ x ðtÞ be the probability of that individual being in the various states at time t, it follows that ~ x satisfies d~ x ¼ G~ x; dt ~; ~ x ð0Þ ¼ Y

ð14Þ

~ . Noting the form of G, we which has solution given by ~ x ðtÞ ¼ eGt Y see that the eigenvalues are all negative for biologically reasonable parameters. Thus, xðtÞ-0 as t-1, as we would expect. We

ARTICLE IN PRESS 554

D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

can then calculate the expected amount of time that an individual will spend in the various states by integrating ~ x with respect to time from 0 to 1. Doing so, we arrive at Z 1 1 ~ dt ¼ G1 eGt Y ~  ¼  G1 Y ~: eGt Y ð15Þ 0

0

Once we have the expected amount of time spent in each state, we simply multiply by the rate of disease transmission for each state and sum to get the basic reproductive number. If we let the vector ~ h represent infectivity in each state, this is done with the ~ . In our case, I is the only infectious h  G1 Y dot product R0 ¼  ~ T state so we have that h ¼ ð0; 0; bpS þ ð1  q1 ÞbpV Þ. Putting this all together, we arrive at our main result. The basic reproductive number for our model is given by  ~¼  0 h  G1 Y R0 ¼  ~ 2

bðm  mcq1 þ oÞ

0



mþo

n  m  t

0

0

31

0

m  n þ q3 n

0

n

ð1  q3 Þn

m  mT  g 3

7 5

6 4 2

ð1  pÞ½ð1  cÞm þ o 7 6 ð1  cq1 Þm þ o 7 6 7 6 ½1  pð1  q2 Þð1  q1 Þmc 7 6 7 6 6 7 Þ m þ o ð1  cq 1 7 6 6 p½mð1  cÞ þ o þð1  q Þð1  q Þmc 7 1 2 5 4 ð1  cq1 Þm þ o  nðp  1Þðm þ rnÞq3 ¼ Kðn  nr þ nrp þ pmÞ þK mc ðm þ oÞ½nð1  q3 Þ þ m  ðm þ nÞðmpq2  mp þ nq3  nÞq1 ðm þ nÞ½mpq2 þ r nðp  1Þ  þ ðm þ oÞ½nð1  q3 Þ þ m ðm þ oÞ½nð1  q3 Þ þ m ð16Þ

b ðm þ mT þ gÞðm þ nÞ

:

In (16), we have replaced t by rðm þ nÞ=ð1  rÞ. With the reproductive number established, finding the threshold for when vaccination is the better strategy is straightforward. Vaccination is advantageous when it lowers R0 . Examining (16), we see that this is exactly the case where the coefficient of c is negative. Noticing that the denominator is strictly positive for biologically feasible parameters, we can restrict our attention to the numerator. Setting this to zero and solving for r, we find that the threshold value for detection/treatment of latent tuberculosis that justifies discontinuing BCG vaccination is r ¼

The modifications for G and ~ h are given by 2 3 n  m  t 0 0 0 6 7 0 m  ð1  q3 Þn 0 0 6 7 7 G¼6 6 7 d n dð1  q Þ n  g  m d=2 3 1 4 5 ð1  dÞð1  q3 Þn ð1  sÞg1 d=2  g2  m  mT ð1  dÞn and  T ~ h ¼ 0

0

bðm  mcq1 þ oÞ

bðm  mcq1 þ oÞ

mþo

mþo

 :

~ , we get the basic reproductive Recalling that R0 ¼  ~ h  G1 Y number for the model which includes treatment of active tuberculosis:

where K¼

still restrict attention to the case of a closed population. The system of differential equations is (3). First, we realize that our current system has four disease states. We will associate the indices i ¼ 1; 2; 3; 4 with disease states E, Ev , IT , IT c , respectively. To calculate R0 for the model including treatment of active TB, we must simply adjust our previous ~ , the probability ~ , G, and ~ h. To adjust Y expressions for Y distribution for a new infection, we need to separate the probability of a new active infection into those associated with the treated and untreated classes. In our model, the parameter d denotes the percentage of new active infections which are ~ follows detected and thus treated. The new expression for Y immediately: 3 2 ð1  pÞ½ð1  cÞm þ o 7 6 Þ m þ o ð1  cq 1 7 6 7 6 ½1  pð1  q Þð1  q Þ m c 2 1 7 6 7 6 ð1  cq1 Þm þ o 7 6 ~ 6 Y ¼ 6 dp½mð1  cÞ þ o þð1  q Þð1  q Þmc 7 7: 1 2 7 6 7 6 Þ m þ o ð1  cq 7 6 1 7 6 4 ð1  dÞp½mð1  cÞ þ o þ ð1  q1 Þð1  q2 Þmc 5 ð1  cq1 Þm þ o

ðm þ nÞðpm  mpq2 þ n  nq3 Þq1 þ ½nð1  pÞq3 þ pq2 ðm þ nÞm : nð1  pÞðn  nq3 þ mÞ



nðp  1Þðm þ nrÞq3 ðm þ oÞðn  nq3 þ mÞ  ðn þ mÞðpmq2  pm  n þ nq3 Þq1 ðn þ mÞ½mq2 p þ nrðp  1Þ  ; þ ðm þ oÞðn  nq3 þ mÞ ðm þ oÞðn  nq3 þ mÞ ð18Þ

R0 ¼ Kðpm þ n  nr þ nrpÞ þ K mc

where K¼

ðg1  g1 dsþ m þ dg2 þ d=2 þ dmT Þb : ðn þ mÞ½ðm þ mT þ g2 þ ds=2Þg1 þ mðd=2 þ m þ mT þ g2 Þ

Again, we find the threshold for the benefit of vaccination by setting the coefficient of c to 0. Doing so, we get ð17Þ

In using the stability of the DFE equilibrium to establish this threshold, we were forced to assume q3 ¼ 0 in order to allow for manageable expressions for the eigenvalues. One can verify that the threshold given by (17) is indeed a generalization of (12) by noting that the two agree if we let q3 ¼ 0 in (17). 4.2. Model with treatment of active tuberculosis As we have seen in Section 4.1, calculating the reproductive number allows us to find the threshold for detection/treatment of latent tuberculosis more efficiently than studying the stability of the DFE. Even under the assumption of no treatment for active TB, the stability analysis proved to be restrictive enough to justify assuming that q3 ¼ 0. Here, we generalize our model to include treatment of active TB and use R0 alone to find our threshold. We

r ¼

ðm þ nÞðpm  mpq2 þ n  nq3 Þq1 þ ½nð1  pÞq3 þpq2 ðm þ nÞm : nð1  pÞðn  nq3 þ mÞ ð19Þ

This is precisely the same as (17), the threshold obtained when treatment was not included in the model. This is only slightly surprising as R0 deals with introducing one infection into a completely susceptible population. In such a scenario, the treatment of active TB should be expected to have a minimal effect on whether BCG vaccination is the best strategy. 4.3. Biologically reasonable parameters Before examining the threshold closer, we establish some biologically appropriate parameters for the epidemiology of tuberculosis. For illustration, we will obtain estimates of the

ARTICLE IN PRESS D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

relevant parameters for the United States. In the discussion that follows, we will then be able to obtain a more quantitative feel for the information that is given by our threshold r  . It is commonly estimated that a third of the world’s population is latently infected with tuberculosis, of which 5–10% will develop active TB during their lifetime (World Health Organization, 2004). Estimates of the percentage of newly infected individuals that display fast progression have been given at 1–5% in Kumar et al. (2007) and a range of 0–30% was used in Blower et al. (1995). For simplicity, we will assume that the probability of fast progression is 0.05 (i.e. p ¼ 0:05) and that approximately 5% of individuals progressing slowly to active TB will do so before death. The average life expectancy for the US is now estimated to be 78 years which immediately gives us a value for m. As our deterministic system of equations implies that the amount of time an individual spends in each class is exponentially distributed, it follows that the probability of progressing to active TB before death is given by n=ðn þ mÞ. Equating this to 0.05, we can obtain an estimate of n. Although d and s are not present in the expression for r , they provide insight to reasonable values for r. WHO data has the treatment success rate ðsÞ for the United States at 83% for 2000– 2003 and at 61% for 2004. The all new case detection rate ðdÞ ranges from 103% to 106% from 2000 to 2005 (World Health Organization, 2008). While such rates are obviously due to imprecision of estimation, we will assume values of d in the range 90–100%. As studies for the overall efficacy of BCG have had widely divergent results, it is not surprising that the degree to which BCG vaccination utilizes each of the three methods of prevention ðq1 ; q2 ; q3 Þ is unknown. We include these parameters which are relevant for the threshold r  , as well as biologically feasible values and ranges for all model parameters in Table 1. Again, we emphasize that a detailed discussion of parameter values and the process used to calibrate the model to TB data from eight countries of varying TB burden is discussed in an upcoming paper (Gerberry and Milner, 2010).

555

With an explicit expression for R0 and biologically reasonable parameter estimates stated, we can now examine the corresponding basic reproductive numbers for our initial simulations for the TB dynamics of the United States in Fig. 5. The parameter values used form the simulations in Fig. 5a give reproductive numbers of 0.504, 0.494, 0.485 and 0.476 for c ¼ 0; 0:33; 0:66 and 1.0, respectively. For Fig. 5b, we have R0 ¼ 1:09; 0:84; 0:63 and 0.45 for d ¼ 0:65; 0:75; 0:85 and 0.95, respectively. In Fig. 5c, we have R0 ¼ 0:93; 0:68; 0:46 and 0.28 for s ¼ 0:65; 0:75; 0:85 and 0.95, respectively. For the simulations of Fig. 5d, reproductive numbers of 0.56, 0.53, 0.50 and 0.48 result from r ¼ 0; 0:1; 0:2 and 0.3, respectively. As with the simulations, we see that the basic reproductive number is most sensitive to changes in the parameters for detection and treatment of active disease, d and s. It is also interesting to note that many values for R0 are less than 1, which mathematically would lead to the eradication of TB in the United States. However, it should be emphasized that immigration had to be neglected in order to study the basic reproductive number, despite the role it plays as a major source of infection in the United States.

4.4. Examining the threshold r  The threshold established in the previous section gives us an explicit criterion for the level of detection/treatment of latent TB that is necessary in order for the discontinuation of vaccination to increase the stability of the DFE. To emphasize, the parameter r is defined as the probability that a latent infection is detected and successfully treated. This is in contrast to the way that treatment of active infection is described in the model. In this case, d gives the probability that an active case of TB is detected and s gives the probability that a case is treated successfully. The discrepancy is due to the fact that the WHO explicitly reports on the treatment and detection of active TB separately, while the corresponding indicators for latent tuberculosis are not reported. Clearly,

Table 1 Biologically reasonable parameter values and ranges for the United States. Parameter

Symbol Range of values Min

Probability of fast progression Rate of slow progression Average duration of treatment for active TB Natural cure rate for active TB Mortality rate due to TB Proportion of immigrants that are vaccinated Proportion of immigrants with LTBI Average life expectancy (years) Overall recruitment rate (millions/year) Proportion of recruitment due to immigration Vaccination coverage Detection rate for active TB Treatment success rate for active TB Transmission coefficient Average duration of vaccine protection Factor of protection due to latent infection Factor of protection/susceptibility due to previous active infection

Mode

p

0.05

n

0.00075 0.0769 0.058 0.12 0.80 0.30 78 6.416 0.22 0.00

1=g1

g2 mT f c 1=m

p a c d s

b 1=o

y1 y2

0.90 0.60 10 10 15 0.50 0.60 0.80 1.00

Source Max

1.0 0.85 12 55 0.75 1.20

Smith and Moss (1994), Blower et al. (1995), Porco and Blower (1998)a Young and Dye (2006), Frieden et al. (2003)b Mayo Foundation for Medical Education and Research (2009)c Blower et al. (1995), Porco and Blower (1998) Blower et al. (1995), Porco and Blower (1998) Assumption Assumption United States Census Bureau (2008) Fitd United States Census Bureau (2008) World Health Organization (2007)e World Health Organization (2008) World Health Organization (2008) Fit Sterne et al. (1998), Aronson et al. (2004) Cohen et al. (2007), Dye et al. (1998) Assumption

a Estimates of the probability of progression to active infection within a few years (i.e. primary disease) range from 5–10% (Smith and Moss, 1994). We assume 5% as our system models fast progression immediately after initial infection. b Young and Dye (2006), and Frieden et al. (2003) estimate a 5–10% lifetime risk of an infected individual developing active TB. o is chosen so that the probability of reactivation before death (i.e. n=ðn þ mÞ) is 0.05, assuming an average life expectancy of 70 years. Given a 5% chance of fast progression and a 5% chance of slow progression, the overall lifetime risk of developing active TB is slightly less than 10%. c TB treatments last 6–12 months, but patient is noninfectious after a few weeks. We assume 4 weeks of infectiousness on average to account for a short period before detection. d Derived by fitting population curve to data from United States Census Bureau (2008). e Used reported coverage from 2006.

ARTICLE IN PRESS 556

D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

1.0

0.8

r

0.6

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

Fig. 7. r  ¼ 0:35 (light), r  ¼ 1:0 (dark).

q Fig. 6. rq1 (dashed), rq2 (solid) and rq3 (dotted) for US TB parameters.

estimating such indicators would be a monumental undertaking as the vast majority of the latently infected population will never even know of their infection. While estimates of treatment rates for latent TB are scarce at best, we would be safe to assume that they are considerably lower than their counterparts for active infection. Mathematically, this translates to saying that we can certainly expect that r o ds. To begin our study of the threshold r  given by (19), we consider the threshold under the assumption that the disease exhibits no latent class. Mathematically, this situation can be achieved by recalling that p gives the proportion of new infections that exhibit fast progression to active TB. Hence, by letting p-1 we eliminate the latent class altogether by having all new infections bypass this stage. A minor computation gives that lim r  ¼ 1:

ð20Þ

p-1

This is as expected because in the absence of the latent stage, our ability to treat latent infection will have absolutely no impact and hence will never justify the discontinuation of vaccination. Another interesting case is that of perfect vaccination. Here, we assume that the vaccine is completely effective at preventing initial infection ðq1 Þ, reducing the probability of fast progression to active TB ðq2 Þ, and reducing the rate of slow progression to active TB ðq3 Þ. Mathematically, this is modeled by setting q1 ¼ q2 ¼ q3 ¼ 1. Doing so, we see that the threshold rate for successfully treated latent infections becomes prohibitively large: r  jðq1 ;q2 ;q3 Þ ¼ ð1;1;1Þ ¼

n þ mp n þ mp 4 41: nð1  pÞ n

This confirms the intuition that vaccination is always the best policy if the vaccine is completely effective. As was mentioned, the way in which BCG vaccination prevents TB infection is unclear. Some suggest that the vaccine prevents slow progression to active TB more than it prevents the initial infection (Gomes et al., 2004). Here, we consider each mode of prevention separately. This is easily achieved by setting the remaining efficacies to zero. We obtain rq1 ¼

q1 ðn þ mpÞ ; nð1  pÞ

For a vaccine that prevents initial infection only (i.e. q1 only), we see that an efficacy of about 50% is sufficient to guarantee that mass vaccination results in a more stable DFE equilibrium, regardless of the rate at which latent TB is detected and successfully treated ðrÞ. For a vaccine that only reduces the probability of fast progression to active TB (i.e. q2 only), the data for the US gives the approximate relation that r   q2 . In the case of a vaccine only preventing slow progression to active TB (i.e. q3 only), we also see that r   q3 . To obtain a feeling for a vaccination which confers all three methods of protection, we examine the surfaces which result if we fix various values of r  . As depicted in Fig. 7, we see that even if we can detect and treat all latent infections, vaccination is still justified for a wide range of parameter space. As any realistic value of r is probably quite small, it follows that vaccination would result in a more stable DFE in most situations.

rq2 ¼

q2 mp ; nð1  pÞ

rq3 ¼

q3 m : m þ nð1  q3 Þ

ð21Þ

For the parameters given in Table 1, we graph the three thresholds in Fig. 6.

5. Worst-case scenario for vaccination In defining the model in Section 2, a few implicit assumptions were made that work in favor of deciding to implement a mass vaccination program. Recalling the dynamical system which represents our model, given by (1) and Figs. 1–4, we describe these issues and pose another model that the authors believe to be the ‘‘worst-case scenario’’ for vaccination. While such a description is almost certainly too pessimistic of vaccination, it will provide the most restrictive thresholds for the levels of detection/treatment of latent tuberculosis ðrÞ which justify a mass vaccination program. Noting the form of the threshold for detection/treatment of latent TB given by (19), we see that r  is independent of the average duration of vaccination protection (i.e. 1=o). To explain this, we note that as vaccination protection wanes an individual becomes susceptible. Thus, a subsequent latent infection will then be considered detectable and benefit from the detection/treatment rate r. As our threshold only weighs the effects of vaccine protection ðq1 ; q2 ; q3 Þ versus treatment of latent TB ðrÞ, a very short duration of protection is not detrimental to the vaccination policy because these individuals can still benefit from treatment of latent TB. To penalize a vaccine with very short duration, we include an additional class of individuals consisting of those who were vaccinated but through the waning of vaccine protection are now fully susceptible. We denote this class VS .

ARTICLE IN PRESS D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

Another aspect of our initial model that is beneficial to the case of vaccination is the lack of vaccine waning for latently infected individuals. As latent infections are typically long-term, it is reasonable to incorporate further waning of vaccine protection during this period. In order to do so, we add a class of latently infected individuals who are undetectable, but do not benefit from the protection of the vaccine, q3 . We denote this class ES as it will also include individuals in VS who become latently infected. Lastly, our initial model dictates that individuals who are reinfected with M. tuberculosis move to EV . The purpose of doing so was that these latent infections would be undetectable due to the previous active infection. However, this also confers the protection associated with vaccination to these individuals. To avoid such a situation, we now divert reinfected individuals into the newly created ES class. It is important to reiterate that these adjustments reflect the worst-case scenario for vaccination and are almost surely overly pessimistic of the effects of vaccination. In most clinical trials, tuberculosis status is tracked for groups of vaccinated and unvaccinated individuals over the course of several years and efficacy is calculated by simply dividing the disease rate of the vaccinated group by that of the unvaccinated group. In such a study, it is difficult to distinguish between waning of protection (o) and pure lack of effectiveness ðq1 ; q2 ; q3 Þ. Thus, a reported BCG efficacy of 80% and a reported duration of protection of 15 years, for example, may not be independent findings. The current model adjustments may double penalize a vaccination policy as both factors are considered separately. The adjustment and additions to the model are depicted in Fig. 8 and the resulting system of ODEs is given by S0 ¼ ð1  cÞp  bSI=N  mS;

557

IT0 ¼ d½nEþ y1 pbEI=N þ pbSI=N þ pð1  q2 Þð1  q1 ÞbVI=N þ y1 pbES I=N þ pbVS I=N þð1  q3 ÞnEv þ y1 pbEv I=N þ py2 bRI=N þ nES  g1 IT  mIT þd=2IT c ; IT0 c ¼ ð1  dÞ½nE þ y1 pbEI=N þ pbSI=N þ pð1  q2 Þð1  q1 ÞbVI=N þð1  q3 ÞnEv þ y1 pbEv I=N þ py2 bRI=N þ nES þ y1 pbES I=N þ pbVS I=N þ ð1  sÞg1 IT  g2 IT c  ðm þ mT ÞIT c  d=2IT c ; R0 ¼ sg1 IT þ g2 IT c þ tE  y2 bRI=N  mR:

ð22Þ

5.1. Basic reproductive number and threshold r  Again, using the method of Diekmann and Heesterbeek (2000) we can compute the reproductive number for the newly expanded model. The model now contains five infected states: E, EV , ES , IT and IT c . After making the necessary adjustments we arrive at 3 2 ð1  pÞð1  cÞðm þ oÞ 7 6 ð1  cq1 Þm þ o 7 6 7 6 ½1  pð1  q2 Þð1  q1 Þmc 7 6 7 6 7 6 ð1  cq1 Þm þ o 7 6 7 6 ð1  pÞoc 7 6 ~ ¼6 Y 7; ð1  cq1 Þm þ o 7 6 7 6 6 dp½mð1  cÞ þ o þð1  q1 Þð1  q2 Þmc 7 7 6 7 6 ð1  cq1 Þm þ o 7 6 6 ð1  dÞp½mð1  cÞ þ o þ ð1  q Þð1  q Þmc 7 5 4 1 2 ð1  cq1 Þm þ o 2 6 6 6 G¼6 6 6 4

0

V ¼ cp  oV  ð1  q1 ÞbVI=N  mV;

n  m  t

0

0 0 dn ð1  dÞn

0

0

0

m  o  ð1  q3 Þn

0

0

0

o

n  m

0

0

dð1  q3 Þn ð1  dÞð1  q3 Þn

dn ð1  dÞn

g1  m ð1  sÞg1

d=2 d=2  g2  m  mT

3 7 7 7 7 7 7 5

and VS0 ¼ oV  bVS I=N  mVS ; E0 ¼ ð1  pÞbSI=N  ðnE þ y1 pbEI=NÞ  tE  mE; Ev0

¼ ð1  pð1  q2 ÞÞð1  q1 ÞbVI=N  oEV ðð1  q3 ÞnEv þ y1 pbEv I=NÞ  mEv ;

ES0 ¼ ð1  pÞbVS I=N þ oEV  nES  y1 pbES I=N þ ð1  pÞy2 bRI=N  mES ;

h T ~ h ¼ 0

0

0

bðmmcq1 þ oÞ

mþo

i :

bðmmcq1 þ oÞ

mþo

~ , we get the basic Recalling once again that R0 ¼  ~ h  G1 Y reproductive number for the ‘‘worst-case’’ model:  ½ðm þ oÞðp  1Þr þ mq1 ðq3  1Þn2 R0 ¼ Kðpm þ n  nr þ nrpÞ þ Kc ðm þ oÞðn  nq3 þ o þ mÞ ½ðq3  r  q1 þ q1 q2  q2 Þp þ r  q1 þ q1 q3  q3 nm2 ðm þ oÞðn  nq3 þ o þ mÞ  pm2 ½ðq1  1Þq2  q1 ðm þ oÞ  noð2rp  2r þ q1 Þm  nr o2 ðp  1Þ þ ðm þ oÞðn  nq3 þ o þ mÞ ð23Þ þ

where K¼

bðg1  g1 dsþ m þdg2 þ d=2þ dmT Þ : ðn þ mÞ½ðm þ mT þ g2 þ ds=2Þg1 þ mðd=2 þ m þ mT þ g2 Þ

To obtain the threshold for detection/treatment of latent tuberculosis under which vaccination increases the stability of the DFE, we set the coefficient of c to 0. Solving for r, we obtain r ¼

Fig. 8. Worst-case scenario for vaccination.

ðm þ nÞðpm  mpq2 þ n  nq3 Þmq1 þ ½nð1  pÞq3 þ pq2 ðm þ nÞm2 nðm þ oÞð1  pÞðm þ o þ n  nq3 Þ m½mpq1 ð1  q2 Þ þ q1 n þ mpq2  þo : ð24Þ nðm þ oÞð1  pÞðm þ o þ n  nq3 Þ

A quick comparison shows that (24) reduces to the previous threshold (19) when o ¼ 0. Using the realistic parameters from TB in the United States, we examine the effect of the model modifications on r . As Fig. 9 depicts, our modifications have a substantial effect on the thresholds. The most significant is that

ARTICLE IN PRESS 558

D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

Fig. 10. r  ¼ 0:15 (white), r  ¼ 0:35 (gray) and r  ¼ 0:65 (black) for o ¼

1 45.

5.2. Comparing vaccine waning with imperfect protection As previously mentioned, distinguishing imperfect vaccination which results from waning of protection and that which results from pure lack of effectiveness can be nontrivial. Hence, it is interesting to examine the effects on the threshold for detection/ treatment of latent TB, r  , that result from considering an imperfect vaccination in the two different manners. For instance, a vaccine that is judged to be 75% effective, can be accomplished in two ways: (1) by offering complete protection but only for a limited time or (2) by offering life-long but imperfect protection. For simplicity, we consider vaccines which only work via one of the three methods described in Section 2.3. Assuming that the vaccine works with efficacy q and provides life-long protection (i.e. o ¼ 0), we wish to find the necessary protective efficacy, which we denote qw , that would result in the same threshold value for r as a vaccine which wanes with time (i.e. o a 0). Using (24), the calculation is straightforward and results in     mþo mþo q3 ðm þ oÞðo þ n þ mÞ w q q2 ; qw : ¼ ; q ¼ qw 1 1 2 3 ¼ m m q3 on þ mðm þ nÞ

Fig. 9. Thresholds plotted for US tuberculosis parameters for various lengths of 1 , vaccine protection: (solid) lifelong protection i.e. o ¼ 0, (spaced dots) o ¼ 50 1 1 1 , (dashed) o ¼ 20 , (dotted) o ¼ 10 . (a) q1 only, (b) q2 only and (c) q3 (dash-dot) o ¼ 30 only.

seen for a vaccine which only protects against slow progression to active TB (q3 only). This is not surprising as two of our implicit assumptions in the original version of the model had the effect of giving additional individuals protection from slow progression. In Fig. 10, we examine the resulting surfaces for various values of r  . Again, we see that the additions to the model have greatly reduced the region in which vaccination results in a more stable DFE.

The expressions for q1 and q2 are intuitive if we recall the expression for the DFE given by (4). For a vaccine with waning rate o, we expect that in the absence of infection a proportion m=ðo þ mÞ of the population will be protected by the vaccine. Thus, in order to preserve the threshold for r, the efficacy of the vaccine must be improved by a factor of ðm þ oÞ=m. This highlights a major drawback in modeling a waning vaccine via our methods which relates to the implication of exponentially distributed sojourns. Under exponential distributions, even a vaccine which last 80 years on average will only be expected to protect slightly more than half of the population in the long term. This is in stark contrast to what one would expect in the realistic setting. In addition to the ideas previously mentioned, this gives another manner in which our ‘‘worst case’’ model for vaccination may overestimate the detrimental effects of universal vaccination.

6. Discussion Nearly 90 years have passed since the Bacille Calmette-Gue rin (BCG) vaccine was first administered to humans (Centers for Disease Control and Prevention, 1996), yet many uncertainties

ARTICLE IN PRESS D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

regarding its use and effectiveness still remain. One of the greatest arguments against mass vaccination has been that BCG can cause positive reactions to tuberculin skin testing and hence interfere with the diagnosis of latent tuberculosis infection. In this work, we have constructed a dynamical system to model the epidemiology of tuberculosis which includes vaccination, treatment of latent TB, treatment of active TB, and immigration in order to study the trade-off that exists between the benefits of vaccination and the ability to detect/treat latent TB. Our model assumes that latent infection in a vaccinated individual is completely undetectable due to the tuberculin reactivity caused by BCG vaccination. This assumption is almost surely too restrictive as evidence suggests that tuberculin reactivity wanes quickly in individuals vaccinated as infants (Menzies, 2000; Menzies and Vissandjee, 1992) and that tuberculin skin testing can still be a useful tool in detecting latent infection in BCGvaccinated individuals (Centers for Disease Control and Prevention, 2000; Garc´ıa-Sancho et al., 2006; Musoke Mudido et al., 1999). Using a stability analysis of the disease free equilibrium and the basic reproductive number, an analytic threshold which dictates when vaccination will increase the stability of the disease free equilibrium was established. The threshold was examined under various scenarios both theoretically and using TB data from the United States. Lastly, adjustments were made to the model which the authors believe represent the strictest possible assumptions against mass vaccination and the resulting threshold was examined. The modeling efforts of this paper suggest that a society may have to detect and treat latent tuberculosis at a fairly high rate in order for the discontinuation of vaccination to result in a more stable DFE. It is, however, very unlikely that this is the case. In the US, the combined detection and treatment success rate ðdsÞ has been reported in the range of 54–83% (see Section 4.3) and one would surely suspect the detection/treatment of latent TB to be much smaller (i.e. r{ds). However, the results of this work are far from conclusive evidence that mass vaccination should be continued in low incidence countries. The policy of vaccinating high-risk individuals only, which is typically the policy that replaces mass vaccination, can not be considered via our dynamical systems approach. In excluding immigration in order to study the DFE, we also lose a key factor in the dynamics of TB in low incidence countries. Also, it is important to note that our estimates of y1 and y2 in Table 1 are such that y2 4 y1 . Hence, by treating latently infected individuals and moving them to the recovered class we actually subject them to higher rates of reinfection. Thus, a strong possibility exists that the treatment of LTBI may actually increase disease prevalence in some settings despite increasing the stability of the DFE. One way in which such behavior could be observed analytically would be by proving the existence of a backwards bifurcation in which the disease persists despite having R0 o1. In the current analysis, the issue is mute as we consider only one infectious individual introduced into a completely susceptible population which negates all reinfection terms. Most importantly, the theoretical threshold we establish only dictates which policy will increase the stability of the DFE and gives no quantitative insight into how much better the policy will perform. Surely, this is a vital component of any major policy decision as a mass vaccination program that requires millions of immunizations to prevent one case of disease would never be a prudent choice. In a subsequent paper (Gerberry and Milner, 2010), we address these issues with a simulation-based analysis of the model. The model is calibrated to WHO for eight countries with widely varying levels of TB burden. A ‘‘best fit’’ analysis as well as

559

extensive uncertainty analysis allow us to better quantify the impact of vaccination in various current TB settings. References Araujo, Z., de Waard, J.H., Ferna ndez de Larrea, C., Borges, R., Convit, J., 2008. The effect of Bacille Calmette-Gue rin vaccine on tuberculin reactivity in indigenous children from communities with high prevalence of tuberculosis. Vaccine 26, 5575–5581. Aronson, N.E., Santosham, M., Comstock, G.W., Howard, R.S., Moulton, L.H., Rhoades, E.R., Harrison, L.H., 2004. Long-term efficacy of BCG vaccine in American Indians and Alaska Natives: a 60-year follow-up study. JAMA 291 (17), 2086–2091. Bhunu, C.P., Garira, W., Mukandavire, Z., Magombedze, G., 2008a. Modelling the effects of pre-exposure and post-exposure vaccines in tuberculosis control. Journal of Theoretical Biology 254, 633–649. Bhunu, C.P., Garira, W., Mukandavire, Z., Zimba, M., 2008b. Tuberculosis transmission model with chemoprophylaxis. Bulletin of Mathematical Biology 70, 1163–1191. Blower, S.M., McLean, A.R., Porco, T.C., Small, P.M., Hopewell, P.C., Sanchez, M.A., Moss, A.R., 1995. The intrinsic transmission dynamics of tuberculosis epidemics. Nature Medicine 1 (8), 815–821. Centers for Disease Control and Prevention, 1996. The role of BCG vaccine in the prevention and control of tuberculosis in the United States. Morbidity and Mortality Weekly Report: Recommendations and Reports, 45(RR-4). Centers for Disease Control and Prevention, 2000. Targeted tuberculin testing and treatment of latent tuberculosis infection. Morbidity and Mortality Weekly Report: Recommendations and Reports, 49(RR-6), pp. 1–54. Cohen, T., Colijn, C., Finklea, B., Murray, M., 2007. Exogenous re-infection and the dynamics of tuberculosis epidemics: local effects in a network model for transmission. Journal of the Royal Society Interface 4, 523–531. Diekmann, O., Heesterbeek, J.A.P., 2000. Mathematical Epidemiology of Infectious Disease. Wiley, Chichester. Dye, C., Garnett, G.P., Sleeman, K., Williams, B.G., 1998. Prospects for worldwide tuberculosis control under the WHO DOTS strategy. Lancet 352, 1886–1891. Frieden, T.R., Sterling, T.R., Munsiff, S.S., Watt, C.J., Dye, C., 2003. Tuberculosis. Lancet 362, 887–899. Garc´ıa-Sancho, F.Ma.C., Garc´ıa-Garc´ıa, L., Jime nez-Corona, Ma.E., PalaciosMart´ınez, M., Ferreyra-Reyes, L.D., Canizales-Quintero, S., Cano-Arellano, B., ´ Ponce de Leon, A., Sifuentes-Osornio, J., Small, P., DeRiemer, K., 2006. Is tuberculin skin testing useful to diagnose latent tuberculosis in BCGvaccinated children?. International Journal of Epidemiology 35, 1447–1454. Gerberry, D.J., Milner, F.A., 2010. A data-based model for tuberculosis: Who should be vaccinating? Preprint. Gomes, M.G.M., Franco, A.O., Gomes, M.C., Medley, G.F., 2004. The reinfection threshold promotes variability in tuberculosis epidemiology and vaccine efficay. Proceedings of the Royal Society of London B 271, 617–623. Hart, P.D., Sutherland, I., 1977. BCG and vole bacillus vaccines in the prevention of tuberculosis in adolescence and early adult life. Final Report of the Medical Research Council. British Medical Journal 2, 293–295. Infuso, A., Falzon, D., 2006. European survey of BCG vaccination policies and surveillance in children, 2005. Eurosurveillance 11, 6–11. Kumar, V., Abbas, A.K., Fausto, N., Mitchell, R.N., 2007. Robbins Basic Pathology, eighth ed Saunders Elsevier. Mayo Foundation for Medical Education and Research. Tuberculosis January. Accessed May 21, 2009 from /http://www.mayoclinic.com/health/tuberculo sis/DS00372/S. Menzies, D., 2000. What does tuberculin reactivity after Bacille Calmette-Gue rin vaccination tell us?. Clinical Infectious Diseases 31 (Suppl. 3), S71–4. Menzies, R.I., Vissandjee, B., 1992. Effect of Bacille Calmette-Gue rin vaccination on tuberculin reactivity. American Review of Respiratory Disease 145, 621–625. Musoke Mudido, P., Guwatudde, D., Nakakeeto, M.K., Bukenya, G.B., Nsamba, D., Johnson, J.L., Mugerwa, R.D., Ellner, J.J., Whalen, C.C., 1999. The effect of bacille Calmette-Gue rin vaccination at birth of tuberculin skin test reactivity in Ugandan children. International Journal of Tuberculosis and Lung Disease 3 (10), 891–895. Porco, T.C., Blower, S.M., 1998. Quantifying the intrinsic transmission dynamics of tuberculosis epidemics. Theoretical Population Biology 54, 117–132. Smith, D., Wiegeshaus, E., Balasubramanian, V., 2000. An analysis of some hypothesis related to the Chingelput Bacille Calmette-Gue rin trial. Clinical Infectious Diseases 31 (Suppl. 3), S77–80. Smith, P.G., Moss, A.R., 1994. Epidemiology of tuberculosis. In: Tuberculosis: Pathogenesis, Protection, and Control. ASM Press. Sterne, J.A.C., Rodrigues, L.C., Guedes, I.N., 1998. Does the efficacy of BCG decline with time since vaccination?. International Journal of Tuberculosis and Lung Disease 2 (3), 200–207. United States Census Bureau, 2008. Population division. International Data Base, Decemeber 15. Accessed January 9, 2009 from /http://www.census.gov/ipc/ www/idb/index.htmlS. Verver, S., Warren, R.M., Beyers, N., Richardson, M., van der Spuy, G.D., Borgdorff, M.W., Enarson, D.A., Behr, M.A., Helden, P.D., 2005. Rate of reinfection tuberculosis after successful treatment is higher than rate of new tuberculosis. American Journal of Respiratory and Critical Care Medicine 171 (12), 1324–1325.

ARTICLE IN PRESS 560

D.J. Gerberry / Journal of Theoretical Biology 261 (2009) 548–560

World Health Organization, 2004. WHO position paper on BCG vaccination. Weekly Epidemiological Record 79 (4), 27–38. Accessed on May 5, 2008 from /http://www.who.int/werS. World Health Organization, 2008. Global Tuberculosis Database. Accessed March 12, 2009 from /http://www.who.int/tb/country/global_tb_database/en/S. World Health Organization, 2007. Department of Immunization, Vaccines and Biologicals. WHO vaccine-preventable diseases: monitoring system. 2007

Global Summary. Accessed on March 12, 2009 from /http://whqlibdoc.who. int/hq/2007/who_ivb_2007_eng.pdfS. Young, D., Dye, C., 2006. The development and impact of tuberculosis vaccines. Cell 124 (4), 683–687. Ziv, E., Daley, C.L., Blower, S.M., 2004. Potential public health impact of new tuberculosis vaccines. Emerging Infectious Diseases 10 (9), 1529–1535.