Malaria elimination on Hainan Island despite climate change.

Background: Rigorous assessment of the effect of malaria control strategies on local malaria dynamics is a complex but vital step in informing future strategies to eliminate malaria. However, the interactions between climate forcing, mass drug administration, mosquito control and their effects on the incidence of malaria remain unclear.
Methods: Here, we analyze the effects of interventions on the transmission dynamics of malaria (Plasmodium vivax and Plasmodium falciparum) on Hainan Island, China, controlling for environmental factors. Mathematical models were fitted to epidemiological data, including confirmed cases and population-wide blood examinations, collected between 1995 and 2010, a period when malaria control interventions were rolled out with positive outcomes.
Results: Prior to the massive scale-up of interventions, malaria incidence shows both interannual variability and seasonality, as well as a strong correlation with climatic patterns linked to the El Nino Southern Oscillation. Based on our mechanistic model, we find that the reduction in malaria is likely due to the large scale rollout of insecticide-treated bed nets, which reduce the infections of P. vivax and P. falciparum malaria by 93.4% and 35.5%, respectively. Mass drug administration has a greater contribution in the control of P. falciparum (54.9%) than P. vivax (5.3%). In a comparison of interventions, indoor residual spraying makes a relatively minor contribution to malaria control (1.3%-9.6%). Conclusions: Although malaria transmission on Hainan Island has been exacerbated by El Nino Southern Oscillation, control methods have eliminated both P. falciparum and P. vivax malaria from this part of China.

Introduction

Malaria is a mosquito-borne infectious disease caused by single-celled microorganisms of the Plasmodium group and spread by Anopheles mosquitoes. Although global malaria incidence and mortality rates have decreased since 2000, it remains one of the world’s most serious public health concerns, with 228 million new cases reported in 2018[1], particularly in lower and middle-income countries. Although the World Health Organization (WHO) has set the goal of reducing malaria burden by 90% by 2030[2], a recent study has suggested that the complete eradication of malaria by 2050 may be feasible if the appropriate strategies and actions are implemented[3]. For endemic areas, the main priority is to control and ultimately eliminate malaria[4-6]. Considering many Southeast Asian countries are still trying to eliminate this disease by 2030 or 2040, there’s an urgent need to conduct studies using real-world data of malaria interventions in the west-pacific region.

A few countries have been successful in eliminating malaria using a combination of control strategies[4,7] and those intervention methods have also been well documented[8-10]. In addition, changes in environmental conditions[11-18] can drive malaria dynamics by influencing the mosquito and parasite life cycles[19-29]. Rainfall is required to establish suitable mosquito habitats, adequate levels of humidity enable high activity and survival of mosquitoes, and temperature affects multiple stages of mosquito and parasite development, as well as biting rates[30-32]. However, assessment of the effects of interventions on local malaria dynamics is complex[33-35], and an evaluation of their effectiveness at reducing incidence is difficult to ascertain due to the presence of multiple factors[36,37]. Little is known about the effectiveness of interventions at reducing the malaria burden when accounting for the potential effects of climate changes[34,38].

The epidemiological and parasitological surveys conducted on Hainan Island between 1995 and 2010 provide a longitudinal and comprehensive dataset for the assessment of the relative impacts of interventions, including mass antimalarial drug administration (MDA), indoor residual spraying (IRS), and insecticide-treated bed nets (ITNs), as well as climatic factors on driving malaria epidemics. Specifically, we quantify multiple exposures, nonlinear feedbacks, and complex interactions between interventions and transmission dynamics, which is an essential component of successful malaria control. A number of factors make the area well suited for analysis in this study; it has undergone both passive and active surveillance of cases, and our entomological survey indicates that mosquito population abundance has remained approximately constant over time[39,40]. In this study, we assess the effectiveness of control methods using mathematical and statistical methods. We find that ITNs are the most effective strategy in controlling malaria transmission in Hainan. It is also shown that MDA has a greater contribution in the control of P. falciparum.

Methods

Description of study sites

The study area Qiongzhong is located in the highland area of Hainan Island, China, with a population of more than 200,000. The county covers 2704.66 km2, and has a population density of over 70/km2. Its altitude ranges between 1200 m above sea level in the south to 1400 m in the north and to 1800 m in the west. The main epidemic season occurs from May to October during the rainy season.

Epidemiological data

A time series of monthly cases of Plasmodium falciparum and Plasmodium vivax infections were obtained from the Hainan Provincial Centre for Disease Control and Prevention. Each case was confirmed through a microscopic examination of blood slides from clinical (febrile) patients seeking a diagnosis and treatment, according to the diagnostic criteria of the Chinese Ministry of Health[39]. Experimental procedures were performed in compliance with guidelines established by the Chinese Centre For Disease Control and Prevention and have been approved by ethics committee of Hainan Centre for Disease Control and Prevention. As the research involved no risk to the subjects and it used data from an anonymized dataset, informed consent was not required in this study.

Between 1995 and 2010, mass blood examinations were conducted annually, comprising 293,117 blood smears. Antimalarial drug campaigns were conducted biannually in both previously infected and uninfected people in April and August, together with vector control (indoor residual spraying, IRS) and prevention (insecticide-treated bed nets, ITNs). A total of 122,510 people took antimalarial medication, with a strategy that consume 8 days of piperaquine with primaquine. An area of 4,231,517 m2 was sprayed on indoor areas and 170,129 insecticide-treated nets were used, with the rate of bed net using ranging from 0.37 to 7.88 per 100 people each year. All these statistics were gathered in Qiongzhong highland. The time series for the population sizes of inhabitants was obtained from the Hainan Statistical Yearbook.

Climate data

Climate data, including temperature and rainfall, were obtained from the local meteorological station in Qiongzhong from 1995 to 2010. The Niño 3.4 index[41] is calculated as the difference between the monthly average sea surface temperatures in the region 5°N–5°S, 120°W–170°W. Only Niño 3.4 index was included as a climatic component in both statistical model and mathematical model.

Statistical model

The associations between climate, antimalarial drug campaigns, IRS, and ITNs with malaria incidence were assessed with a statistical model, Generalized Additive Model (GAM). We fitted separate negative binomial regression models to the monthly incidence of P. vivax and P. falciparum.where Y is the monthly malaria incidence, t is the time, α is the intercept, and β, γ, φ and λ are regression coefficients. fSEAS is a nonlinear smoothed functions of seasonality. CLIM is Niño 3.4 index. IRS is the area of IRS and ITNs is the number of insecticide-treated bed nets. The statistical model was used for exploring initial relationships between climate, intervention factors, and malaria incidence. Only significant variables were included in the epidemic model.

Epidemic modelling

To investigate the impact of climate condition and specific control intervention on malaria mitigation, we established a mechanistic model[20,42,43] by fitting to the number of new confirmed cases reported each month using Bayesian Markov Chain Monte Carlo methods[44]. Comparing to these previous models, we included new parameters representing for ITNs and IRS in mosquito classes. To estimate the impact of MDA, the class T (antimalarial drug treatment with temporary immunity) has also been included. We used the fitted model, with posterior estimates of parameters (SI Appendix, Fig. S7, Table S2), to simulate malaria epidemic, with and without control measures. The model is:

P. falciparum model

P. vivax modelwhere human classes consist of S1 (susceptible), E (exposed), I (infected), S2 (recovered subpatent status, with partial immunity to reinfection), H (dormant liver stage), and T (antimalarial drug treatment with temporary immunity). Multiple H classes were used to partition the dormant stage into a sequence of identical stages (i = 2, …, n)[42]. T and T represent preventative antimalarial drug treatments administered to the general population and patients with previous malarial infections before the epidemic season in April and August every year, respectively. S1 + E + I + H + S2 + T = N, and b and d are birth and death rates of the population, both of which were obtained from statistical yearbook. λ is the force of infection at the current time. μ represents transition rate between classes. t is the treatment success probability. s represents superinfection from S2 to I and r is the rate at which the infected population transits into dormant liver stage. The observed monthly cases equal the number of new cases predicted by the model multiplied by the reporting rate r. The annual number of reported cases is strongly correlated with positive rate of blood examinations in our study area, a constant reporting rate over time is then assumed.

The role of mosquitoes in transmission is represented through a delayed equation between the latent f and current λ force of infection, taking into account the extrinsic incubation period τ[20].where IRS (t) = θIRS × area of IRS (m2/per person) (t), ITNs (t) = θITN × coverage of insecticide-treated nets (t), CLIM (t) = θCLIM × Niño 3.4. Seasonality, βseas, is modelled nonparametrically through the coefficients β of the periodic cubic B-spline basis φ (t). Δ is a vector of dummy variables of length 12. We developed the model shown in Fig. S1. We use k = 2 in the model.

We fit the malaria model using a Bayesian state space framework. Model fitting was performed using Metropolis–Hastings Markov chain Monte Carlo (MCMC) algorithm with the MATLAB (version R2016b) toolbox DRAM (Delayed Rejection Adaptive Metropolis). In model parameterization, we chose a Gaussian prior for the distribution of parameters, with a mean value of 0 and a variance of 102. After a burn-in of 1 million iterations, we ran the chain for 10 million iterations sampled at every 1000th step to avoid auto-correlation (SI Appendix, Fig. S7). We also fitted our model and estimated the parameters using 1995–2004 data, and then predicted malaria cases in 2005–2010. We found that our model is able to predict and capture the trend of malaria, which reflects the effect of interventions on malaria transmission (SI Appendix, Fig. S8).

Wavelet

Wavelet analyses were performed to investigate periodicity of ecological time series[45-47]. Wavelet analysis estimates the spectral characteristics of time series as function of time[45,46]. This method with the wavelet power spectrum (WPS), allows us to determine the contribution of variance at different times and periods. Here for the wavelet decomposition we adopted Morlet wavelet to analyze the P. vivax and P. falciparum time series. Additionally, a bootstrapping method was used to evaluate the statistical significance of the WPS (more details are found in ref. [47]).

Table 1

Parameter estimates included in the negative binomial regression model.

Covariates	Coefficient (standard error)
	P. vivax	P. falciparum
Intercept	−3.76 (2.27)	−5.10 (6.71)
Seasonality	EDFa 3.75**	EFD 2.52
CLIM	0.17 (0.08)*	0.13 (0.24)
MDA	−0.09 (0.11)	0.09 (0.29)
IRS	0.11 (0.22)	−0.10 (0.64)
ITNs	−0.30 (0.04)**	−0.29 (0.12)*

CLIM climatic component, MDA mass antimalarial drug administration, IRS indoor residual spraying, ITNs insecticide-treated bed nets.

**P < 0.01, *P < 0.05.

aEFD effective degrees of freedom.

Table 2

Comparison of the dynamic models including the seasonal component, climatic component, vector control, and prevention of transmission of observed malaria epidemics.

Model	P. vivax		P. falciparum
	R2	*DIC	R2	DIC
Seas	0.25	897.64	0.36	117.45
Seas + CLIM	0.34	893.27	0.47	117.72
Seas + CLIM + MDA	0.34	882.90	0.48	113.33
Seas + CLIM + MDA + IRS	0.33	832.35	0.48	114.88
Seas + CLIM + MDA + IRS + ITNs	0.68	500.14	0.49	109.44

Seas seasonal component, CLIM climatic component, MDA mass antimalarial drug administration, IRS indoor residual spraying, ITNs insecticide-treated bed nets. *DIC deviance information criterion.