Tuberculosis Epidemiology at the Country Scale: Self-Limiting Process and the HIV Effects

Background The global spread of the human immunodeficiency virus (HIV) is the main hypothesis behind tuberculosis (TB) positive trends in the last decades, according to modeling studies and World Health Organization Reports (WHO). On one hand, TB (WHO) reports do not explicitly consider a modeling approach, but cover country and global TB trends. On the other hand, modeling studies usually do not cover the scale of WHO reports, because of the amount of parameters estimated to describe TB natural history. Here we combined these two principal sources of TB studies covering TB High Burden Countries (HBCs) dynamics. Our main goals were: (i) to detect the endogenous component of TB dynamics since 1974 for TB HBCs; and (ii) to explore the HIV exogenous effects on TB models`parameters. Methods and Findings We explored the relationship between the TB per capita population rate of change (RI) and the infectious class size following an endogenous/exogenous framework. RI can be affected by intra-population processes (i.e. competition, predation) and exogenous variables like HIV. We found that TB dynamics had always a strong endogenous component, represented by a negative correlation between TB population size and RI, which was captured by the discrete logistic model. Moreover, we explored the HIV exogenous effects on TB models`parameters. We found that overall the TB+HIV logistic model was more parsimonious than TB model alone, principally in the African region. Our results showed that HIV affected principally TB carrying capacity, as expected by the known HIV effects on TB natural-history. We also tested if DOTS (Directly Observed Treatment Short-Course Strategy), poverty levels and BCG (Bacillus Calmette-Guérin) coverage explained the models´ residuals variances, but they did not. Conclusions Since 1974, TB dynamics were categorized in distinct chronological domains, with different dynamics but nearly the same underlying mechanism: a negative relationship between RI and infected class size (i.e. self-limiting). In the last decades, not only HIV spread represented a new TB chronological domain, but it also has been pushing TB carrying capacity (K) to higher levels. TB has a complex natural-history and imposes real challenges to model its dynamics. Yet, we were able to explore and reveal the main drivers of TB dynamics for HBCs since 1974, through a simple approach. Based on our results, we suggest that the endogenous view should be considered as a plausible hypothesis to model and explain TB dynamics and that future World Health Organization reports could include the endogenous/exogenous framework as a supplement to help to decipher the main drivers of TB dynamics and other diseases.


Introduction
Tuberculosis (TB) is an ancient air-borne disease [1], with high incidence and mortality rates among human populations [2][3][4]. Despite the advances as BCG vaccine (Bacillus Calmette-Guérin) and antibiotics, TB started to rise again, mainly in low income countries. Although the appearance of multi-drug resistance TB strains may have also played a role, the global spread of the human immunodeficiency virus (HIV) has been the main driver of recent TB positive trends, especially in Asia and Africa [3][4][5][6][7].
The HIV and TB levels peaked in the 1990s and several countries joined in a global effort to reduce TB burden. The programs Directly Observed Treatment Short-Course (DOTS) strategy, the Stop TB strategy and TB-related Millennium Development Goals established general goals aiming to reduce TB levels by half until 2015. Some countries are closer to these goals than others and a major challenge seems to be the link between HIV and TB, principally in Africa where the level of co-infections can reach 75% of TB patients [4].
The knowledge about TB burden comes from two key sources of information: modeling studies and TB World Health Organization reports (WHO). The usual approach to model TB diseases is the SEIR (Susceptible->Expose->Infect->Recover) type models [8]. Due to its inherent variable natural-history, authors have been modeling TB dynamics with distinct assumptions and types of models that range from unstructured ordinary differential equations to spatially structured models, with many compartments [8][9][10][11][12][13][14][15][16][17][18][19][20][21] (S1 Text). WHO reports usually present global and country scale TB trends, which are avoided by modeling studies, because of the amount of SEIR parameters needed to be estimated to capture TB dynamics.
Here we hope to combine these two principal sources of TB studies, by exploring the TB High Burden Countries dynamics from 1974 until 2012, using simple population dynamic models. Population ecologists have developed a simple framework using principles adopted from system dynamics and biological interactions, which helps to explore the contributions of endogenous and exogenous processes even in potential complex systems. Briefly, density-dependency (endogenous process) occurs when the state of a variable (i.e. population size) is determined, positively or negatively, by its own previous states. Density-independency (exogenous process) occurs when a variable affect another but is not affected by the changes it causes. There is no feedback structure between them [22][23][24][25][26][27][28][29][30]. Climate variables and governmental disease control policies are known examples of exogenous pulse/press perturbations regarding diseases [22]. Both processes may operate to determine the dynamic behavior of a particular population through time. Therefore, a more complete understanding of the dynamics of a population is achieved when both endogenous and exogenous processes are taken into account.
It is widely known that HIV may affect TB transmission rate, re-infection risk, the herd immunity effect of BCG vaccine and progression rates from susceptible to infected class [12,13,19,20]. Hence, as suggested in WHO reports and HIV-TB studies, HIV may be acting as a positive press exogenous force on TB dynamics. This study had two general goals: (i) to detect the endogenous component for TB dynamics since 1974 for TB Higher Burden Countries (HBCs); and (ii) to explore and the HIV exogenous effects on TB model parameters. We also explored the effects of other variables that could affect TB levels mentioned in the WHO reports, including poverty levels, vaccination coverage and DOTS (see below).

Materials and Methods
The complex TB natural-history would require the estimation of a large amount of parameters for deciphering TB dynamics for the HBCs following the SEIR approach (S1 Text). Hence, we deviated from the usual SEIR type models. Here we focused only on TB annually reported cases (i.e. diagnosed, new and relapsed) aiming to reveal general patterns of TB infectious class dynamics. Annual TB reported cases were gathered from WHO reports available at WHO website (www.who.int/en/).
The realized per capita population rate of change is the central stone of the framework used here and is estimated as: where R I is the realized rate of change of the average TB patient in the infectious class over a time interval and I t and I t − 1 are the number of annual new and relapse infectious cases at consecutive years. R I can be affected by intra-population processes (including self-limiting and predator-prey interactions) and exogenous variables, including HIV (22)(23)(24)(25)(26)(27)(28)(29)(30). Based on TB literature, system dynamics principles and disease dynamics theory, the R I of TB infectious class can be expressed as: where R I is an unidentified function (f) of the population endogenous processes (I t − d ), exogenous variables (DOTS t -1 , P t -1 , BCG t -1 , HIV t -1 ) and other variables (V) that may also determine the environment where TB dynamics takes place. Assuming a constant environment (exogenous variables), a negative relationship between R I and I t − 1 (i.e. d = 1) suggests that as the infected portion of the population increases, infected individuals face more difficulties to recruit susceptible individuals, reducing the realized per capita rate of transmission (22). This process resembles the law of diminishing returns (i.e. self-limiting), which has the intra-specific competition as its ecological analogous. It represents a negative first-order feedback between R I and I t − 1 . A negative relationship between R I and I t − 2 suggests that TB population size affects itself with a delay, probably through another variable (i.e. second-order feedback) as the susceptible class, a predator-prey interaction. In this case, a circular behavior is observed between R I and I t − 1 . R I can also be a positive function of I t − 1 . Population size may enhance individual probabilities of finding resources and refuges, which is the "Allee effect". Populations under this phenomenon describe a positive correlation between R I and population size. Also, R I may be independent of the infectious population size. Under a constant environment, population may increase unbounded at constant R levels [22][23][24][25][26][27][28][29].
Since data visual inspection revealed abrupt changes in TB dynamics since 1974, it was necessary to cut TB dynamics in distinct chronological periods. Although the cut procedure is subjective, it is based on data behavior (TB and R I levels and trends). Briefly, if R I was declining and suddenly increased and then started to decline again, we assumed that the time series should be cut prior to the increase and the periods should be analyzed separately (S3 Text).
Despite distinct cut-points may produce distinct results, it would probably not change our general findings, principally because we focused on TB last period of growth (see below).
The TB chronological periods and positive trends imposed a challenge to select the model to fit the data. Model selection was based on the visual inspection of relationship between R I and I t − 1 . Population Dynamic Theory states that: R I should increase with population size when some kind of cooperation exists and affects reproduction and survival and that R I should decline with population size when an intra-specific competition is operating and affecting negatively the reproduction and positively mortality [25]. Whenever a decrease or increase in R I with I t − 1 was observed, we fitted the following discrete models, respectively: where R max is the maximum per capita rate of change, in the previous year, K is the stable population size at equilibrium (carrying capacity) and Q is a non-linear term that reflects the intensity of intra-specific competition around K. E is an unstable equilibrium population size for the discrete Allee model (Eq 4), where, above it, population always grows and below it, population tends to extinction. T represents how fast population deviates from the unstable equilibrium size [25,31]. Most countries exhibited a negative relationship between R I and I t − 1 , suggesting the logistic model (Eq 3) as the first candidate to model and capture TB dynamics [22]. Nevertheless, two countries showed an increase followed by a decrease in R I levels (UR Tanzania and Bangladesh. See below). In these cases we fitted the following Allee-Logistic model: where most parameters are the same as in Eq 3 and Eq 4, except for the absence of the Q and T.
In this formula, we assumed that R I behaved as a quadratic function of the I t − 1 .
The R I -function (Eq 2) shows that R I can also be affected by exogenous variables like the DOTS (DOTS t − 1 ), poverty levels (P t − 1 ), BCG coverage (BCG t − 1 ), HIV (HIV t − 1 ) and other variables (V) which may also determine the environment where TB dynamics takes place. Here, we were specifically interested in exploring and estimating the HIV impact in TB dynamics. As already mentioned, most countries showed periods of logistic growth, and because the logistic model has simple biological interpretations, we can explore the HIV effects on TB model parameters. We took the last period of growth of TB time series for each country and tested three types of HIV exogenous effects [24]: Lateral effect, where HIV affects only TB carrying capacity (Eq 6); Additive effect, when HIV affects both R max and K (Eq 7); and Non-linear effect, where HIV affects the Q parameter (Eq 8).
where g is the linear coefficient that measures the effect of HIV on TB model parameters.
We tested HIV effects on TB time series for distinct time lags (t) to obtain the best model, which was selected based on the coefficient of determination (R 2 ) and the Akaike Information Criterion (AIC). Models with higher R 2 and lower AICs are more parsimonious and hence, were selected. We FItted Eqs 3-8 using the nls library in R using non-linear regression analysis [32]. We used the coefficient of determination (R 2 ) from the above TB models (without HIV) as an estimation of the contribution of the endogenous processes to TB dynamics.
Due to the low degrees of freedom for the last periods of TB growth (see below), it was not always possible to properly estimate the additive (Eq 7) and non-linear HIV effects (Eq 8). Hence, we took the residuals from TB logistic model without HIV (Eq 3) and TB+HIV lateral model (Eq 6) and explored whether the non-HIV variables (DOTS, poverty levels and BCG coverage) individually explained the remaining variance, with Pearson´s product moment correlation coefficient. Detection of TB cases (DT) and treatment success (TS) are two key components of the DOTS strategy. WHO adopted a 70% detection rate and treatment success of 85%, expecting that these levels would lead to a 5% year reduction in the TB incidence rates. We collected from WHO web page the annual countries percentage of TB detected cases (DT) by sputum-smear microscopy and the percentage of these cases successfully treated (TS) From 1974 and 2012, DT and TS for most of all HBCs were below the WHO levels [4,5]. We used per capita GDP as a proxy of poverty levels. Per capita GDP was estimated by the ratio between country GDP and human population size obtained from the World Bank web site (www.worldbank. org). BCG coverage and HIV reported cases were obtained from Gapminder web site (www. gapminder.org/data/). DT and TS time series usually were available from 1990 on, and HIV and BCG coverage since 1980.
In the main text, we present the TB dynamics and analyses for South Africa, Nigeria, Kenya, Mozambique, UR Tanzania, Zimbabwe, Brazil, Bangladesh, Cambodia and China, as a sample of TB dynamic behaviors and statistical results observed for the HBCs countries. For a complete appreciation of the results found in this study, see the supplementary material (S1 Table, S1 and S2 Figs).

Results
Model fits revealed that TB dynamics was always dominated by a first order negative correlation between R I and population size, which was captured by the discrete intra-specific competition model (Figs 1 and 2, S1 and S2 Figs, Table 1, S1 Table).
We also found other behaviors including R I increase trend (e.g. UR Tanzania and Bangladesh), exponential growth (e.g. Mozambique), and circular R I behavior (e.g. Myanmar and Nigeria), which suggest a predator prey interaction (S1 and S2 Figs, S1 Table). These behaviors were mostly restricted to the initial phases of TB time series, while the last periods were dominated by the logistic growth (Fig 2, S2 Fig, S1 Table). Myanmar and perhaps Nigeria showed the most drastic changes in R I behavior (Figs 1 and 2, Table 1, S1 and S2 Figs, S1 Table). They started with a circular R I behavior, while in the last years these countries displayed a logistic growth pattern. Curiously, both circular behaviors had a lag of approximately 7 years.
R max , K, Q, E and g are the logistic models parameters as in the text. R 2 is the coefficient of determination which measures the endogenous component for TB dynamics. AICs is the Akaike Information Criteria. Lag refers to the year in which TB+HIV AICs showed the lowest value. R max was estimated directly from the data and fixed, except for Vertical and Non-linear    (Figs 1 and 2, S1 and S2 Figs). The coefficient of determination (R 2 ) measures the endogenous component of a dynamic (i.e. how much R I is determined by population size). In most periods of logistic growth, R 2 was high suggesting that TB dynamics had always a strong endogenous component (Table 1 and S1  Table). In Africa, the average increase in K between the first and the last TB logistic growth periods was of 878%. Outside Africa it was of only 45.8%. R max increased 21.9% in Africa, whereas it decreased 173.2% outside Africa. These general trends suggest an overall enhancement of the carrying capacity and an usually negative change in the maximum per capita rate of change, R max ( Table 2). Introducing HIV as an exogenous variable in the last TB logistic period resulted in more parsimonious models than TB alone, based on AICs values (Table 1, S1 Table), principally in the African region. In general, R 2 did not improve after including HIV into the TB model (Table 1, S1 Table), which can be explained by TB strong endogenous component. We also found some unexpected results. For example, based only in AICs results, in Uganda and Tanzania the lateral TB+HIV model was less parsimonious than TB model alone (Table 1, S1 Table), which can be explained by HIV shorter time series (not showed here). We were not able to test for time-lags greater than five years, which might have resulted in better AICs for these countries, where the link between HIV and TB is well known [7]. On the other hand, the TB+HIV model improved R 2 in about 20% in both Uganda and Tanzania (S3 Text) supporting the inclusion of HIV lateral effect on TB logistic model.
Outside Africa, the TB+HIV models did not improve R 2 , neither decreased AICs, as observed for Bangladesh, Pakistan, Philippines and Vietnam (Table 1, S1 Table). For China, including HIV produced higher R 2 and lower AIC than TB model alone. Yet, it produced a negative relationship (g, the HIV coefficient) between HIV levels and TB carrying capacity (K). We suggest that the statistical link between these diseases in China should be taken with caution. Overall, DT (detection rate), TS (treatment success), BCG coverage (Bacillus Calmette-Guérin vaccine) and per capita GDP (as a proxy of poverty) had a poor explanation power regarding TB and TB+HIV (lateral) model residuals (S2 Table).

Discussion
Two major findings deserve to be highlighted: (i) besides some exceptions, since 1974, TB dynamics presented a strong endogenous component, a first order feedback between the TB population size and its population per capita rate of change; and (ii) TB models with HIV exogenous effects were more parsimonious than TB model alone, especially in the African region, These findings suggest a strong endogenous component and supports the known HIV effects on TB natural-history and dynamics.
We found that the simple logistic discrete model (for the infectious class) can capture the TB behavior for a particular period of time. Although the R I decline trend may be also the result of improvements in health programs and life standards, the sequential logistic periods, the R 2 levels between R I and I t − 1 , together with the weak correlations between non-HIV variables and model residuals, suggest that the endogenous view should be considered as a plausible hypothesis to explain and model TB infectious class dynamics.
The sequence of TB logistic periods of growth suggests a succession of population thresholds levels, from which TB is freed by its exogenous control strategies. Hence, the sequence of logistic growth patterns can be interpreted as a series of density-dependency chronological domains. Since 1974, TB has suffered changes on its equilibrium surface, changing only the parameter values (principally K), not the biological process underpinning its dynamics (i.e. diminishing returns, self-limiting). This resembles the fifth principle of population dynamics suggested by Berryman [23,25], which is based on the "Minimum law" of Leibig. At a particular period of time, the population (i.e. the infectious class), at equilibrium, is limited by the self-limitting process for the lower resource available in the environment (i.e. susceptible class). During the last decades, HIV appearance not only represented a new TB period, but it has also been continuously pushing TB model parameters to a new globally stable landscape. Oncoming HIV and TB vaccines, in addition to new prophylaxes, combined with control strategies will probably change TB dynamic properties again, leading to a new period of a global decline [3][4][5][6][7].
The changes between the first and last TB model (Table 2) and the results of our TB+HIV models suggest that HIV probably acted as an exogenous forcing on TB dynamics. These results are in agreement with the expected HIV effects on TB dynamics. HIV simplifies the complex TB natural history. As HIV load increases, individuals are more susceptible to be infected by TB. When individuals are co-infected, HIV accelerates the progression rates from TB exposed to infectious stage. It has been estimated that one third of the global population is in the Exposed stage of TB [3], representing a large pool of individuals that can rapidly progress to active TB [9-13, 19, 20]. Yet, a co-infected individual is not necessarily more contagious than a HIV negative individual with active TB. These known HIV effects on TB may explain the disproportional K changes compared to R max inside Africa, where the level of co-infection may reach 75% of TB cases [3].
Despite the main effect of HIV and the lack of statistical relationship between residuals and non-HIV variables, we cannot rule completely out the importance of these variables. Both TB and HIV are related to poverty and the other variables (DT, TS and BCG) may be correlated to poverty levels through governmental funding [33][34][35][36][37][38][39][40] (S4 Text). Almost all countries analyzed here were until recently European colonies. They share low per capita income levels, have deficient health infrastructures and suffer from internal political instabilities, which create a perfect environment to the spread of HIV. Also, most of the countries have suffered from armed conflicts, which may increase disease burdens due to their impacts on health services, food production and distribution, migrations and refugees overcrowding [41,42]. For example, the TB initial peaks in Afghanistan and in Viet Nam may be correlated to the war conflicts (S1 and S2 Figs).
TB reported cases may not exactly measure TB infected class and may be also a function of changes in diagnosis, reporting practices and case findings [12]. Until the spread of HIV, some positive trends could be explained by improvements in detection rate and not by TB threshold levels and the endogenous process, as suggested above. Despite these drawbacks, we used TB reported cases as a proxy of TB burden, which allowed describing empirically and theoretically TB time series since 1974, and we were also able to capture the known regionally HIV effect on TB burden in Africa. A finer scale analysis would reveal details about the contribution of case finding and reporting practices vs TB endogenous behavior to the periods prior the HIV spread. A finer scale analyses may also describe more complex scenarios and the effects of other variables including MDR-TB strains, HIV strains, the role of susceptible, exposed and recovered pools, urbanization process, ART availability, BCG waning immunity with human increase life expectancy, colinearities between DT and TS, among others [12][13][14][15][16][17][18]43]. A finer scale analysis was not the scope of this study though. Our goal was to apply population ecology theory to explore TB trends and estimate HIV effects. We found distinct TB density-dependency domains captured by the logistic discrete model, and we found that HIV had an exogenous effect on TB carrying capacity, principally in Africa. The framework is based on the realized per capita population rate of change (R I ), which is surrounded by plausible ecological principles, and hence, a good starting point to explore TB dynamics. Any country interested in its TB epidemic may use our results as hypotheses and move forward. For example, outside Africa, governments can substitute HIV lateral effect for DT, TS, BCG or per capita GDP, compare and select the best model based on statistical results and their knowledge about TB epidemic. In addition, it is also possible to test the endogenous and exogenous contributions at distinct scales inside the country (Regions, cities, etc), principally for large territory countries like Brazil. Any government may disentangle R I in its components (new per capita infections and per capita mortality), explore which of them are more important for R I trends and also explore the contributions of the endogenous and exogenous processes.
This study suggests that governments' efforts in controlling TB epidemic were not enough to avoid the impact of HIV in TB, principally in African region, in agreement with many other studies and reports [3,4]. This does not mean that control policies are completely inefficient, since many TB related deaths are annually avoided by improvements in health care and the impact of HIV would be even greater than observed. This study also shows that governments can use simple tools to explore and identify the mechanisms underlying TB dynamics and other diseases. Hence, we suggest that this framework could be included as a supplement into WHO reports with minimal cost and time demanding efforts, which could provide insights, hypotheses and may allow testing and estimating the relationship between TB and other variables.   Table. Pearson´s product moment correlation coefficient between TB alone, TB+HIV lateral models`residuals and non-HIV variables.

Supporting Information
(DOC) S1 Text. TB natural-history, SEIR models and interactions between TB infectious and susceptible classes.