INTRODUCTION
Recently, the World Health Organization (WHO) documented the diagnosis of nearly 8·7 million new cases of tuberculosis (TB) in 2011 with an estimated 1·4 million deaths from TB in the same year [1]. Although the absolute number of TB cases has been falling since 2006, TB remains a leading cause of global morbidity and mortality in that one-third of the human population is infected. The bacterium Mycobacterium tuberculosis that causes TB is generally spread in airborne droplets when people with active disease cough or sneeze.
Multidrug-resistant (MDR) TB has been documented in 114 countries and regions worldwide and has emerged as a global public health issue [2]. MDR TB is caused by strains resistant to at least isoniazid and rifampicin, the two principal first-line drugs used in combination chemotherapy [Reference Dye3]. Treatment for MDR TB patients requires use of second-line drugs for at least 24 months [Reference Iseman4]. Therefore, MDR TB is increasingly becoming a serious threat to TB control, and the recognition of extensively drug-resistant TB has further highlighted this threat [Reference Shah5].
In Taiwan, TB has always had the highest incidence rate in all communicable diseases. The incidence rate of TB in Taiwan was 62·0–74·6/100 000 population for the period 2002–2008, respectively [6]. About 149–164 new MDR TB cases were reported in Taiwan during 2007–2010 [7]. Although MDR TB represents only 1·2% of total new TB cases in Taiwan, controlling MDR TB is challenging because it is becoming increasingly hard to diagnose and treat.
Previous studies have reported a strong association between seasonality and TB in temperate, tropical, and subtropical regions [Reference Luquero8–Reference Nagayama and Ohmori12]. Their findings indicated that seasonal peaks of TB cases generally occurred at the end of winter and at the beginning of the spring [Reference Luquero8–Reference Willis10] or summer seasons [Reference Leung11, Reference Nagayama and Ohmori12]. They also implied that indoor activity in winter, diagnostic delays, seasonal change in human immunity, and an association with sunlight and vitamin D levels might play the crucial factors in affecting TB seasonality [Reference Willis10, Reference Nagayama and Ohmori12]. Thus, there are strong interactions between the effects of seasonality and the effects of weather because seasonality might influence the season of emergence of TB, making both sensitive to additional stress such as climate change.
The simplest mathematical model for modelling drug resistant (DR) and MDR TB epidemics came from Blower & Gerberding [Reference Blower and Gerberding13]. However, their model did not consider the contribution of exogenous re-infection to the overall disease incidence. Over the past two decades many expanded and sophisticated models have been used to predict the future burden of DR and MDR TB. Dye & Williams [Reference Dye and Williams14] and Cohen & Murray [Reference Cohen and Murray15] developed the MDR TB model, incorporating re-infection at a reduced rate by partial immunity applying to latent or recovered individuals. In view of these models, it is recognized that the assumptions about the relative fitness (RF) of MDR strains play a crucial role in describing MDR dynamics [Reference Cohen and Murray15, Reference Luciani16].
The seasonality and meteorological factors might influence the dynamics of TB; however, little is known about the time-series dynamics of TB in the (sub)tropical Taiwan region taking into account seasonal patterns and weather effects. The purpose of this study was to examine seasonal population dynamics of TB and to estimate implicitly the MDR TB infection risk in the Taiwan region. We sought to develop a RF-based MDR TB model built on previous models, seasonal transmission to predict the potential impact of seasonality on TB and MDR TB transmission. A probabilistic model was also developed to estimate the site-specific infection risks of TB and MDR TB in Taiwan regions. We anticipated that our study could provide insights into the dynamics of TB/MDR TB transmission, which would inform the uncertainty over the different regions for populations at high-burden risk.
MATERIALS AND METHODS
Study data
Monthly-based disease burden of TB data in Taiwan were obtained from the Centers for Disease Control of Taiwan (Taiwan CDC) (http://www.cdc.gov.tw/) for 2005–2008. We geographically divided Taiwan into four regions (northern, central, southern, eastern) in order to estimate TB incidence rates. The annual disease burden of MDR TB for each county was adopted from the Taiwan tuberculosis control report [17] and Taiwan CDC national notified disease surveillance system [7] for 2006–2010 to calculate MDR TB incidence rates. The monthly-based mean temperature, relative humidity, and rainfall intensity were adopted from the Taiwan Central Weather Bureau (http://www.cwb.gov.tw/) for 2005–2008.
We found that the average incidence rates during 2005–2008 in northern, central, southern, and eastern regions were 57·55 ± 3·88/100 000 67·32 ± 6·27/100 000 81·52 ± 5·18/100 000 and 111·55 ± 19·45/100 000 population, respectively. Moreover, we also found that the higher TB epidemic areas appeared in Hwalien and Taitung counties in east Taiwan and Pingtung county in southern Taiwan.
The annual incidence rates of TB and MDR TB in Hwalien county were 119·72 ± 13·84/100 000 and 4·91 ± 2·37/100 000 respectively, whereas the annual TB and MDR TB incidence rates were 103·39 ± 9·82/100 000 and 2·82 ± 1·79/100 000 respectively, in Taitung county. The annual incidence rates of TB and MDR TB in Pingtung county in south Taiwan were 105·98 ± 10·36/100 000 and 0·97 ± 0·34/100 000 respectively, whereas Taipei city in north Taiwan experienced the lowest incidences of TB at 48·46 ± 3·39/100 000 and MDR TB at 0·43 ± 0·11/100 000. We thus used TB epidemic data of Hwalien, Taitung, Pingtung counties, and Taipei City to investigate the seasonal transmission dynamics of TB and to estimate implicitly the MDR TB infection risk.
To model drug resistance dynamics, the data of RF of DR strains have to be determined. One of the methods to measure RF of resistant strains is based on the results of genotype clustering studies in that a cluster can be defined as two or more cases having an identical DNA fingerprint [Reference Dye3]. Based on the genotype clustering method, RF can be estimated by calculating the odds ratio as RF = (C R/G R)/(C S/G S), where C R, C S, G R, and G S are the numbers of resistant (R) and sensitive (S) cases that appear singly (G) or in clusters (C) [Reference Dye3].
García-García et al. [Reference García-García18] have recently provided a valuable genotype clustering study that can be used to estimate RF of MDR TB. They found that the overall rate of resistance was 28·4%, in that 10·8% had MDR TB. DNA fingerprinting was performed on 188 (81%) of 232 culture-confirmed cases for whom bacterial DNA was available. Of the 188 patients with a DNA fingerprint, 120 (64%) were infected with unique genotypes and 68 (36%) with identical DNA fingerprints were grouped into a total of 20 clusters. Based on the genotype clustering analysis, Garcia-Garcia et al. [Reference García-García18] estimated the odds ration of MDR TB compared to drug-sensitive (DS) TB was 0·16 [95% confidence interval (CI) 0·04–0·6]. Therefore, we adopted this odds ratio and used the Monte Carlo (MC) technique to obtain the best-fitted distribution of RF to capture the uncertainty.
Seasonal transmission dynamic model
The previous well-developed DR TB transmission models [Reference Blower and Gerberding13–Reference Luciani16, Reference Dye and Espinal19, Reference Rodrigues, Gomes and Rebelo20] were adopted and modified to describe parsimoniously the population dynamics of MDR TB in Taiwan. Thus a two-strain TB model was used (see Supplementary Fig. S1). Briefly, (i) susceptible individuals may be infected with either DS or MDR strains, (ii) two certain types of primary progressive TB (i.e. fast TB) and latently infected TB caused by endogenous reactivation or exogenous re-infection (i.e. slow TB) were included, (iii) a case may be spontaneously cured at a cure rate and move into the latent non-infectious state, and (iv) the MDR TB epidemic can be composed of two main subepidemics in which individuals are primarily infected/re-infected with MDR TB (i.e. primary resistance) and have DS TB treatment failure resulting in resistance (i.e. acquired resistance). The system equations of the present two-strain TB model incorporated with seasonal transmission are listed in Table 1.
Table 1. Equations for the present proposed two-strain tuberculosis (TB) model
π, Recruitment rate (person yr−1); β S(t), seasonal transmission rate for DS TB at time t (person−1 yr−1); β R(t) = RF × β S(t), seasonal transmission rate for MDR TB at time t (person−1 yr−1), where RF is the relative fitness of MDR strains; μ, background mortality rate (yr−1); p probability of new infections that develop progressive primary active TB within 1 year; c s, cure rate of active DS TB (yr−1); ν, progression rate from latency to active TB (yr−1); σ, partial immunity that decreases probability of fast progression after re-infection; c R, cure rate of active MDR TB (yr−1); μ S, DS TB-caused mortality rate (yr−1); c F, rate of treatment failure (yr−1); r proportion of DS TB treatment failure acquiring resistant; μ R, MDR TB-caused mortality rate (yr−1); N T, total population size (person).
* See Supplementary Figure S1.
This study used a regression model incorporating seasonality, temperature, relative humidity, and rainfall as potential predictors of TB to assess the characteristics of TB epidemics in Taiwan during 2005–2008. The model was fitted to data to estimate TB trends as
$$\eqalign{Y_{{\rm RM}} (t) = & \beta _{\rm 0} + \sum\limits_{n = {\rm 1}}^{\rm 5} {\left( {\beta _{{\rm 1,}n} {\rm sin}\left( {\displaystyle{{{\rm 2}n\pi {\kern 1pt} t} \over {{\rm 12}}}} \right) + \beta _{{\rm 2,}n} {\rm cos}\left( {\displaystyle{{{\rm 2}n\pi {\kern 1pt} t} \over {{\rm 12}}}} \right)} \right)} \cr & + \beta _{\rm 3} {\rm Temp} + \beta _{\rm 4} {\rm RH} + \beta _{\rm 5} {\rm Rain,}} $
where Y RM(t) is the expected monthly number of TB new cases at time t estimated based on the regression model, β 0 is the intercept, β 1,n , β 2,n , and β 3 to β 5 represent the fitted coefficients, sin(2nπt/12) and cos(2nπt/12) represent seasonality in which n is the number of seasonal cycles per year, and Temp, RH, and Rain are the mean temperature (°C), relative humidity (%), and rainfall (mm), respectively.
However, the expected monthly number of new TB cases at time t can also be calculated based on the two-strain TB model by incorporating equations (T4) and (T5) and is designated as Y TM(t):
$$\eqalign{Y_{{\rm TM}} (t) = \ & p\beta _{\rm S} (t)T_{\rm S} S + (\nu + p\sigma \beta _{\rm S} (t)T_{\rm S} )L_{\rm S} \cr & + p\beta _{\rm R} (t)T_{\rm R} S + p\sigma \beta _{\rm R} (t)T_{\rm R} L_{\rm S} + (\nu + p\sigma \beta _{\rm S} (t)T_{\rm S} \cr & + p\sigma \beta _{\rm R} (t)T_{\rm R} )L_{\rm R} + c_{\rm F} rT_{\rm S}.} $$
β S(t) can then be written as
$$ \beta _{\rm S} (t) = \displaystyle{{Y_{{\rm TM}} (t) - \nu L_{\rm S} - \nu L_{\rm R} - c_{\rm F} rT_{\rm S}} \over \eqalign {& pT_{\rm S} S + p\sigma T_{\rm S} L_{\rm S} + pRFT_{\rm R} S\cr & +p\sigma RFT_{\rm R} L_{\rm S} + p\sigma T_{\rm S} L_{\rm R} + p\sigma RFT_{\rm R} L_{\rm R}}}. $$
To estimate β S(t), we linked the two-strain TB model and regression models based on the relationship between β S(t) and expected monthly number of new TB cases. We combined equations (1) and (2) to solve β S(t).
The expressions for basic reproduction number (R 0) [Reference Rodrigues, Gomes and Rebelo20, Reference van den Driessche and Watmough21] quantifying the transmission potential of M. tuberculosis due to the subepidemics driven by DS TB (R 0S) and MDR TB (R 0R) are summarized in Table 1. R 0 is defined as the average number of successful secondary infectious cases generated by a typical primary infected case in an entirely susceptible population [Reference Anderson and May22]. When R 0 > 1 it implies that the epidemic is spreading within a population and incidence is increasing, whereas R 0 < 1 means that the disease is dying out. An average R 0 of 1 means the disease is in endemic equilibrium within the population.
Probabilistic TB risk model
To develop a probabilistic DS/MDR TB risk model, a dose–response model describing the relationships between transmission potential of DS/MDR M. tuberculosis quantified by R 0 and the total proportion of infected population must be constructed. Generally, the probability of infection for each susceptible person each day is based on the transmission probabilities for each potentially infected contact. According to Anderson & May [Reference Anderson and May22], in a homogeneous and unstructured population, the total proportion of infected population during the epidemic (I) depends only on R 0, and can be theoretically expressed as
$$I = 1 - \exp ( - R_0 I).$$
Equation (3) cannot be solved analytically. Thus, we solved equation (3) numerically by using a nonlinear regression model to best fit the profile describing the relationship between I and R 0 [equation (3)] for R 0, ranging from 1 to 10. We found that I can be expressed as a function of R 0 only,
$$I(R_0 ) = 1 - \exp (1{\cdot}63 - 1{\cdot}66\;R_0 ),\;1 \lt R_0 \lt 10,\quad (r^2 = 0{\cdot}99)$$
Equation (4) can be seen as a conditional response distribution describing the dose–response relationship between I and R 0 and can be expressed as: P(I|R 0). Thus, following Bayesian inference, DS/MDR TB infection risk (the posterior probability) can be calculated as the product of the probability distribution of R 0 (the prior probability) and the conditional response probability of proportion of the population expected to be infected, given R 0 [the likelihood P(I|R 0)]. This results in a joint probability distribution or a risk profile as
$$R(I) = P(R_0 ) \times P(I\left| {R_0} \right.),$$
where R(I) is the cumulative distribution function describing the probabilistic infection risk of a TB epidemic in a susceptible population at specific R 0, and P(R 0) is the probability density function of R 0. The exceedance risk profile can be obtained by 1 − R(I).
Model parameterization and validation
The parameter values used in the two-strain TB model [Table 1, equations (T1)–(T5)] can be estimated based on site-specific TB data obtained from Taiwan CDC, Department of Statistics, Ministry of the Interior, ROC (Taiwan) [23], and published data adopted from the literature [Reference García-García18, Reference Dye and Espinal19, Reference Dye24–Reference Yeh26]. We used the two-strain TB model to predict the seasonal incidence dynamics of TB in our four study sites for 2006–2016 with 95% credible intervals.
We validated the model performance by using the root mean squared error (RMSE) to compare the MDR TB incidence rates between predictions and observed data obtained from Taiwan CDC for 2006–2010. The RMSE is a measure of the differences between the predictive values by the two-strain TB model and observation data obtained from Taiwan CDC. The RMSE is calculated as
$${\rm RMSE} = \sqrt {\sum\nolimits_{k = 1}^K {(I_{o,n} - I_{s,n} )^2 /K}}, $$
where K denotes the number of observations, I o,n is the observed incidence rate, and I s,n is the simulation result corresponding to data point k [Reference Hyndman and Koehler27].
Uncertainty analyses and simulation scheme
TableCurve 2D package (AISN Software Inc., USA) and Statistica (version 9, Statsoft Inc., USA) were used to perform model-fitting techniques and statistical analyses. A MC technique was implemented to quantify the uncertainty and its impact on the estimation of expected risk. A MC simulation was also performed with 10 000 iterations to generate 2·5 and 97·5 percentiles as the 95% CI for all fitted models. Crystal Ball software (version 2000·2, Decisioneering Inc., USA) was employed to implement MC simulation. Model simulations were performed using Berkeley Madonna 8·0·1 (Berkeley Madonna was developed by Robert Macey and George Oster, University of California at Berkeley).
Akaike's Information Criterion (AIC) was used to measure the relative goodness of fit for the regression model. AIC is an assessment for regression model selection and can be expressed as AIC = 2k – 2ln(L), where k is the number of predictors in the regression model and L is the maximized value of the likelihood function for the estimated model [Reference Liao28]. In the regression model selection procedure, we used equation (1a ) as reference model to compare another model and then calculated the differences in AIC (∆AIC). A model was considered more likely when ∆AIC ⩾ 4. When two models had similar values with ∆AIC < 4, the model with the lowest AIC value was selected [Reference Burnham and Anderson29]. A P value <0·05 was taken as significant.
RESULTS
Optimal regression model determination
We first used a sine wave function to optimally fit the meteorological data to obtain likelihood distribution of mean temperature (r 2 = 0·91–0·96), relative humidity (r 2 = 0·18–0·62), and rainfall (r 2 = 0·32–0·45) in four study sites for 2005–2008 (see Supplementary Table S1 and Supplementary Fig. S2). We then incorporated sine wave functions of mean temperature, relative humidity, and rainfall into the regression model, Y RM(t) [equation (1a)] to capture the TB trends, allowing both the magnitude (quantitative) (Table 2 and Supplementary Fig. S2) and shape (qualitative) (Fig. 1) of trends to be estimated. We found that meteorological factors showed no significant effects on TB incidence (Table 2). The best-fitting models with ∆AIC ⩾ 4 were in Pingtung county (∆AIC = 5·92, r = 0·36) and Taipei city (∆AIC = 5·10, r = 0·52), whereas the smallest AIC values were in Hwalien (AIC = 206·28, r = 0·63) and Taitung (AIC = 182·43, r = 0·46) counties (Table 2).
Fig. 1. Comparison of monthly number of new tuberculosis (TB) cases between regression model-fitting outcomes with 95% confidence intervals (CI) and observed data for July (2005–2008) in (a) Hwalien county, (b) Taitung county, (c) Pingtung county, (d) Taipei city.
Table 2. Fitting of regression models for Hwalien, Taitung, Pingtung counties, and Taipei city, respectively, in the period 2005–2008
AIC, Akaike's Information Criterion.
* Ref. = reference model [equation (1a)].
The best fitting models are highlighted in bold font.
Seasonal population dynamics of TB/MDR TB
The results of model parameterization were listed in Table 3. We optimally fitted the published data [Reference García-García18] to obtain the likelihood distribution of RF by MC simulation. This resulted in a normal (N) distribution of RF with a mean of 0·32 and a standard deviation (s.d.) of 0·14 (Table 3). The fitted coefficients in Y RM(t) [equation (1a)] for Hwalien, Taitung, Pingtung counties, and Taipei city during 2005–2008 are listed in Supplementary Table S2 and were used in equation (2) to estimate β S(t). β R(t) can then be estimated by β S(t) × RF. We incorporated β S(t) [equation (2)], estimated probability distributions of the model parameters, and site-specific initial population size in 2006 (Table 3) into the two-strain TB model [Table 1, equations (T1)–(T5)] to project future site-specific population dynamics of TB and MDR TB incidence for 2006–2016 (Figs 2 and 3).
Fig. 2. Modelling seasonal tuberculosis (TB) incidence rates (per 100 000 population) with 95% credible intervals from 2006 to 2016 (for July) based on the two-strain TB model and the comparison of incidence data with model simulation outcomes for 2006–2008 in (a) Hwalien county, (b) Taitung county, (c) Pingtung county, (d) Taipei city.
Fig. 3. Annual incidence rates (per 100 000 population) of multidrug-resistant tuberculosis (MDR TB) estimated by the two-strain TB model varying with different percentile estimates of β R during 2006–2016 and the comparison of incidence rates between predictions and observed data for 2006–2010 in (a) Hwalien county, (b) Taitung county, (c) Pingtung county, (d) Taipei city.
Table 3. Probability distributions (N = normal, LN = lognormal) of parameter values and initial population sizes used in the two-strain TB model and basic reproduction number (R0) estimations a
a See Table 1 for explanation of symbols.
b Estimated by 0·04 (0·015−0·14) for <15-year-olds and 0·14 (0·08−0·25) for >15-year-olds [Reference Dye24].
c Estimated from Lin et al. [Reference Lin25].
d Estimated based on data from Department of Statistics, Ministry of the Interior, Taiwan [23].
e Estimated based on Taiwan CDC data.
f Estimated based on Dye & Espinal [Reference Dye and Espinal19].
h β R = β S × RF.
i RF is the relative fitness that is estimated based on Garcia-Garcia et al. [Reference García-García18].
j The initial population size in 2006 of N T, T S, and T R are adopted from Taiwan Tuberculosis Control Report (2007) [17]. S = N T – L S – L R – T S – T R. L S = 0·004 × 0·92 × 0·99 × N T and L R = 0·004 × 0·92 × 0·01 × N T, where 0·004 is adopted from Yeh et al. [Reference Yeh26], 0·92 = (1–0·08) adopted from Dye et al. [Reference Dye24], and the proportions of infections that develop L S (0·99) and L R (0·01) are assumed.
Figure 2 demonstrates the comparison of the TB incidence data with our model simulation output with 95% credible intervals, indicating that the predictions are in apparent agreement with the observed data for 2006−2008. Figure 3 shows the comparison of MDR TB incidence rates between predictions varying with different percentile estimates of β R by the two-strain TB model for 2006–2016, indicating that the predictions were consistent with observed data for 2007–2010.
The model was also extended to project the TB and MDR TB incidence rate for periods 2009–2016 and 2011–2016, respectively. The RMSE of MDR TB simulations (Fig. 3) ranging from 0·38 to 3·85 were comparable to the data's average standard deviation (s.d.) of 0·11–2·37. However, our model had the lowest RMSE values for the predictions with 50 (RMSE = 2·29), 50 (RMSE = 1·65), 2·5 (RMSE = 0·28), and 25 (RMSE = 0·10) percentiles in Hwalien, Taitung, Pingtung counties, and Taipei city, respectively, indicating that all RMSE values were less than the s.d. of observed data. Overall, the model captures the seasonal transmission and population dynamics of TB incidence for 2006–2008 and MDR TB in high TB incidence areas in Taiwan during 2006–2010.
DS/MDR TB infection risk estimates
To estimate the probabilities of DS and MDR TB infection risk, the transmission potential quantified by R 0S and R 0R had to be determined. The site-specific seasonal R 0S (Fig. 4 a, c, e, g) and R 0R (Fig. 4 b, d, f, h) were calculated based on equations (3) and (4) with input parameter values listed in Table 3. The MC simulation results showed that the optimized log-normal distribution was the most suitable fitted model for seasonal values of R 0S and R 0R (Fig. 4 i). We found that, for instance, in the highest TB epidemic area of Hwalien county, the R 0S and R 0R estimates were 0·97 (95% CI 0·26–2·03) and 0·44 (95% CI 0·06–1·45), respectively, whereas R 0S and R 0R values were 0·96 (95% CI 0·26–1·97) and 0·44 (95% CI 0·06–1·38), respectively, in Taitung county. The R 0S and R 0R estimates in Pingtung county were 0·97 (95% CI 0·26–1·92) and 0·41 (95% CI 0·05–1·25), respectively, whereas in Taipei city, the R 0S value was 0·98 (95% CI 0·27–1·77) with a lowest R 0R value of 0·37 (95% CI 0·05–1·02).
Fig. 4. Site-specific seasonal basic reproduction numbers of (a, c, e, g) drug-sensitive tuberculosis (TB) (R 0S) and (b, d, f, h) multidrug-resistant (MDR) TB (R 0R) in Hwalien, Taitung, and Pingtung counties, and Taipei city. (i) The box-and-whisker plot illustrates the overall R 0S and R 0R.
Figure 5 a presents the conditional dose–response profile showing the estimate for the total proportion of TB-infected population (I) that depended only on R 0 based on equation (4). Given the site-specific R 0S and R 0R distributions (Fig. 4 i) and conditional dose–response relationship P(I|R 0) (Fig. 5 a), the site-specific exceedance risk probability of DS and MDR TB infection can then be estimated by using equation (5) (Fig. 5 b, c). We found that DS TB in Hwalien, Taitung, Pingtung counties, and Taipei city, respectively, had nearly 17·6%, 15·5%, 14·5%, and 12·8% probabilities of the total proportion of infected population exceeding 50%, whereas there was a 29–32% probability of having exceeded 20% of the total proportion of infected population (Fig. 5 b). Our results also indicated that the selected four regions had only ∼3% probability of having exceeded 50% of the population infected attributed to MDR TB (Fig. 5 c).
Fig. 5. (a) The conditional dose–response profile representing the relationship between total proportion of tuberculosis (TB)-infected population (I) and R 0. The site-specific exceedance risks of total proportions of TB infections estimated for (b) drug-sensitive (DS) TB and (c) multidrug-resistant (MDR) TB in Hwalien, Taitung, and Pingtung counties, and Taipei city.
DISCUSSION
Seasonal dynamics of TB/MDR TB transmission
This study linked the two-strain TB and the statistical regression models to describe the seasonal transmission and population dynamics of TB and MDR TB in Taiwan. Our study showed that the seasonal peak of new TB cases generally occurred during late spring to early summer seasons in Taiwan. Thus dynamic consequences of seasonal variation in TB appeared in Taiwan. Several studies have found variable periods of peak seasonality of TB at the end of winter to late spring in Spain [Reference Luquero8], in spring to late autumn in the USA [Reference Willis10], during summer in Hong Kong [Reference Leung11], and during spring and summer in Japan [Reference Nagayama and Ohmori12].
Our results revealed that meteorological factors showed weak effects on new TB cases based on the regression model analysis. In previous analysis, we also found that there was weak relationship between (lagged) temperature and TB incidence in a Poisson regression model; however, it is interesting to note that seasonality, aborigines, gender, and age showed stronger association with TB trends in Taiwan [Reference Liao28].
Greenman et al. [Reference Greenman, Kamo and Boots30] indicated that external forcing can cause some cycles, oscillations, and even chaotic phenomenon in disease dynamics of population. This study applied five cycles of seasonality to examine periodic trends of monthly new TB cases, which was similar to those of in China [Reference Liu, Zhao and Zhou9] and USA [Reference Willis10] and the seasonality was significant in the regression model, expect in Pingtung county. Therefore, we conclude that seasonality is a major factor in shaping the TB dynamics and can also be seen as a predictor to better understanding the TB incidence patterns in Taiwan.
The causes of TB seasonality are not well understood. Some studies indicated that the transmission of TB may be intensified by increased time spent in overcrowded, poorly ventilated places, and by an increased frequency of viral infections like flu in winter [Reference Liu, Zhao and Zhou9]. Marais et al. [Reference Marais31] also indicated that viral respiratory co-infection may increase susceptibility to TB infection and progression to disease. Many works suggest that vitamin D supplementation induces antimycobacterial immunity which can decrease the risk for TB infection [Reference Martineau32, Reference Fares33]. The seasonality of TB may reflect the seasonality of human immunity to TB [Reference Willis10, Reference Nagayama and Ohmori12].
Furthermore, it is also important to consider that the seasonal peak of TB in spring may be due in part to delay in diagnosis of wintertime disease [Reference Nagayama and Ohmori12]. Missed early diagnosis of TB may occur more often in winter as coughs and fever are common seasonal symptoms, reflecting the viral respiratory infection [Reference Willis10]. The interval between observable immunological response and infection may be in excess of 7 weeks [Reference Fares33]. This could lead to a longer infectious period and increased transmission during winter, as well as more new TB cases will be detected in spring.
We used RMSE to validate the model fitness and test model predictive power in this study. The RMSE is a measure of the differences between the observed data and predictive values by optimal model. The RMSE values were comparable to the s.d. of the data. Our simulation results by the two-strain TB model showed that all RMSE values were less than the s.d. of observed data, demonstrating that the model-simulated values agreed with the observed data. Overall, our present model captures the seasonal transmission and population dynamics of total TB incidence and MDR TB in high TB incidence areas in Taiwan for the periods 2006–2008 and 2006–2010, respectively.
Given the high frequency of MDR TB in eastern Taiwan, our simulation showed that the incidence of MDR TB seems to be falling by 2013–2016, which is consistent with the downward trend of observed MDR TB incidence rate due to the effective TB control programmes. In particular, the implementation of the directly observed treatment short-course (DOTS) strategy since 2006 to prevent the occurrence of new MDR TB resulted from treatment failure and the MDR TB programme (DOTS-plus) since 2007 in providing MDR TB patients with complete and non-stop care [Reference Huang34]. Our finding also implicitly provides information that the ongoing control programmes implemented in Taiwan may succeed in curing patients with TB and MDR TB and reduce TB incidence nationwide.
Infection risk estimates of DS/MDR TB
The median values of the R 0R estimate were higher in Hwalien and Taitung counties in east Taiwan, whereas Taipei city had the lowest R 0R value. Our results also show that all R 0S were larger than R 0R in the four study sites. The persistence of both DS and DR TB (i.e. co-existence) occurs if R 0S > 1 and R 0S > R 0R. Under these conditions, co-existence can even occur when R 0R < 1 [Reference Dye3, Reference Blower and Gerberding13].
The results for infection risk assessment show that the incidences of DS TB in Hwalien, Taitung, Pingtung counties, and Taipei city had nearly 17·6%, 15·5%, 14·5%, and 12·8% probabilities, respectively, of the total proportion of infected population exceeding 50%, whereas there was only ∼3% probability of having exceeded 50% of the population infected attributable to MDR TB. Although MDR TB seems unlikely to result in an emergence, yet the case reproduction numbers of DS TB are alarming from a conservative point of view. As long as patients are carrying sensitive strains, there will always be some MDR TB cases, because MDR TB occurs from treatment failure, mutation at some constant frequency, and is transmitted occasionally for MDR strains.
Limitations and implications
The quality of the results depend on the input data. In our study, the limited datasets used in the two-strain TB model and regression model may pose the greatest sources of uncertainty. In the present study, the most unknown important parameter is RF of MDR strains. The MDR TB incidence rate of Mexico in 2008 was estimated to be 0·6/100 000 population (range 0·3–0·9/100 000) [35]. The average MDR TB incidence rate of Taiwan during 2007–2010 (0·7/100 000 population) was similar to Mexico. Thus, we estimated RF for MDR strains based on Garcia-Garcia et al. [Reference García-García18].
Several researchers used low values of RF for mathematical modelling of MDR TB such as values ranging from 0·7−1·0 [Reference Dye and Williams14] and 0·04–0·6 [Reference Dye and Espinal19]. Cohen & Murray [Reference Cohen and Murray15] assumed that the unfit MDR strains have a low fitness (0·3), whereas the fit MDR strain has RF ranging from 0·8−1·2. Burgos et al. [Reference Burgos36] investigated the relative secondary-case ratio of MDR TB to DS TB based on genotype clustering method, indicating that there were no secondary cases associated with MDR strains. Therefore, MDR strains may have lower transmissibility than DS strains, resulting in a low RF value.
Moreover, TB epidemiology values in the two-strain TB model are not easy to parameterize in Taiwan due to data limitation. To compensate for these shortcomings, we adopted published TB data [Reference Dye24, Reference Lin25] together with our study results to improve model predictability. Although the data sources were adopted from England and Wales and The Netherlands, theses parameter values could be used to examine the population dynamics of TB in different countries, such as Cambodia, China, India, Russia, South African, and Brazil [Reference Cohen and Murray15, Reference Legrand37].
In addition to seasonality and meteorological factors, we did not consider age, gender, socioeconomic factors, and subpopulation effects at the local level, where data remain scarce, and has limited MDR TB burden and time-series data. In Taiwan, the incidence rate of TB rises with age; of all new patients, 53% were aged ⩾65 years [17]. The subpopulation effects such as number of aborigines, alcoholics, smokers, and diabetics that have proven critical to TB trend variability [Reference Hsueh38]. Aborigines in east Taiwan are a subpopulation with high alcohol and cigarette consumption. Notably, the immune system is affected by age, co-infections, and past therapeutic history [Reference Keeler39]. We thus anticipated that those factors may be incorporated into our models to improve their predictability and to reflect important dynamic consequences that are worth exploring in the future.
Our integrated-level model provided an opportunity to understand the complex interactions between seasonal dynamics of TB and environmental factors and could predict the seasonal patterns of TB incidence and estimate implicitly the MDR TB infection risk based on previous data information. Our Taiwan-based analysis could be applied to assess the DS/MDR TB infection potential in high-burden countries, where M. tuberculosis remains a substantial cause of morbidity and mortality.
A better understanding of the mechanisms of TB seasonality will help to develop more effective public control programmes. The quality of the local data allows us a rare opportunity to estimate the site-specific R 0 values that could be used to establish an early warning system. We could then evaluate the site-specific probabilistic infection risk that depends only on R 0 values. The practical implications of infection risk might be initiated for risk management. We suggest that more preventive implementation in TB control should be initiated during the period of higher risk infection such as actively seeking TB cases in high-risk populations, accurately diagnosing TB and MDR TB, and preventing transmission in clinics and hospitals [Reference Dye and Williams40].
SUPPLEMENTARY MATERIAL
For supplementary material accompanying this paper visit http://dx.doi.org/10.1017/S0950268813001040.
ACKNOWLEDGEMENTS
The authors acknowledge the financial support of the National Science Council of the Republic of China under Grant NSC 100-2313-B-002-012-MY3.
DECLARATION OF INTEREST
None.
