Association between climate variables and pulmonary tuberculosis incidence in Brunei Darussalam, 2001 - 2018

Seasonality of tuberculosis is a long known but less understood phenomenon, particularly in equatorial countries. This study investigated the association between climate variables and pulmonary tuberculosis (PTB) incidence in Brunei Darussalam. Weekly data on PTB case counts and climate variables from January 2001 to December 2018 were analysed using the distributed lag non-linear model (DLNM) framework. After adjusting for long term trend and seasonality, multivariable analysis showed inverse association between PTB incidence and minimum temperature (relative risk [RR] signi�cant up till the �rst 12 lagged weeks), total rainfall (RR signi�cant after lag 12), and total sunshine hours (RR signi�cant after lag 12). We also found a positive association between PTB incidence and mean relative humidity (RR signi�cant up till the �rst 8 lagged weeks). No signi�cant results were observed for average wind speed. Our �ndings reveal evidence of seasonal variation in PTB incidence in Brunei, but with varying degrees of magnitude, direction and timing. Though explainable by environmental and social factors, further studies on the relative contribution of recent (through primary human-to-human transmission) and remote (through reactivation of latent TB) TB infection in equatorial settings is warranted.


Introduction
Seasonality of tuberculosis (TB) is long known but less understood phenomenon 1 .Previous ecological studies on TB seasonality were conducted in temperate 2,3 , subtropical 4,5 ,and tropical 6,7 countries, and investigated the relationship between TB incidence and climate variables (relating to temperature, humidity, rainfall, sunlight and wind speed).However, their ndings were not consistent, likely due to known substantial heterogeneity in disease distribution 8 .Reasons for such heterogeneity include geographic variation of disease burden, environmental factors, individual-level factors for the infectious and susceptible host and social mixing patterns in a particular setting.Time-series analyses 2,3,9 , as well as techniques incorporating spatial data and regression 7,10 , have been used to study TB seasonality.As the incubation period for TB can range from several months to 2 years 11 , it is also important to take in account the delayed and non-linear effects of climate variables when studying TB seasonality.This can be done using the distributed lag non-linear model (DLNM) framework, which was previously applied to TB 4 and other diseases such as acute respiratory infections 12 , Zika 13 , and hand, foot and mouth disease 14 .
Brunei Darussalam (pop.453,600 15 ) is a small Southeast Asian country that is located at about 4° north of the equator.The country has a tropical equatorial climate where it is generally hot and wet throughout the year, and with little seasonal variation.Most of the population live close to the coastal plain regions.
The country has also has an intermediate TB burden, with an incidence of 57 per 100,000 population in 2017 16 .
TB seasonality studies in equatorial settings are limited.This study thus aims to investigate any association between climate variables and PTB incidence in Brunei.We decided to focus on PTB due to its infectiousness and that it can be caused by human-to-human transmission.Our ndings could also be helpful to identify any patterns to inform future public health interventions, or for future modelling studies.

Data collection
Case counts of all diagnosed PTB cases in Brunei between January 2001 and December 2018 (19 years)   were compiled from the National TB Coordinating Centre (NTCC).NTCC was established as part of the National TB programme in Brunei, and has implemented TB surveillance, treatment and control programmes since 2000.All patients suspected to have any form of TB across the whole country are often referred to NTCC, or any respective district directly observed treatment, short course (DOTS) centre, for diagnosis, treatment and follow-up 17 .All modes of diagnosis for the PTB cases were included (such as smear-positive, smear-negative, and through chest X-ray and/or clinician's decision).These case counts were summed up by epidemiological week and year, based on treatment start date.In cases where the treatment start date is missing, the NTCC registration date was used.
Daily data on climate variables for the same period were obtained from a local meteorological station, located in the Brunei-Muara district where 69.7% of the country's population reside 15 .The variables provided includes total sunshine hours, total rainfall (in millimeters), average wind speed (in knots), relative humidity (RH) in percentage (minimum, mean and maximum), and temperature in degree celsius (minimum, mean and maximum).These daily data were averaged by epidemiological week and year.Any missing daily values (n= 5) were replaced with the mean value for that particular month and year.Vapour pressure (a measure of absolute humidity) was calculated using the Clausius-Clapeyron equation 18 , by inputting the mean RH values and the standard temperature and pressure conditions.

Statistical analysis
Spearman's rank correlation test was used to explore the correlation between each climate variable, and with PTB case counts.Stationarity of the time series for weekly PTB case counts and each climate variable were checked using the augmented Dickey-Fuller test.
We used distributed lag non-linear model (DLNM) framework to investigate the association between climate variables and PTB incidence.Under this model framework, negative binomial distribution was assumed to account for overdispersion, and crossbasis terms were constructed for each climate variables.These terms comprise of 2 dimensions: one specifying the conventional exposure-response relationship and the other specifying the lag-response relationship 19 .To account for long-term trends and seasonality, natural cubic splines with 7 degrees of freedom (df) per calendar year were used.Natural cubic splines with 3 df were used to describe both the lagged and non-linear effects of each climate variable.
The incubation period for TB varies considerably 11 , and there is often a delay in diagnosing TB, by about 5 months 20 .Considering as well the potential diagnosing delay of 6 months in Brunei 17 , we decided to specify lags of up to 24 weeks to capture the delayed effects of climate variables.The resulting model formula is as follows: Where E(Y t ) is the expected number of PTB cases at week t, alpha is the intercept, CB is the cross-basis function used for each climate variable to be assessed (M), and ns is the natural cubic spline function applied to account for long-term trend and seasonality.The presence of any residual auto-correlation were assessed using partial autocorrelation function plots (PACF).Any remaining autocorrelation detected was accounted for by adding lags of the model's deviance residuals into the nal model.
Although not all variables give signi cant results during univariate analysis, we decided to include 5 crossbasis terms that represent different aspects of climate variables and that are also previously known to be associated with TB incidence.The rationale here is to include these variables to control for potential confounding.The Akaike's Information Criterion (AIC) value was used to assess which variables to be included in nal model.This resulted in the choice of the following 5 variables in the nal model: average wind speed, total sunshine hours, total rainfall, mean RH and minimum temperature.To ensure minimal issues with multi-collinearity and/or correlation (due to the use of multiple crossbasis terms in a single model), consistency in results obtained between univariate and multivariate were checked using visual analysis and referring to the AIC value.
For each climate variable, we reported the cumulative relative risk (RR) and 95% con dence intervals (95% CI) at the 10th and 90th percentile, compared to their reference values (centred at their medians).Overall relationship patterns were also described using three-dimensional (3D) and lag plots.
Sensitivity analyses were conducted by repeating the analysis using natural cubic splines of 5 and 9 df for long-term trend, and also repeating the multivariate analysis with the only 4 crossbasis terms that return signi cant results in the nal model.All analyses were done in R (ver.4.1), using tseries, splines and dlnm packages 21,22 .Ethical approval was given by the Medical and Health Research and Ethics Committee, Ministry of Health, Brunei (Ref: MHREC/UBD/2019/2).

Results
A total of 3170 PTB cases were reported from January 2001 to December 2018.The median number of weekly PTB cases was 4 (IQR), ranging between 0 and 13 (Fig. 1).Table 1 shows the summary statistics of each climate variable.
[Figure 1] Figure 2 shows the overall exposure-lag-response relationship of PTB and climate variables included in the nal model.It is apparent that different variables affect PTB incidence differently, with the shape of the interaction changing with increasing lags.
[Figures 2 and 3] After adjusting for long term trend and seasonality, multivariable analysis revealed signi cant results for all climate variables, except for average wind speed.For the latter, no signi cant relationship was observed; lag plots show a relatively at pattern across most lags, with wide 95% CI at the tail end (Fig. 3, 1st row).On the contrary, we observed an inverse relationship between PTB incidence and minimum temperature, with signi cantly low cumulative relative risk (RR) at the 90th percentile when compared to the median (RR = 0.15 [95% CI: 0.03, 0.82], Table 2).This relationship was signi cant up until the rst 12 lagged weeks (Fig. 3, 2nd row).An inverted U-shaped pattern was observed for the rst 16 lags, after which an inverted S-shape emerged.
An inverse association was also observed between PTB incidence and total rainfall, with signi cantly low cumulative RR at the 90th percentile when compared to the median (RR = 0.03 [95% CI: 0.002, 0.33], Table 2).A U-shaped pattern can be observed across lags, with signi cant RR seen from lag 12 onwards (Fig. 3, 3rd row).
Conversely, low total sunshine hours was found to be associated with low PTB incidence, with signi cantly low cumulative RR at the 10th percentile when compared to the median (RR = 0.012 [95% CI: 0.02, 0.78], Table 2).This association was signi cant after 12 lagged weeks (Fig. 3, 4th row).An inverted U-shaped pattern is observed for the rst 16 lags, after which an inverted S-shape emerged.
Lastly, high mean RH was found to be associated with signi cantly high cumulative RR results at the 90th percentile when compared to the median (RR = 54.0 [95% CI: 5.67, 513.70],Table 2).This association was signi cant from lag 8 onwards (Fig. 3, 5th row).Lag plots show a at relationship for mean RH values below the median (82.7%), after which an upward trend is observed.).Though statistical signi cance was lost for the remaining 3 variables (minimum temperature in ns 5df, and total rainfall and total sunshine in ns 9df), the direction of association remained unchanged.

Discussion
In this study, we used a multivariable negative binomial model with DLNM framework, to investigate the association between climate variables and PTB incidence in Brunei.We found signi cant results for minimum temperature, total rainfall, total sunshine hours and mean relative humidity, but with varying degrees of magnitude, direction and timing, in terms of the number of lagged weeks.
Firstly, we found that high minimum temperature was associated with low PTB incidence, and that this association was signi cant until the rst 12 lagged weeks.This means that an increase in minimum temperature within the rst 3 months prior could lower PTB incidence in Brunei.On the other hand, we found low total rainfall was also associated with low PTB incidence, with signi cant results observed from lag 12 onwards.Both ndings were consistent with previous studies in China 4 and Bangladesh 23 , which both used DLNM approach but with monthly and quarterly time series data, respectively.Our study ndings were also consistent with another Chinese study that uses spatial panel modelling 10 .A plausible explanation for this nding relates to changes in human activities in both indoor and outdoor settings in response to temperature changes and rainfall, which in turn could lead to differences in the human-tohuman transmission of Mycobacterium tuberculosis 4 .Such effects can be notably observed in Ethiopia, where the rainy season negatively affects healthcare-seeking behaviour for chronic illnesses 6 .
Secondly, we also found an inverse relationship between PTB incidence and low sunshine hours, and that this association was signi cant from 12 to 24 lagged weeks.This means that decrease in total sunshine hours 4 -6 months prior could lower PTB incidence in Brunei.This nding contrasts to those found in countries with particularly temperate and subtropical climates, where TB incidence peak in spring/summer, and trough in autumn/winter 1,7 .Other studies suggest a relationship between vitamin D de ciency in winter months, a proxy for the lack of sunshine hours, with high TB noti cations 3 -6 months later 24,25 .Our nding could be due to the fact that Brunei lies near to the equator, where variations in daily sunlight exposure throughout the year are minimal.Indeed, multiple studies have shown stronger seasonality patterns in areas further from the Equator 1 .Further, although a Vietnamese study 7 reported similar results to that in temperate and subtropical countries, the authors also noted that individuals in Vietnam often take measures to prevent sun exposure; a point that was not considered in their analysis, and is also a relevant limitation to the Bruneian setting.Hence, a plausible explanation for our nding could be again related to changes in human activities at both indoor and outdoor settings.
Lastly, we found that high mean RH was positively associated with high PTB incidence, with signi cant results observed from lag 8 onwards.This means that an increase in mean RH within the last 2 months could lead to increase PTB incidence.Although the reported cumulative RR at the 90th percentile for mean RH was very high (RR = 54.0 [95% CI: 5.67, 513.70], this result should be treated with caution due to the very wide 95% CI ranges.Despite this, the positive relationship of mean RH and PTB incidence remains robust, based on results from univariate and sensitivity analyses (S1 -S4 Table ).
Our nding for mean RH could be explained by the unfavourable conditions for small droplets containing M. tuberculosis to evaporate slowly, thus allowing it to remain suspended in the air for a longer period of time 26 .Experiments also suggested that the quantity of small droplets or aerosolized particles containing M. tuberculosis produced by a TB patient is predictive of infection among household contacts 27 .However, it should be noted that this explanation could be more applicable to outdoor humidity conditions.Due to the daily hot and humid conditions in Brunei, the use of indoor air-conditioning is very common albeit for households who can afford it.Further studies that could incorporate data on airconditioning usage in households or collecting socio-economic status data, or instead focusing on indoor climate conditions would help to determine the role played by mean RH in the human-to-human transmission of M. tuberculosis.
Understanding PTB seasonality is complex due to the ways it could manifest, that is either recently through exogenous infection, or remotely through endogenous reactivation of latent TB infection.It is still undetermined whether seasonal variation affects either or both types of infection, though recent infection is likely 1 .While TB cases in high and low TB burden areas tend to be driven, respectively, by recent and remote infection 11,28 , it is unclear where areas with intermediate TB burden, such as Brunei, stands.While our study ndings and their possible explanations point more towards climate factors affecting humanto-human transmission, and thus leading to recent infections, it is still premature to rule out its contribution to remote infection.Future simulation studies on the driving factor in Brunei and/or equatorial countries could help to shed light on this.
Our study has several limitations.Firstly, our ndings cannot be translated at the individual level as this is an ecological study.Secondly, climate data was only collected from a single meteorological station.While it is plausible that weather patterns could vary across Brunei (particularly rainfall), we chose to use data from this station alone as it was the only complete set of data available locally, and the station is based in the most populated district in the country.Thirdly, we did not incorporate non-climatic factors that are also known risk factors of PTB, such as age, presence of co-morbidities and socio-economic status, which could confound the study results.The main strength of our study is the availability of a 19-year long dataset with shorter weekly intervals.Analysing using time series with shorter time intervals could possibly lead to more accurate results.Also, we opted to report the multivariate model results for 5 crossbasis terms encompassing the different aspects of climate, which would theoretically control for possible confounding.
In conclusion, our study ndings indicate that there is seasonal variation of PTB incidence in Brunei.PTB incidence is associated with minimum temperature, total sunshine hours, total rainfall and mean relative humidity; but with varying degrees of magnitude, direction and timing.These variations could be explained by environmental and social factors, mainly affecting human-to-human transmission.To better understand TB seasonality, future studies on the relative contribution of recent and remote TB infection in equatorial settings is warranted.

Declarations
Figures

Figure 1 Time
Figure 1

Table 1
Summary values for all climate variables and PTB case numbers, Brunei.
PTB and humidity-related variables (speci cally mean RH, maximum RH and vapour pressure).Moderate and signi cant association were generally observed between climate variables, with stronger association seen among related variables (such as for temperature and humidity).

Table 2
Cumulative adjusted relative risk (RR) for each climate variable at the 10th and 90th percentile, derived from the nal model.In general, univariate results for each climate variable follow a similar trend and direction of association (S1 Fig andS2Table), when compared to that from the nal multivariable model.Sensitivity analysis revealed similar trends after excluding average wind speed (S3 Table), indicating minimal multicollinearity issues.Results for total rainfall and mean RH remained robust when natural cubic splines (ns) of 5 and 9 df were used for long-term trend (S4Table and S5 Table