A spatial predictive model for malaria resurgence in central Greece integrating entomological, environmental and social data

Malaria constitutes an important cause of human mortality. After 2009 Greece experienced a resurgence of malaria. Here, we develop a model-based framework that integrates entomological, geographical, social and environmental evidence in order to guide the mosquito control efforts and apply this framework to data from an entomological survey study conducted in Central Greece. Our results indicate that malaria transmission risk in Greece is potentially substantial. In addition, specific districts such as seaside, lakeside and rice field regions appear to represent potential malaria hotspots in Central Greece. We found that appropriate maps depicting the basic reproduction number, R0, are useful tools for informing policy makers on the risk of malaria resurgence and can serve as a guide to inform recommendations regarding control measures.


Introduction
Malaria is one of the well-studied vector-borne diseases in terms of its transmission and the potential for transmission's change. Malaria is endemic in about 100 countries around the world, mainly in sub-Saharan Africa and Asia ( [1]). It is a mosquito-borne parasitic infectious disease, transmitted through the bite of the infected female Anopheles mosquito. Five types of plasmodia cause disease to humans: Plasmodium falciparum, Plasmodium vivax, Plasmodium ovale, Plasmodium malariae and Plasmodium knowlesi, of which P. falciparum and P. vivax are the most prevalent and P. falciparum the most dangerous ( [2]). In 2015 the number of malaria deaths globally were 438000 whereas 3.2 billion people being at risk, with more deaths occurring in Africa (90%), followed by South-East Asia (7%) and the Eastern Mediterranean region (2%) ( [1]). (Climatic determinants are considered particularly important for malaria, since both the disease agent (Plasmodium) and vectors (Anopheles mosquitoes) are strongly affected by climate ( [3]). Global warming has been considered as a potential risk for malaria resurgence in northern hemisphere areas ( [4] between malaria transmission or vector occurrence and climate in order to project the potential future distribution of malaria transmission areas ( [5]). Epidemics in general, can be caused, except from climate conditions, by many other factors, for instance movement and displacement of human populations. Recently, an increasing number of studies point towards a shift from climate as the main driver for malaria occurrence, to additional non-climate factors as responsible for malaria occurrence and transmission ( [6]; [7]; [8]). Others suggested a more realistic approach, proposing that this complex system of malaria transmission is based on a combination of climate and non-climate factors ( [9]). A notable such factor relates to movement and displacement of human populations. Human migrants infected elsewhere can move, temporarily or permanently, into a new non-infected region or country. There is a significant body of literature that concentrates on mobile populations and their importance as a malaria transmission factor (see, e.g., [10]; [11]; [12]). Understanding the spatio-temporal distribution of risk for mosquito-borne infections, such as malaria, is important ( [13]).
Free movement within the European Union has been largely facilitated by an environment essentially free of important diseases. However several European regions have recently seen an influx of new-coming populations, including migrants and refugees. This fact, combined with natural habitats which can serve as hot spots for malaria resurgence due to their suitability in terms of mosquito populations, suggest that increased vigilance is required in order to avoid the (re-)establishment of previously endemic diseases ( [14]). Already in the years 2011 and 2012, in Southeast Europe renewed malaria transmission was observed, with Greece and Turkey counting localized outbreaks as a result of malaria importation from endemic countries ( [1]).
Malaria has been eradicated from Greece in 1974 after coordinated efforts by the WHO and the local authorities ( [15]). Over the last decade a large population of migrants from countries were malaria is endemic came to Greece. As a result, in the recent years, a small number of imported cases of malaria have been reported annually ( [16]). Although these cases are sporadic, they raise the chance of local transmission ( [17]). Specifically, in 2009 and 2010 some cases of P. vivax malaria (51 and 44, respectively) were laboratory confirmed in Greece (mainly in the agricultural area of Evrotas, Laconia in Peloponnese in Southern Greece). In 2011, an outbreak reached a total number of 96 laboratory confirmed P. vivax cases of which the 54 were imported cases of migrants from malaria endemic countries, and 42 were domestically transmitted cases ( [16]). In 2012, a total of 93 laboratory confirmed cases of malaria were reported in Greece, of which 73 were imported (64 in migrants from malaria endemic countries). After the 2011 outbreak, a multidisciplinary strategy with a variety of intensive response activities, was adopted and implemented in Evrotas ( [15]). Despite the fact that during 2014 Greece had no domestically transmitted cases of malaria, there were 6 new cases in 2015, indicating the need for constant vigilance in the region ( [1]). See Table 1 below for an analytical description of malaria cases during 2009-2015.
From the information presented above, we may recognize that malaria represents an important emerging public health issue in Greece, particularly in the Central and Southern regions of mainland Greece ( [15]; [18]). A malaria resurgence risk in Greece varies dramatically with space and depends, among other factors, on environmental changes as well as the introduction of parasite carriers. Therefore, it is imperative that malaria early-warning systems are put in place to enhance public health decision making for prevention and control of malaria epidemics. This aim can be achieved via the development of a prediction model which incorporates all the relevant sources of evidence, including mosquito abundance and potential disease transmitters.
In the current paper, in order to determine the potential transmission risk of malaria resurgence in Greece, we develop and present an appropriate host-vector spatially explicit model which integrates entomological, geographical, social and environmental evidence in order to guide the mosquito control efforts. This model, which is a suitably tailored version of the classical Ross-Macdonald model ( [19]) combines entomological, environmental and social data. We present spatial maps with three distinct but complementary severity scales; (i) the basic reproduction number, which is a combination of human to mosquito and mosquito to human transmission and represents the (mosquito-driven) disease potential from one (infected) human host to other humans, (ii) the probability of getting infected, which incorporates the initial number of infected humans and (iii) the expected number of infected cases.

GIS database and data collection
Study area. The entomological evidence for this study is informed by a large-scale mosquito abatement program conducted by Bioapplications Ltd during the last three years in central Greece. The study took place in 8 municipalities of the Prefecture of central Greece (see S1 Fig). The study area covers about 406.000 ha.
We established "MALGDB" as a Malaria Geographical Database which gathers together geographical, environmental, population, entomological and epidemiological data. We used ArcGIS 10.3 software to create, develop and populate that database. Nomad rugged handheld computers running Arcpad 7.1 software were used to gather field data and information. All data collected were sent via GPRS connection to the main MALGDB. The database that was build contains administrative units, populated places, digital elevation model, land cover, surface temperature and humidity, population density and composition, # of migrants from malaria endemic countries, larvae and adult anopheles positions and population, human cases, etc.
As concerns the species involved in the study, it has been found after morphological examination that the species Anopheles sacharovi, An. maculipennis, An. superpictus, An. claviger and An. hyrcanus are the most important species in the study region (account for approximately 90% of all Anopheles species in the region). It should be noted that all these types of identified species are potential malaria vectors. According to the above, in our study we suppose that all the positive samplings for Anopheles mosquito larvae are larvae belonging to species that can transmit malaria.
As regards the larval sampling, net mesh has been utilized. The diameter of the net mesh was 20 cm and for each sample we scanned 5 different lines of 1 meter long each. Based on the morphology of the natural (relief, streams etc.) and the human environment (lakes, rice fields, irrigation system, settlements) we selected 10 areas where we had established sample-collecting stations. The sampling stations were situated in the outskirts of towns and villages. Each station consisted of CO 2 mosquito traps. Breeding sites included canals, rice pads, tanks etc. Mosquito breeding sites were monitored on a weekly basis from the end of April to the end of autumn in order to record the larval dynamics and some of the physicochemical parameters of breeding sites, such as pH temperature conductivity. In total, there were more than 4000 different mosquito breeding sites in the study areas. Up to 200 different sites were positive at least once for Anopheles larvae. Following collection, mosquito samples were transported to the laboratory, where they were identified and separated by species using morphological identification. The samples provided evidence regarding the variation of the population of Anopheles mosquitoes between the different seasons of the year and highlighted the areas where Anopheles occurred in high density.
The main malaria reservoirs in Greece relates to migrants from countries where malaria is endemic. Thus, the exact locations of these migrants' habitats were marked up to a GIS database by using field computers. We also kept records of their population in each site. This migrant monitoring was repeated in 3 different periods (in June, early August and mid-September respectively) in order to obtain as possible as accurate data.
The larvae sampling and monitoring was conducted by using digital sampling protocols. We created a user-friendly database (domains in geodatabases), for data collected from the field and loaded in field rugged computers (see Fig 1a and 1b). In total, the larvae data from the study area comprised of more than 4000 records. Every 10 to 15 days we were placing CO 2 traps in 10 places. The results from the adult trapping were embodied in the MALGDB. Temperature data were recorded from the meteorological stations of National Observatory of Athens (NOA) (see Fig 2a and 2b for examples of temperature collection data at varying time periods). There were 10 meteorological stations of the NOA recording temperature data at the study area. For our modeling, we utilized the average temperature values with a 30-day interval, starting from the first Saturday of each month, between the April and November of each year (2012 & 2013). The inverse distance weighting (IDW) interpolation method was additionally used to estimate the average temperatures at the locations where no measurements were available.

Methods
The main disease severity measure utilized in this paper is the basic reproduction number, R 0 , which can be interpreted as the expected number of (human) hosts that would be infected by the introduction of an infected host into a large fully susceptible population. The prospects for the success of disease control depend crucially on this measure ( [20]). Hence, R 0 appears to be an appropriate measure and our control strategies are mostly concerned with proportionally reducing R 0 . We could also calculate the vectorial capacity, a purely entomological equivalent of R 0 , defined as the expected number of mosquito bites that would eventually arise from all the mosquitoes which would bite an infectious individual on a single day. However, since the mobile populations appear to be an important potential driver of transmission we mostly focus upon R 0 since it incorporates the vectorial capacity as well as the exposure due to potentially infected human hosts within a single measure. Some of the parameters necessary for R 0 estimation were simulated based upon values reported in the literature and these results were fused with entomological information collected in the field by Bioapplications Staff (see S1

The spatial predictive model
In this paper we utilize a host-vector model that combines entomological, environmental and geographical data to provide estimates on the average infection number due to malaria in central Greece. To achieve this, the methodology proposed here takes into account the potential host population in the region, related to migrants from countries where malaria is endemic, in addition to the standard entomological parameters which are based on the well- we seek to estimate the probability of infection, τ = Pr(infection). If we denote witht this estimate, let us assume that the latter is a function of two different measures, one describing the disease potential due to mosquito abundance and the other the component attributed to host infections. The reproduction number R 0 is utilized for the former whereas for the latter transmission route we use the proportion of initially infected human hosts, denoted by μ 0 .
For the calculations concerningR 0 , which describes the expected number of hosts that would be infected by a single infected host in a large susceptible population, we are obtaining a different estimate of R 0 for each region, calculated as: where V denotes the vectorial capacity, i.e. the expected number of infective mosquito bites that would eventually arise from all the mosquitoes that would bite a single fully infectious person on a single day ( [22]), and is given by: In Eqs 1 and 2, m i denotes the ratio of mosquitoes to humans in each region i; α i the biting rate, i.e. the proportion of mosquitoes that feed on humans each day; b i the probability a bite produces infection to a human; r i the average recovery rate per day; v i the mosquito latent period, i.e. the number of days from infection to infectiousness; g i the mosquito instantaneous death rate per day. Finally, with c we denote the probability a mosquito becomes infected after biting an infected human, which for our analysis is set to the constant value 0.5.
The parameters α i , b i , r i and v i were sampled from suitable distributions according to the relevant literature (see S1 Table for a detailed description of model parameters and corresponding values and distributions assigned).
For the purposes of the current study, we modify Eq 2 to account for the more realistic assumptions of temperature-dependence of the mosquito latent period (v i ) and the mosquito instantaneous death rate per day (g i ). Hence, since v i is known to be strongly dependent on temperature, T i , (see [23]; [24]; [25]), we incorporate the latter dependence into our modeling via the following equation: where the pre-blood meal period (or infection delay), δ i , is a function of temperature through: The above relationship has been obtained through experiments on Anopheles mosquitoes ( [24]). Note also that in order to calculate g i we are assuming that it changes with the temperature levels ( [13]), since that temperature is known to influence the mosquito life cycle and in particular the development rate of larvae and adult survival ([26]), hence: where each g i corresponds to the maximum of the calculated values for each one of the different temperatures, i.e. g i = max {g ij }, for temperatures T i recorded at i = 1, 2, ‥, 30 (see also [27] for a study on how to modify the mosquito daily mortality rate g i to account for the immature mosquito stages).
As regards the estimation of the external host component due to the migration, denoted by μ 0 , expressing the proportion of initially infected hosts in each one of the 710 and 183 subregions for 2012 and 2013 respectively, we use the following approximation: where W ik is a simple exponential kernel function of the form: that is used to model the spatial component of migrant transmission, through the distances d ik from the larvae areas, measured during the three periods of migrant monitoring k = 1, 2, 3. The prevalence of asymptomatic malaria varies widely between and within country (see e.g. [28]; [29]; [30]) and also depends upon the detection method. Therefore, we used three distinct scenarios where the baseline prevalence was set to 10% and we also examined a prevalence of 5% and 20% in deterministic sensitivity analyses. Hence, the estimated proportion of initially infected hosts,m 0i , is multiplied with the latter pre-specified incidence rates.
Having obtained theR 0i andm 0i estimates for each region, we proceed to the estimation oft by solving the following non-linear equation forR 0i ! 1: 1 þm 0i Àt i À expðÀt i ÁR 0i Þ ¼ 0; whereas forR 0i < 1 we sett i ¼ 0: Accordingly, we can estimate the E(infections). We have used the R package ( [31]) for data manipulation and in order to fit the above described model to central Greece malaria resurgence data for the years 2012 & 2013. The R code of the spatiotemporal model is available upon request by the corresponding author. The data for the years 2012 & 2013 can be obtained from S1 and S2 Files, respectively.
In order to aid the illustration of our results, we construct R 0 maps that can be used for informing policy makers on the risks of malaria resurgence. Similar maps, in the case of other infectious viruses, have been presented elsewhere; see for instance [32] for an application of An altitude zone that includes the areas between 0 and 300 m was applied as a threshold in order to avoid the depiction of the density in areas which are unlikely to be affected from the represented phenomenon. Therefore, we decided to effectively preclude interpolation methods since those may or may not be appropriate in our setting due to the presence of intermediate areas of high altitude which could have been depicted as high risk when, in fact, this is not the case.
The method used to map the results was "Kernel Density" (Spatial Analyst of ArcGIS-Arc-Map v.10.3). The Density tool distributes a measured quantity of an input point layer throughout the landscape to produce a continuous surface. These geographical maps can be used to identify the areas of higher risk for a malaria outbreak after the introduction (see a detailed explanation of the method in S3 File). As we observe, for both years there exists a geographical similarity as regards the potential malaria transmission hotspots.

Results
Approximately ten distinct geographical units in malaria resurgence can be identified by inspecting the two graphs. Specifically, the areas of higher resurgence risk can be identified around the cities of Lamia, Levadeia and Chalkida, as well as in the Northern and Southern parts of Evoia. Generally, risk of malaria resurgence is increased in seaside areas, areas near lakes where the Anopheles mosquito population density peaks (e.g. Chalkida), or at the Anopheles breeding habitats such as areas with paddies (e.g. Lamia).
The findings point to the fact that resurgence is possible in both areas of sparse human populations (rural areas) and areas with dense human populations (urban areas). In case of disease resurgence or outbreak, interventions could be targeted in these areas identified by our analysis as regions of high transmission intensity.
Another important-complementary to the R 0 -factor in disease modeling is the probability of infection, τ, and accordingly the amount of people that are expected to be infected in case of the disease resurgence in the areas under investigation. These two outcomes are complementary but add extra source of information, hence we additionally present the results of the two outcomes in the supporting information. S2 to S5 Figs show the latter measures, for the years 2012 and 2013, respectively. We can see from the inspection of thet estimates (S2 & S3 Figs) that the human populations with the higher risk for infection are located mainly in the wider region of the city of Lamia, and in the region of Northern Evoia. This is a result that seems to be time invariant for both years of research. In addition, a high risk for infection is Next, as regards the model estimates of the numbers of expected infections from malaria (see S4 & S5 Figs), we observe similar results between the years 2012 and 2013. The highest amounts of people expected to be infected are widely scattered in the regions of central Greece, with higher at the city of Lamia and neighboring regions. We also observe high values for northern Evoia and other seaside areas.
As expected, malaria transmission is highly seasonal, with transmission limited to the warm and rainy summer months. Cases in general increase by the end of June, peak in late summer and decline by the end of November. The average seasonal pattern in malaria incidence follows the periodicity in rainfall and temperature, with a 3 to 4 month lag ( [3]). In the following table (Table 2) the summary values of median R 0 are presented, for various periods of the two years as were estimated from our model, thus giving a picture of the seasonality pattern. Besides the total median R 0 estimates, we also report the median R 0 for both the urban and rural regions of Central Greece. As we observe, there is a seasonal trend of R 0 increase in both years, occurring during the summer months when the temperatures are obviously higher. The R 0 values are showing a decline by October-November, where the transmission potential is lower. We note that the major peaks in potential risk resurgence are in July, August and September.
It is also noteworthy that malaria resurgence is more likely to occur in rural regions of Cen- Finally, S4 Table shows the results of sensitivity analysis where we are varying the value of the baseline prevalence. As seen in S4 Table, the three outcomes do not alter significantly from the changes in the baseline prevalence.

Discussion and conclusions
We adopted a model-based framework that integrates entomological, geographical, social and environmental evidence in order to examine the potential for malaria resurgence in Greece. In addition, disease risk maps were generated in order to assist the interpretation of the results. The results clearly indicate a potential malaria transmission risk in Greece. Our results characterize the higher and lower risk areas of malaria resurgence in parts of central Greece. We term as hotspots for malaria resurgence those areas with R 0 > 1. It has been shown that there are spatial variations in the risk for malaria resurgence. Similar results have been reported in another recent study modeling the risk of malaria resurgence using a mathematical model in Portugal ( [10]). Malaria resurgence in Greece may occur, even in areas that are currently free from the disease. Among the main findings is that at elevated temperatures the new-emergence regions are more susceptible to potential outbreaks. Therefore, authorities should be vigilant in order to avoid future outbreaks in areas of high disease potential. We also found that the potential for malaria transmission resurgence, especially in the regions identified by our study, is affected not only by the potential of the virus or the specific climate conditions but also by the number of potential disease hosts who are the main drivers of disease risk. Our model highlighted areas of central Greece as being particularly suitable for disease resurgence. One important question relates to the stability of potential regions for disease resurgence. Our study is indicative of a positive answer to this question, since that the comparative analysis between the two years under study did not showed significant variations. The spatial clusters of high probability for disease resurgence generally do not seem to vary significantly over time. Further research is warranted in this direction since a large time-horizon will give stronger evidence regarding this issue.
The results emphasize the importance of maintaining population and health systems awareness on the potential resurgence of malaria in Greece. Enhanced prevention and control strategies should be planned for rapid implementation in the case of future pathogen importation to prevent malaria resurgence. An approach of this kind ought to be used in order to plan, control, assess and manage a malaria prevention and control program in Greece. In particular, this study suggests a prediction mechanism which should be put in place so that mosquito control programs be efficient in the case of malaria resurgence along with rapid public health response to treat any infection. These two types of action operate in a complementary manner. Therefore, continuous data update and new data monitoring is critical in order to estimate the relative merit of each and assist public health authorities to respond to potential risks of malaria in a cost-effective manner.
One immediately apparent outcome of this work relates to the need to strengthen surveillance and perform integrated mosquito control programs that will help to eliminate the potential risk of malaria reintroduction and reestablishment. We used a contemporary technique which offers control over a massive amount of field data in real time. Consequently, we maximize the effectiveness of the abatement programs in terms of public health risk. This approach suggests a mechanism for efficient and reliable mosquito control programs: use GIS technology and telematics in order to plan, control, assess and modify abatement applications in-situ and in a real-time manner.
Our study has a number of drawbacks. The selected regions were based upon field studies and were not selected within a large randomized experiment. Therefore these results are representative of the studied (and similar) region but further studies are required in order to generalize these findings to wider regions. Also, as with all such field studies, we obtained approximate larvae numbers and a rough estimate of their survival based upon its relationship with temperature. In addition, weekly-and not daily-samples were collected so these numbers are subject to sampling bias. Future studies will mark selected mosquitoes and use capture-recapture methodology for more accurate vector estimates. Also, while validation against an external/independent source of evidence is desirable whenever an analysis with potential for policy action is performed, this seemed nearly impossible in our case since we are essentially looking at counterfactual scenarios because what actually happened is the result of different interventions of distinct nature for each area. However, when interest focuses on highlighting the disease potential and the lack of evidence regarding several important aspects of the disease, this can be thought of an acceptable drawback. Finally, we had limited data on the mobility patterns of potential hosts. Further evidence regarding their movement habits would inform a more realistic dynamic movement network leading to more powerful results. Such information, if available, can easily be fused within our model-based framework.
The current study represents an attempt towards an integrated method for predicting risk resurgence of malaria, based on a spatio-temporal mathematical modeling approach, combined with a proposed framework for constructing R 0 maps. Future work should include larger geographical areas, both at national and international level, leading to more widely applicable conclusions.

Supervision: PP YT ND AT.
Validation: PP YT ND AT.