Modeling the potential impact of host population survival on the evolution of M. tuberculosis latency.

Tuberculosis (TB) is an infectious disease with a peculiar feature: Upon infection with the causative agent, Mycobacterium Tuberculosis (MTB), most hosts enter a latent state during which no transmission of MTB to new hosts occurs. Only a fraction of latently infected hosts develop TB disease and can potentially infect new hosts. At first glance, this seems like a waste of transmission potential and therefore an evolutionary suboptimal strategy for MTB. It might be that the human immune response keeps MTB in check in most hosts, thereby preventing it from achieving its evolutionary optimum. Another possible explanation is that long latency and progression to disease in only a fraction of hosts are evolutionary beneficial to MTB by allowing it to persist better in small host populations. Given that MTB has co-evolved with human hosts for millenia or longer, it likely encountered small host populations for a large share of its evolutionary history and had to evolve strategies of persistence. Here, we use a mathematical model to show that indeed, MTB persistence is optimal for an intermediate duration of latency and level of activation. The predicted optimal level of activation is above the observed value, suggesting that human co-evolution has lead to host immunity, which keeps MTB below its evolutionary optimum.

Introduction

Tuberculosis (TB) is an infectious disease caused by the bacterium Mycobacterium tuberculosis (MTB). MTB has been infecting humans for a long time, at least millenia and likely even longer [1]. As such, MTB had to evolve evolutionary strategies that allowed it to persist in small groups of human hosts. An interesting question is if one can see signatures of this potential evolutionary adaptation to small groups of human hosts in MTB's “life history” today.

Upon infection with MTB, most hosts enter a latent state. Those hosts do not show signs of disease but do harbor MTB, which can activate and lead to disease at any future time [2]–[4]. It might take a long time before activation occurs, and the majority of TB infected hosts die from other causes besides TB without ever developing disease [5]. Hosts latently infected with MTB cannot infect others. Therefore, at first glance, latency does not seem to be beneficial for MTB. One possible explanation for the long latency and the fact that activation to disease only occurs in a fraction of hosts is based on the human immune response [6]–[8]. It is known that a competent immune response is needed to contain infection and avoid disease, as dramatically illustrated in HIV infected hosts with weakened immune responses, who activate at much higher rates [9]–[11]. This suggests that the co-evolutionary arms race between MTB and its human hosts has lead to a host immune response that can successfully contain MTB in a state of suboptimal fitness [12], [13]. The fact that there seems to be local adaptation on both the pathogen and host side lends support to this idea [13]–[17].

Although host immune pressure is a plausible explanation for reduced activation rates, there are also arguments against this idea. Pathogens like MTB replicate rapidly (compared to their human hosts) and often reach large population sizes; both features are known to foster rapid evolution, especially under strong selective pressure [18], [19]. This is at vivid display in the evolution of drug resistance for MTB and many other pathogens [20], [21]. Therefore, one could argue that if long latency and low activation rates were evolutionary strongly suboptimal for MTB, evolution would have led to higher rates of activation. Instead, as has been suggested previously [22], [23], long latency and low activation rates might be strategies that are evolutionary beneficial to MTB by increasing its fitness.

There is no single way to quantify the fitness of an organism. For infectious diseases, an often used measure of fitness is transmissibility, usually defined by the reproductive number, R 0 [24], [25]. It can be shown that in direct competition of two strains (under equilibrium conditions), the strain with the higher R 0 outcompetes the one with the lower R 0 [25]. However, in the absence of direct competition, strains with lower R 0 might sometimes be more advantaged as they can better avoid local extinction and therefore win through indirect competition against strains with higher transmissibility [26]–[29]. Therefore, an alternative way to capture fitness of pathogens is by quantifying their ability to persist in a host population.

The impact of transmissibility and persistence on overall fitness has been an active area of research [30]–[34]. The importance of different measures of fitness such as these depends on the situation. Since MTB is an ancient organism that has infected humans for millenia [12], [35]–[38], it had to evolve strategies that allowed it to persist in relatively small, spatially structured host populations. Under such circumstances, fitness benefits due to non-extinction might have had an important impact on overall evolutionary dynamics [28], [39], [40].

It has been previously proposed that a prolonged latency might have been one persistence strategy [23], [41]. Here, we use a mathematical model to explore this idea. We introduce a measure of persistence of MTB in a population and investigate how well MTB can persist as a function of the latent period. We find that TB persistence is optimal for an intermediate duration of latency and level of activation. We also find that the optimal level of activation is above the observed value, suggesting that host immunity plays some role in keeping MTB below its optimal level of fitness.

Methods

Mathematical model description

We use a compartmental mathematical model formulated as a set of ordinary differential equations to describe the population-level infection dynamics of TB. The model is shown schematically in figure 1. Our model is similar to other recently studied TB models [42], [43]. We consider 3 types of hosts (compartments) in our model: susceptible hosts, S, latently infected hosts, L, that harbor MTB but are not infectious and show no signs of disease, and infectious hosts with active disease, I. We keep the model simple and do not distinguish between features such as age-related differences (e.g. children versus adults). Such additional details could be included in further more detailed models. The total population size is N = S+L+I. New hosts enter the system at a maximum rate per person, this rate saturates once the population reaches some maximum level, N. All hosts die due to causes other than TB after an average lifespan of 1/m years. Uninfected hosts can become infected through contact with an infectious host at rate b. The infection process is modeled in a density dependent manner [43], [44]. After infection, a small fraction f of hosts rapidly develop disease (fast progression) [11], [45], [46] and move into the active disease compartment, I. The majority of hosts enter the latent state, L (slow progression). Latently infected hosts can convert to infectious, diseased hosts later in their life at rate a (activation) or through reinfection, with the chance of disease due to reinfection reduced by a factor k [47]–[49]. Infectious, diseased hosts either die due to disease after on average 1/m years, or regress and return to the latent stage at rate w. Following previous models, we assume that recovered individuals do not fully clear the infection but instead return to the latent stage [43], [50]–[52]. Table 1 summarizes the model variables and parameters and provides references for the parameter values used for most of our analysis. The model equations are

Figure 1

Flow diagram for the TB transmission model.

Susceptible hosts (S) are born at rate with birth rate saturating at high population density, N. Hosts in each compartment die at a background mortality rate m. Infection of susceptibles occurs at rate b upon contact with TB infectious hosts (I). After TB infection, a fraction f of hosts develop active TB in a short time (fast progression), the remainder enter the latent stage. Latent hosts (L) can activate and develop TB disease sometime later at rate a (slow progression). Latent hosts can also develop active TB after re-infection with TB, with probability of this occurring compared to fast progression of susceptibles reduced by a factor k. Infectious, diseased hosts (I) might regress to the latent stage at rate w. Hosts with active disease die due to diseased induced mortality at rate m.

Table 1

Initial conditions and parameter values.

Symbol	Interpretation	Value	Source and Comments
Nm	maximum population size	1	population size normalized to 1
S 0	initial susceptible hosts	0.6	calculated as , to obtain steady population size in the absence of TB
L 0	initial latent hosts	0	arbitrary choice
I 0	initial infectious, diseased hosts	S 0/1000	one infected per 1,000 hosts
	maximum birth rate	0.05 per year	50 births per 1,000 hosts, representing a high birth rate scenario
mn	natural mortality rate	0.02 per year	assuming a life-span of 50 years for healthy hosts
md	disease-induced mortality rate	0.33 per year	assuming 3 years life-span for untreated diseased hosts [81]
b	rate of transmission	10 per year	[82]–[84]
w	rate of regression	0.25 per year	return from active TB stage to latency [81]
f	fraction of TB infections that lead to disease via fast progression	0.1	[11], [46], [50], [84]
k	reduction of fast progression rate upon reinfection of latent hosts	0.5	[45], [46], [49]
a	rate of TB activation in latent hosts	varied

Initial conditions of model variables and values of model parameters. These values are chosen for all simulations unless indicated otherwise.

We implemented all model simulations in R [53]. The code to reproduce the simulations described here is available from the author's webpage at http://handelgroup.uga.edu/resources.htm.

Results

Persistence of MTB as a function of latency

Our main outcome of interest is the potential of MTB to persist in a host population and how persistence potential depends on the duration of the latent period and fraction of latent TB hosts that activate. While pathogens that have an environmental stage as an important component of their transmission cycle can persist even in the absence of any infected hosts [54]–[58], for pathogens where direct transmission is the main important component, persistence (non-extinction) requires the continued presence of infected hosts. For MTB, extinction occurs if no more latently infected hosts, L, and infectious, diseased hosts, I, are present. Persistence is more likely (i.e. extinction is less likely) as the number of infected hosts in the population increases [59]–[61]. Since for MTB, only a fraction of latently infected hosts will become infectious and contribute to transmission, a better measure of persistence is given by the quantity P, defined as , where I and L are the number of infectious and latently infected hosts, and represents the fraction of latently infected hosts that will develop TB disease, become infectious and are able to transmit. For our model, is given by the ratio of the rate of latently infected hosts that move on to develop active TB, , divided by the total rate at which latently infected hosts leave the latent stage, a+kfbI+m, i.e. .

Persistence as defined by P, especially at the steady state, provides a useful measure for the ability of MTB to not go extinct in the population (see our comparison with a stochastic model below). We run simulations of our model for different values of the activation rate, a, and record I and L at steady state and use this to compute P as function of a. Figure 2A shows that optimal persistence is achieved at intermediate rates of activation. Figure 2B shows individually the three components that make up P. The number of infectious hosts at steady state, as well as the fraction of hosts activating, , increase with increasing activation rate. The number of latent hosts first increases and then decreases with larger activation rate. The combination of these three quantities leads to a maximum for persistence P at intermediate values.

Figure 2

Persistence as function of activation.

A) Persistence, P, as a function of activation rate. B) Latent and infectious hosts, and , and the fraction of activation, , as a function of activation rate. The left axis applies to L and I, the right axis to . C) Persistence as function of fraction of hosts that activate. Also indicated in the figures is the optimal level of persistence, P 0, and the values for activation rate and fraction of activators at which this optimum occurs, and . All parameters are as given in table 1.

Optimal persistence at an intermediate rate of activation also implies that instead of having every latent host activate and become infectious, it is beneficial for the pathogen to let some infections “go to waste” by way of latent hosts dying before they become infectious. This helps with overall persistence and is worth the “loss” of a fraction of latent hosts due to natural death before they activate and are able to transmit. Figure 2C shows persistence as a function of the fraction of hosts that activate, . The figure also illustrates another interesting point: The optimal fraction of hosts that activate is slightly above 50% given the chosen model parameters. This is higher than the ≈10–20% observed [62], [63], suggesting that MTB is not able to achieve the activation rate that would optimize its persistence. This might be attributable to the host immune response playing a role at reducing activation.

The impact of parameter value uncertainty on optimal persistence and activation

To investigate how sensitive results reported in the previous section are to changes in parameter values, we performed an uncertainty analysis and sampled the model parameters using Latin Hypercube Sampling [64]–[67]. For each parameter, we considered uniform distributions ranging between 0.5 and 2 times the base parameter value shown in table 1. For each parameter sample we computed the values P 0, a 0 and , i.e. the optimal level of persistence and the activation rate and fraction at which the optimum occurs. Figure 3 shows the distribution of those values. Of most interest, for almost all parameter samples the fraction of hosts that activate remains well above the 10% found for MTB in the field.

Figure 3

Uncertainty analysis of model results.

Impact of model parameter variations on A) optimal persistence, P 0, B) optimal rate of activation, a 0, and C) optimal fraction of activation, . Model parameters were sampled in a range of 0.5–2 times the base parameters listed in table 1. Sampling was done using a latin hypercube approach with 1000 samples. Samples for which parameter combinations lead to a biologically unreasonable scenario, specifically natural mortality rate above birth rate, were discarded. For each sample, persistence as function of activation rate and activation fraction was determined (as shown in figure 2) and from this the optimal values for persistence and activation were obtained.

In text S1, we show further results by exploring how changes in each model parameter individually affect optimal persistence and activation. As expected, increased population size (though increased carrying capacity parameter or increased birth rate) leads to improved persistence, while decreased population size (through higher death rates) reduces persistence. Other model parameters have less impact on changes in persistence. The optimal fraction of activators as any of these parameters are changed individually is between 40%–70%, again above the value observed experimentally.

We also explore in text S1 how changes in the model structure, specifically the assumption of exponentially distributed natural lifespans, affect our results. We find that it shifts the persistence curve shown in figure 2 slightly, without affecting the overall results and still giving an optimal fraction of activators around 50%.

The deterministic persistence measure is well reproduced with a stochastic model

Our persistence measure, P, is derived from a deterministic model. Of course, non-persistence, i.e. extinction, is an inherently stochastic process. While it is generally well known that larger population sizes lead to less chance of extinction [59]–[61], it is useful to directly test our deterministic measure with a stochastic model. To that end, we implemented the differential equations as a compartmental stochastic model, with transition rates of the deterministic model becoming transition probabilities [25], [68]. We simulated the stochastic model using an efficient form of the Gillespie algorithm (the adaptive- leap method as implemented in the R package adaptivetau [69], [70]). Starting at the deterministic equilibrium state, we simulated the stochastic model for a fixed number of years and record the fraction of simulations for which at least one infectious or latent individual was still present at the end of the simulation. For the stochastic model, there is no discounting of the latent hosts by a factor . Instead, persistence or extinction is a binary event, with extinction defined as no more latent or diseased infected hosts present and persistence if at least one of these hosts was still present. Figure 4 shows that despite this difference between P and the stochastic simulation, the fraction of simulations for which persistence was found in the stochastic model has a very similar functional shape as our deterministic persistence measure, P, with the optimum, P 0, occurring at more or less the same rate of activation.

Figure 4

Comparison of persistence with results from a stochastic model.

The solid line shows persistence, P, as determined from the deterministic model, symbols show result from a stochastic model. For the stochastic simulation, we started at the deterministic steady state, simulated the model 10,000 times for 1000 years each and counted the numbers of simulations for which at least one latent or infectious host was still present after 1000 years. A population of size 100 was used, all parameter values are as reported in table 1. Note that because the absolute magnitude of P is arbitrary and scales with population size, N, to allow comparison with the stochastic model we rescaled P to be between 0 and 1.

In Text S1, we use the stochastic model to investigate how population growth or decline influence P. We find that persistence improves in the presence of a growing population and worsens if the population declines. The shape of the persistence curve does not change in any important way.

Persistence during epidemic cycles

We have so far focused on MTB persistence at the endemic state. Equally important for pathogens is the ability to persist after introduction to a newly susceptible population. During TB's evolutionary history, it likely got introduced and re-introduced repeatedly into small susceptible subpopulations (e.g. new tribes/villages). In the evolutionary context, persistence upon introduction into a fully susceptible population might therefore have played an important role. Nowadays, many areas of the world where TB is very rare consist of largely susceptible individuals – though for those groups the persistence idea explored here is likely not too important. Upon entering a fully susceptible population, pathogens usually cause an epidemic outbreak, depleting the number of susceptibles. The outbreak is often followed by a fade-out of the disease once most susceptibles have been depleted. Extinction often occurs during this fade-out. For the pathogen to not go extinct, it needs to persist long enough until the number of susceptibles has built up again, usually leading to consecutive smaller outbreaks until the endemic state is reached [71], [72]. We can quantify persistence during epidemic cycles by evaluating our expression for persistence not at the steady state but instead at the overall minimum occurring between introduction of the disease in a susceptible population and eventual attainment of the endemic equilibrium, i.e. we determine the overall minimum . Figure 5 shows the distribution of optimal persistence, optimal activation rate and optimal fraction of activating hosts for P, using the same parameter sampling approach as for figure 3. Comparing the results with those found earlier for P at the steady state (figure 3), one sees that the results are very similar. The main reason for the similarity is that TB has a relatively “slow” disease dynamics [73], without pronounced outbreak peaks and minima. Therefore, for most parameter values, the disease reaches the steady state without a large contraction after the first outbreak, leading to essentially the same results as for the steady state.

Figure 5

Optimal persistence, P, activation rate and fraction for the non-steady state model.

Everything as explained in caption for figure 3.

Again, as for P above, we show additional results for P as function of individual parameters in Text S1. The results are almost identical to those obtained for P, for the reason just explained.

Discussion

Recent in-depth studies of MTB genetic sequences have shown a wide diversity between strains [17]. While harder to determine, there is also accumulating evidence that this genetic diversity results in phenotypic diversity [17], [74], suggesting that MTB evolution is shaped by local selection pressures. A recent study indicates that disease activation rate differs between lineages, suggesting that this phenotype is under evolutionary selection [75]. We used a mathematical model to investigate the role of activation rate and latency duration on the ability of MTB to persist in a host population. Our results support the previously proposed idea that the prolonged latency observed for TB infections might provide MTB an evolutionary advantage, namely improved persistence in a host population [23], [41]. We found that an intermediate rate of disease activation is optimal for persistence.

Interestingly, our model suggests that for optimal persistence, the fraction of hosts that eventually become cases is in the 20%–80% range with approximately 50% as the most likely value. This range is higher than the ≈10% generally seen for TB, suggesting that the host immune response plays some role in keeping TB disease in check, lowering activation rate below what would be the evolutionary optimal level for the pathogen. It is likely that some level of co-evolution between pathogen and host occurred and that humans who have been exposed to MTB for a long time evolved some level of resistance that potentially prevents MTB from reaching its evolutionary optimal activation rate [12]. This also agrees with the observation that upon contact with novel MTB strains, some populations have been shown to experience much higher rates of disease than the usually observed 10% [76]. It is further interesting to note that the optimal rate of activation we find in our study is similar to values reported for HIV positive TB patients [77]–[79]. While the simplicity of our model is a caveat in interpreting this agreement too quantitatively, we believe our results provide another suggestive indication that the host immune response is responsible for keeping MTB activation rates below a value that would be optimal for MTB, and once the immune protection fails, as in HIV infected hosts, MTB activates at rates close to its optimum value.

In summary, our results suggest that an intermediate level of activation from latency to disease is optimal for MTB persistence, that the optimal level depends on the detailed pathogen, host and environment characteristics, and that it tends to be higher than the observed value, suggesting an important role for the immune response to keep MTB in check. While increasing activation rates beyond the optimum to reduce MTB persistence is not a suitable goal from a public health perspective, a reduction in activation rate is much more promising. This would lower the number of hosts with disease, and thereby reduce incidence and prevalence for TB cases and at the same time reduce persistence potential. Potential TB vaccines currently under consideration might help us to achieve such a reduction in activation rate [80].

Additional results and explanations.

(PDF)

Click here for additional data file.