Modelling population dynamics and seasonal movement to assess and predict the burden of melioidosis

Background Melioidosis is an infectious disease that is transmitted mainly through contact with contaminated soil or water, and exhibits marked seasonality in most settings, including Southeast Asia. In this study, we used mathematical modelling to examine the impacts of such demographic changes on melioidosis incidence, and to predict the disease burden in a developing country such as Thailand. Methodology/Principal findings A melioidosis infection model was constructed which included demographic data, diabetes mellitus (DM) prevalence, and melioidosis disease processes. The model was fitted to reported melioidosis incidence in Thailand by age, sex, and geographical area, between 2008 and 2015, using a Bayesian Markov Chain Monte Carlo (MCMC) approach. The model was then used to predict the disease burden and future trends of melioidosis incidence in Thailand. Our model predicted two-fold higher incidence rates of melioidosis compared with national surveillance data from 2015. The estimated incidence rates among males were two-fold greater than those in females. Furthermore, the melioidosis incidence rates in the Northeast region population, and among the transient population, were more than double compared to the non-Northeast region population. The highest incidence rates occurred in males aged 45–59 years old for all regions. The average incidence rate of melioidosis between 2005 and 2035 was predicted to be 11.42 to 12.78 per 100,000 population per year, with a slightly increasing trend. Overall, it was estimated that about half of all cases of melioidosis were symptomatic. In addition, the model suggested a greater susceptibility to melioidosis in diabetic compared with non-diabetic individuals. Conclusions/Significance The increasing trend of melioidosis incidence rates was significantly higher among working-age Northeast and transient populations, males aged ≥45 years old, and diabetic individuals. Targeted intervention strategies, such as health education and awareness raising initiatives, should be implemented on high-risk groups, such as those living in the Northeast region, and the seasonally transient population.


Introduction
Melioidosis is an infection caused by the Gram-negative bacillus Burkholderia pseudomallei, which exhibits marked seasonality in most settings where it is endemic, including Southeast Asia and Northern Australia [1]. Melioidosis is a communicable disease that is usually transmitted via contaminated soil or water, and is highly prevalent in Northeast Thailand [2]. Most of the population at risk of melioidosis lives in rural areas, especially those people who frequently come into contact with soil or water, such as rice farmers [3,4]. In Thailand, the highest number of melioidosis reported cases are often in January and October [5]. Infection with B. pseudomallei shows great clinical diversity, spanning asymptomatic infections, localized skin ulcers or abscesses, chronic pneumonia mimicking tuberculosis, and fulminant septic shock with abscesses in multiple internal organs [6]. Both humans and animals are susceptible to B. pseudomallei, and may be infected by percutaneous inoculation, inhalation, or ingestion. Person-to-person spread and zoonotic infections of humans are very rare [7]. The incubation period is between 1-21 days (average 9 days) [8], and is believed to be influenced by the inoculation dose, mode of infection, host risk factors, and probably differential virulence of the infecting organisms. Most cases result from recent infections, although latency with reactivation has been described up to 62 years following exposure [8], while the median times to relapse and reinfection are 21 weeks and 111 weeks, respectively. The risk of relapse is related to a patient's adherence to treatment and the initial extent of disease, but not to any underlying conditions [9][10][11]. Melioidosis seems to be more severe in older people with lower immunity or chronic underlying conditions, such as diabetes [12]. The risk of contracting melioidosis in diabetic individuals is 12 times higher than for non-diabetic individuals [13,14]. Currently, the global burden of melioidosis is estimated to be 165,000 cases per year (95% credible interval 68,000-412,000), with 89,000 deaths (36,000-227,000) [15].
Thailand's Bureau of Epidemiology (BoE) launched a melioidosis surveillance system in 2001 (Report 506) [5]. Approximately 80% of reported melioidosis cases were from Northeast Thailand [5]. In the past, the number of cases shown in the surveillance system was heavily relied on provincial and regional hospitals voluntarily report, very few were reported from private hospitals [16]. In general, melioidosis is diagnosed by testing for antibodies to B. pseudomallei using an indirect hemagglutination (IHA) technique, which has been found to have low sensitivity and specificity [17]. This surveillance system was revised in 2010 in order to capture more health data items. There has been an increase in usage of bacterial culture [16] which could give rise to an increase in total number of culture-confirmed cases. In addition, there has been an improvement to access to healthcare. Nevertheless, the true number of cases is still under-reported because of diverse clinical manifestations and inadequate bacterial identification methods. A previous estimation suggested cases in Thailand were in excess of 7,000 cases per year [15], while the BoE reported just 3,242 cases in 2015 [5].
B. pseudomallei is resistant to a wide range of antimicrobials, and ineffective treatment may result in death in 70% of cases [18]. The treatment for melioidosis consists of an intensive phase of at least 10-14 days of ceftazidime, meropenem, or imipenem, administered intravenously, followed by oral eradication therapy, usually with trimethoprim-sulfamethoxazole (TMP-SMX) for 3-6 months [19]. There is currently no vaccine against melioidosis [20,21].
The demographics of Thailand are currently in a transition phase, becoming more like those of developed countries, with rapid changes in population structure, reductions in birth and mortality rates, and a low rate of population growth. Urbanization is accelerating, and there are large annual population movements. These types of changes have been shown to have important impacts on public health and the disease burden of both non-communicable [22] and communicable diseases [23]. The population at highest risk of contracting melioidosis is the working age group. There is appreciable seasonal movement among this group as they go about their working lives. The internal migration of Thai people involves a number of distinct forms of movement within each year. Three forms have been identified in previous research [24]: a single movement, seasonal movement, and repeated movement. Seasonal migration involves people moving from the North and Northeast regions of Thailand towards the Bangkok metropolis and the Central region during the dry season (from March through to June), and in the reverse direction during the wet season (June to September) [24]. 40% of people from the Northeast are classified as seasonal migrants (a transient population) [25]. It is obvious that for person-to-person transmissible infections, there are significantly more infections when such transient individuals are considered [26][27][28][29]. However, very few studies were trying to look at the effect of transient populations on an infectious disease from a primarily environmental source which will help better describe the temporal and spatial changes of the incidence of such a disease [30]. Developed countries are also observing an emergence of melioidosis related to travelling and importation of cases [1].
To date, only a few approaches have been applied to determine the melioidosis burden, including simple maps of melioidosis [1], maps of the global distribution of B. pseudomallei, and estimates of the total incidence and mortality due to melioidosis worldwide using a statistical model [15]. Only one study has used mathematical modelling, exploring the use of childhood seroprevalence data as a marker of intensity of exposure [31]. In this study, we used mathematical modelling to predict the incidence of melioidosis in the Thai population, taking account of population changes, seasonal movement, and incidence of diabetes. The model provides multi-dimensional forecasting of melioidosis, which could be useful for targeting intervention strategies in this setting.

Demographic and seasonal movement sub-models
We generated a deterministic demographic sub-model to predict the size of the total population (see S1 Figure A). We stratified the population by age and sex into 100 annual interval classes, from 0 to 100 years old. The population in each class followed the actual population structure of Thailand between 1980 and 2015, based on birth, death, and migration rate data from the Population and Housing Census [32,33], and using the 1980 census data as the initial condition. All females in the age classes between 15 to 50 years old were considered to be capable of reproduction, with the fertility rate (fr) [34], while the death rate was age-related [35]. Members of the population were assumed to die upon reaching 100 years of age. Crude net migration rates (immigrant minus emigrant per 1,000 population) for each year had an impact on all age and sex compartments [36]. Most of the at-risk population for melioidosis lives in rural areas, especially in Northeast Thailand, so we modelled internal migration by classifying the population of Thailand into three independent groups. These were: those from the Northeast region who live at home for more than 6 months in a year (NE), the transient group or the people from the Northeast region who move seasonally between home and other parts of the country and spend less than 6 months in a year at home (T) and lastly the non-Northeast group, who live somewhere other than the Northeast (Non_NE). We created the seasonal movement sub-model to overlay with the demographic sub-model to estimate the rates of movement among them (see S1 Figure B). We solved a large set of ordinary differential equations (ODE) for the deterministic demographic sub-model and the seasonal movement submodel, defined in S1 Information on Demographic sub-model and Seasonal movement submodel, respectively.
Melioidosis case data stratified by age, sex, and geographical area were obtained from the Thai annual epidemiological surveillance reports from 2008 to 2015 [5]. Key assumptions for our model were as follows. First, the transient population data used within this model referred only to the movement of the Thai population. The movement of migrant workers from other countries could be significant but was omitted in this study for simplicity [24]. Second, diabetes progression was assumed to be irreversible, i.e. people could not move from diabetic to non-diabetic. Third, we did not consider pre-diabetes or impaired glucose tolerance. Fourth, we assumed that incidence rates of diabetes were constant over time but varied by age. Fifth, we did not focus on chronic symptoms (those of duration greater than two months), including such presentations as chronic skin infections, chronic lung nodules, or pneumonia, which only accounted for around 10% of melioidosis patients [12]. Finally, we did not focus on any behavioral factors such as excessive alcohol use.
We used R software version 3.3.3 to run and analyze the model outputs, and the deSolve package to solve the differential equations [37]. The initial parameter values were calculated from population data and disease burden. Model fitting was carried out using the Markov Chain Monte Carlo (MCMC) method, implemented with the Bayesian Tools R package as defined in S2 Information on the Bayesian framework [38]. The demographic and seasonal movement sub-models were run from 1980 (see S2 Figure A) to calibrate the model by fitting to the average migration data, including the population in the Northeast moving to non-Northeast, and the reverse direction from 2005 to 2015 [25]. We estimated seasonal movement parameters from the transient population model (see S1 Table A) and used them to run the melioidosis infection model from 2005. The model was run and fitted to the annual incidence of melioidosis by age, sex, and area by year, and seasonally by month, from 2008 to 2015 [5]. For model fitting, the DEzs method in the Bayesian Tools package allowed automatic parallelization on three cores to be used for sampling. This method allowed fewer chains to be used for estimated a large number of parameters and thus optimized the computational time [39]. Number of iterations and burn in were decided upon the model convergence by analyzing the differences between multiple Markov chains. The convergence was assessed by several measures including the standard procedure of Gelmal-Rubin [40,41] and the target acceptance rates [42]. Thirty-three parameters were estimated and the median values and credibility intervals were reported. These parameters were those representing the infection rates among both sexes in the Northeast, transient, and non-Northeast populations, (b NE a , b T a , b N a ) respectively, proportion of symptomatic cases (p E ), recovery rate from asymptomatic (σ), recovery rate from symptomatic (γ), Relative susceptibility to melioidosis among diabetic individuals when compare with non-diabetic (q), mortality/death rate for melioidosis (μ M ), amplitude (A inc ), phase angle (φ inc ) and proportion of reporting (Report) (see S1 Table A). Note that the proportion (1-Report) was defined as "Under-reporting" i.e. those symptomatic melioidosis patients that have been seen by a physician, but the physician did not report them to the public health authority for some reasons e.g. improperly diagnosis or missing report. The model was further used to predict the 20-year age-specific incidence of melioidosis among males and females in Thailand, sampling all 33 parameters from the posterior chains. The model predictions were reported as age, gender, and area-specific incidence rates over time.

Results
The demographic sub-model was able to reproduce the past population structure of Thailand from 1980 to the present (see S2 Figure A). The parameters that characterized seasonal movement were estimated by fitting the model to the population movement data (see S2 Figure B). The model showed that majority of movements were made by Northeast individuals who moved to non-Northeast areas, approximately 13,600 persons per 100,000 population per month, or 34% of all movements within a month (see S1 Table A). Moreover, the majority of movements were among those aged between 15 and 60 years old, about 19,000 persons per 100,000 population per month, or 51% of all movements within a month (see S2 Figure C).
The fitting performance is shown in Fig 2. Melioidosis cases occurred seasonally, with a peak in the wet season that lasted from May to October. The infection parameters that minimized the fit statistic, using the Bayesian method, are shown in Table 1. The highest infection rate was estimated to be 6 cases per 100,000 population per month among males aged 45-59 years old in the Northeast. The lowest rate was 0.4 cases per 100,000 population per month among females aged 15-44 years old in the non-Northeast region. Surprisingly, we found that the infection rate among the transient male population aged 15-44 years was higher than the non-Northeast population (0.8 compared with 0.08 per 100,000 persons per month). Overall 46% of melioidosis cases were symptomatic. Recovery rates for untreated symptomatic cases and asymptomatic patients were estimated by the model, with the average period of infection estimated at around 9 and 5 months, respectively. The susceptibility to melioidosis among DM population is 10.84 [95% CI 8.42-12.23] times greater than the non-DM population. If patients' treatment failed and they developed severe melioidosis, they could die within two weeks. We estimated 80% and 50% under-reporting of cases in 2008-2009 and 2010-2015, respectively.
Projections of the numbers of melioidosis cases between 2015 and 2035 are given in Fig 3. Total melioidosis incidence per year was projected to increase by 70%, from 6,569 (4,834-8,701) in 2015 to 11,173 (8,207-14,773) in 2035. The largest increase of melioidosis was projected to occur among the population aged 45-59 years old. The predicted incidence among males was two-fold greater than that of females. The majority of melioidosis cases were seen to occur in the population from the Northeast region of Thailand. The predicted incidence among non-diabetic was two-fold greater than that of diabetic population.
In Fig 4, total melioidosis incidence rates were projected to increase by approximately 10% by 2035, from 11.42 (8.5-13.4) in 2015 to 12.78 (9.6-14.9) per 100,000 population in 2035 (see Table 2). The highest incidence rates were predicted to be among those aged between 45-59 years old, followed by those age 60 years old and above. The incidence was almost double among males compared with females in both Northeast and other regions. The incidence rate among the Northeast population was more than double compared with the transient population, and almost ten times higher when compared with the other regions. This study also highlighted the importance of melioidosis among the transient population who temporally live in the risk area but had almost six times higher incidence compared with other regional populations. From diabetes prospective, the incidence of melioidosis among diabetes was predicted to be as high as 60 per 100,000 population. To summary, the risk of melioidosis among the aging population with some chronic diseases such as diabetes is presenting an increasing trend. The risk of infection among transient population, who temporary get some disease exposure during the agricultural seasons, cannot be ignored.

Discussion
Few models have been used to predict the incidence of melioidosis on either a national or global scale [14,15,43]. We applied population dynamics, seasonal movement, and the impact of diabetes to study melioidosis epidemiology in Thailand. Other approaches such as decision tree or Markov model can also be used to study melioidosis epidemiology given that the rate of transmission is constant and the system is linear. Our model fits a dynamic model for a nontransmissible disease to data on notifications only, which should allow reasonable predictions to be made as to the future course of the epidemic. However, drawing strong inferences regarding parameter values that pertain to transitions through infection/disease states after the point of infection is less safe, such that particular caution should be exercised in regards to parameters, especially those discussed here. Our findings have some similarities and differences compared with previously published work. Limmathurotsakul and colleagues used a negative binomial model to derive estimates of 7,572 (3,396-17,685) cases for global melioidosis incidence in the year 2015 [15]. This figure was similar to the incidence of melioidosis estimated in our study, which was 7,569 cases (4,834-8,701) for 2015. Both studies reached similar estimates of approximately 50% for case reporting and 40-45% for mortality/death rate for melioidosis (see S2 Table A), while the risk factors identified for melioidosis were also in agreement, i.e. being male, aged more than 44 years old, and having diabetes [13]. Two previous studies by both Thai and Australian researchers consistently showed that type 2 diabetes increased the risk of melioidosis by more than 10 times when compared with those non-diabetes [13,44], this figure is similar to our model estimates. Buckee and colleagues pointed out that seasonal disease incidence could be driven by the mobility and aggregation of human populations, which spark outbreaks and sustain transmission [30]. Northeast and transient males aged more than 45 years old were also predicted by our model to be a highest risk groups for melioidosis. Apart from reporting and mortality/death rates for melioidosis, the model also gave the estimates for some natural history of disease parameters which would be hard to  measure i.e. proportion and duration of asymptomatic infections, duration of untreated symptomatic infections, and host susceptibility to melioidosis [45]. With regard to asymptomatic infections, a few studies have tried to characterize and estimate the number of these hidden infections [44][45][46]. Our model suggested that there could be a significant proportion of asymptomatic infections (54%). Although the parameter values used in the fitting and prediction process are based on annual incidence data only, the variation in the parameter values is included to reflect the uncertainty in the predictions, and the posterior distributions represent sets of collinear parameters that reproduce the observed data. Further studies and more appropriate data are required to refine these parameters. Our model can provide many benefits for health policy planning. First, the model, in common with previous studies, estimated that only about half of all melioidosis cases were being reported. Under-reporting results in melioidosis being neglected, even more than other neglected diseases such as dengue and leptospirosis [16]. Previous study suggested that melioidosis was the third most frequent cause of death from infectious diseases in northeast Thailand, after HIV/AIDS and tuberculosis [13]. By regarding melioidosis as being less important disease has made it being further under-recognized by healthcare professionals, low health budgets to invest in intensive prevention and control, poor disease knowledge and practices among the population at risk, and finally a lack of research that would enable the development of concrete strategies to improve standards of care. Second, the model can be used to guide the design of targeted interventions i.e. predicting and identifying populations at high risk for morbidity and mortality. In line with the model's predictions, targeted intervention strategies could be concentrated among the male population of working age who live in the Northeast, as well as the transient population. These strategies could include providing health education to increase protective practices while engaging in agricultural activities, washing after work, and seeking appropriate health advice when feeling sick. To prevent deaths due to infections in older age groups, i.e. 45 years or older (65%) (see S2  Table A), national strategies could focus on early diagnosis and appropriate treatment, as well as improving diabetes screening programmes, since elderly people with diabetes may be prone to severe melioidosis. Our model has some limitations. For simplicity we assumed that diabetes influenced the likelihood of melioidosis infection alone and therefore once the person is infected with melioidosis, the diagnosis and disease progression are independent of diabetes status. Although diabetes has been shown to play roles in increasing severity and/or that medication of diabetes may also affect susceptibility and presentation of melioidosis [47,48]. We also assumed that mortality/death rates for melioidosis, incidence rates of diabetes, and seasonal movement rates were constant over time, although they varied by age. It has been suggested that mortality/death rate for melioidosis due to diabetes have decreased over time because of improved access to hospitals [49], and lifestyle changes might also have affected incidence rates [50]. In this model we classified the population into those living in the Northeast and non-Northeast, which meant that the model was unable to predict the incidence of melioidosis in locations more specifically than non-Northeast. It is important to keep in mind that melioidosis is probably prevalent in all regions of Thailand, the lack of knowledge, disease awareness and diagnosis tools led to heavily report of cases by the Northeast region only [16]. We assumed that the estimates of reporting among both males and females were the same. The annual epidemiological surveillance reports of melioidosis data used in the model included cases from all provinces around Thailand, except for those cases seen in private hospitals, which account for about 30% of hospital provision, although there is no information on the relative likelihood of melioidosis being diagnosed in different sectors. Melioidosis diagnoses reported annually by the BoE are made using an indirect hemagglutination (IHA) technique to test for antibodies to B. pseudomallei, which has been found to have low sensitivity and specificity. It could therefore potentially under-predict the number of cases.

Conclusion
Population dynamics, seasonal movement, melioidosis infection rates, and under-reporting are important components of melioidosis incidence patterns. The increases seen in melioidosis cases are partly attributable to demographic changes as working, transient, and aging population groups are more prone to develop melioidosis. The key findings from our study are firstly, the increasing trend of melioidosis incidence, especially among males aged 45-59 years old, is predicted to continue; and secondly, the male, Northeast, and transient populations aged 45-59 years old were at the highest risk of melioidosis infection. We anticipate that the modelling methods described here could be used in similar settings, especially those with reliable census data, to estimate the future melioidosis burden, as well as the potential effects of under-reporting. In addition, this modelling approach could be adapted to study other infectious diseases, behavioral changes, and seasonal movements, where demographic factors are important drivers of a population's disease burden.
Supporting information S1 File. Figure (A). Schematic representation of the deterministic demographic model.