Identifying Plasmodium falciparum transmission patterns through parasite prevalence and entomological inoculation rate
Abstract
Background:
Monitoring malaria transmission is a critical component of efforts to achieve targets for elimination and eradication. Two commonly monitored metrics of transmission intensity are parasite prevalence (PR) and the entomological inoculation rate (EIR). Comparing the spatial and temporal variations in the PR and EIR of a given geographical region and modelling the relationship between the two metrics may provide a fuller picture of the malaria epidemiology of the region to inform control activities.
Methods:
Using geostatistical methods, we compare the spatial and temporal patterns of Plasmodium falciparum EIR and PR using data collected over 38 months in a rural area of Malawi. We then quantify the relationship between EIR and PR by using empirical and mechanistic statistical models.
Results:
Hotspots identified through the EIR and PR partly overlapped during high transmission seasons but not during low transmission seasons. The estimated relationship showed a 1month delayed effect of EIR on PR such that at lower levels of EIR, increases in EIR are associated with rapid rise in PR, whereas at higher levels of EIR, changes in EIR do not translate into notable changes in PR.
Conclusions:
Our study emphasises the need for integrated malaria control strategies that combine vector and human host managements monitored by both entomological and parasitaemia indices.
Funding:
This work was supported by Stichting Dioraphte grant number 13050800.
Introduction
National malaria control programmes, working in collaboration with global stakeholders, have achieved extensive intervention coverage over the last two decades, leading to significant reductions in morbidity and mortality due to malaria (Bhatt et al., 2015a). However, malaria is still a leading global health problem. The previous successes and current challenges have motivated ambitious, yet feasible, global and national targets towards malaria elimination. A key component of efforts to achieve these targets is surveillance and monitoring, which is critical for continued assessment of intervention effectiveness, identification of areas or groups at the highest risk, and guiding the development and implementation of new intervention strategies (World Health Organization, 2015).
A wide range of metrics exists for monitoring malaria parasite transmission. The strengths and limitations of each metric are related, in part, to the step of the parasite transmission cycle it measures (Tusting et al., 2014). These strengths and weaknesses, including the sensitivity of each metric, vary across epidemiological settings and as parasite transmission declines within a given setting (The malERA Refresh Consultative Panel on Characterising the Reservoir and Measuring Transmission, 2017). Two of the most commonly monitored metrics are the prevalence of Plasmodium parasites and the entomological inoculation rate (EIR), especially in moderate to high transmission settings.
The prevalence of Plasmodium parasites in the human population at a given time point (i.e. the parasite rate; PR) approximates the reservoir of hosts potentially available to transmit the parasite from humans to mosquitoes. While only the gametocyte stage of the parasite contributes to transmission, it remains relatively expensive to detect this stage of the parasite. Conversely, rapid diagnostic tests (RDTs) primarily detect asexualstage antigens, yet they are inexpensive and easily deployed in largescale communitybased surveys (Poti et al., 2020). Still, the limit of detection (50–200 parasites/$\mathrm{\mu}\mathrm{l}$) for RDTs is higher than that of expert microscopy or PCR (Chiodini, 2014), so that RDTbased estimates of PR are biased by excluding lowdensity infections. Despite these limitations, RDTbased crosssectional surveys to measure PR capture both symptomatic and asymptomatic infections, which is important because both are likely to contribute to transmission (Bousema et al., 2014; Slater et al., 2019), and changes in PR over time can indicate changes in transmission.
EIR provides an estimate of the intensity of parasite transmission from mosquitoes to humans, expressed as the number of infectious bites received per person per unit time. EIR is calculated by multiplying the number of malaria vector bites per person per unit time, also known as the human biting rate (HBR), by the proportion of vectors carrying the infectious sporozoite stage of malaria parasites, referred to as the sporozoite rate (SR) (Onori and Grab, 1980). The accuracy and precision of EIR estimates, therefore, depends on the accuracy and precision with which HBR and SR can be measured (Tusting et al., 2014). Two common methods for measuring HBR are the human landing catch and the Centers for Disease Control and Prevention Light Trap, but interindividual variation in attractiveness to mosquitoes restricts standardisation across sampling points for both of these methods (Knols et al., 1995; Qiu et al., 2006). Alternative methods for estimating HBR include the Suna trap, which uses a synthetic blend of volatiles found on human skin and carbon dioxide to attract hostseeking Anopheles mosquitoes (Mukabana et al., 2012; Menger et al., 2014; Hiscox et al., 2014). The standardised odour blend allows for reliable comparisons among trapping locations (Mburu et al., 2019). Regardless of the method used to estimate HBR, the precision of SR decreases as the number of mosquitoes collected decreases. Despite these limitations, EIR is a vital metric of malaria parasite transmission because it directly describes human exposure to malaria parasites before postinoculation factors such as immunity, nutrition, and access to health care (Killeen et al., 2000). Moreover, EIR provides information about the relative contributions of different vector species to transmission, which can impact malaria intervention effectiveness based on interspecies differences in biting behaviours related to time and location, nonhuman bloodmeal hosts, larval ecology, and insecticide resistance profiles (Ferguson et al., 2010).
Malaria parasite transmission is heterogeneous in space and time at fine resolution due to several factors, including the availability of larval mosquito habitat, socioeconomics, human behaviour and genetics, and malaria intervention coverage (Carter et al., 2000; Bousema et al., 2012; McCann et al., 2017a). Repeated crosssectional surveys continuously carried out in communities can reveal this fineresolution heterogeneity (RocaFeltrer et al., 2012), providing timely estimates of malaria control progress at the subdistrict level and potentially identifying hotspots of malaria parasite transmission for targeted intervention (Kabaghe et al., 2017; Bousema et al., 2016). However, understanding this heterogeneity and identifying hotspots in a way that is meaningful for control programmes remains challenging (Stresman et al., 2019), in part because hotspot location and size can depend on which metric is used (Stresman et al., 2017). Given that PR and EIR are indicative of components of the parasite transmission cycle that are separated by multiple complex steps, each metric provides partial but useful information about the underlying risk of transmission. Therefore, measuring and mapping both metrics can provide a fuller picture of parasite transmission (Cohen et al., 2017).
Additionally, modelling the functional relationship between EIR and PR can provide further insights into the underlying malaria epidemiology. Previous studies have demonstrated that this relationship is nonlinear, such that small changes in EIR are associated with large changes in PR when EIR is low, but PR saturates rather than changing at a constant rate when EIR is high (Beier et al., 1999; Smith et al., 2005). These previous studies were metaanalyses using paired estimates of EIR and PR, with one estimate of each outcome per site, from sites representing a wide range of EIR and PR in Africa. Their findings had a number of important implications, which included providing estimated ranges for the change in PR that may be expected for a given change in EIR. However, these estimates implicitly assumed that the relationship is constant across space on a continental scale, such that differences in EIR and PR between sites would be indicative of differences over time within a site. Yet no previous study has explicitly examined this relationship over time within a single geographical region.
In the current study, we use a series of repeated crosssectional surveys conducted over 38 months in one region of southern Malawi to map the finescale spatiotemporal dynamics of P. falciparum entomological inoculation rate (PfEIR) and P. falciparum parasite prevalence (PfPR). The joint monitoring of these two outcomes in space and time allows us to identify and compare the spatial heterogeneities and temporal patterns of PfEIR and PfPR in a region with moderately intense, seasonally variable malaria parasite transmission. We then investigate the PfEIRPfPR relationship based on changes in these outcomes observed at both annual and subannual scales within our study site using several statistical models, which can be distinguished as follows: mechanistic models that are based on different epidemiological assumptions and empirical models where the data inform the PfEIRPfPR relationship. These approaches allows us to address the following questions. (1) How do spatiotemporal patterns of EIR and PR compare? (2) Do EIR and PR lead to the identification of the same malaria hotspots? (3) As EIR changes over time, how do those changes in EIR affect PR? (4) Does EIR have a lagged effect on PR? (5) Does the EIRPR relationship vary between women of reproductive age and children between 6 and 59 months of age?
Materials and methods
Study site
Request a detailed protocolThis study was part of the Majete Malaria Project (MMP), an integrated malaria control project in Chikwawa District, Malawi. The catchment area of MMP consisted of three distinct geographical regions, referred to as Focal Areas A, B, and C (Figure 1), with a total population of about 25,000 people living in 6600 households in 65 villages.
Chikwawa experiences highly variable rainfall during its single rainy season, which spans November/December to April/May. Temperatures are generally high, with daily maximum temperatures in December averaging 37.6°C, and in July averaging 27.6°C (Joshua et al., 2016). A wide range of both permanent and temporary water bodies create suitable larval habitats in the region for Anopheles funestus s.s., Anopheles arabiensis, and Anopheles gambiae s.s., including dams, swamps, ponds, borehole runoffs and drainage channels (Gowelo et al., 2020).
Malaria control in the district is implemented through the Chikwawa District Health Office. During the study period, interventions applied throughout the study area included the continuous provision of insecticidetreated nets (ITNs) to pregnant women and children under five years old, mass distribution campaigns of ITNs targeting universal coverage, intermittent preventative therapy for pregnant women, and malaria case diagnosis and treatment with artemisininbased combination therapy. The only mass distribution of ITNs in the district during the study period occurred in April 2016. As part of the MMP, a randomised trial was conducted to assess the effectiveness of additional, communityimplemented malaria interventions between May 2016 and May 2018 (McCann et al., 2017b). The trial interventions were implemented at the village level, with villages assigned to one of four groups: (a) no additional interventions; (b) larval source management; (c) house improvement; and (d) both larval source management and house improvement (McCann et al., 2017b; van den Berg et al., 2018).
Data
To quantify PfPR and PfEIR over the course of the study, a rolling malaria indicator survey (rMIS) (RocaFeltrer et al., 2012) was conducted in conjunction with mosquito sampling, forming a series of repeated cross sectional surveys. Sampling was carried out over 17 rounds, with each round spanning a period of 2 or 3 months. In the first two rounds of data collection (baseline, from April through August 2015), an inhibitory geostatistical sampling design (IGSD) was used to select 300 and 270 households, respectively, for the rMIS from an enumeration database of all households in the catchment area (Chipeta et al., 2017). The IGSD helped to ensure that randomly sampled households are relatively uniformly spaced over the study region by requiring each pair of sampled households to be separated by a distance of at least 0.1 km, which increases the efficiency of hotspot detection (Kabaghe et al., 2017). In the three subsequent rounds of data collection during the baseline, an adaptive geostatistical sampling design (AGSD) was used to select 270 households per round (Chipeta et al., 2016). With AGSD, new households for the current round of rMIS were chosen from regions with high standard errors of estimated prevalence, based on data from all previous rounds. In the baseline period, previously sampled households were not eligible for sampling in subsequent rounds. For the trial period (starting May 2016), IGSD was again used to select households from the enumeration database of all households. All households were eligible for selection in each round of the trial period regardless of whether they were selected in a previous round. At each round of rMIS data collection in the baseline and trial phases, respectively, 75% and 72% of the households chosen at each round of the rMIS were then randomly selected for mosquito sampling.
In each sampled household, children under five (0.5–5 y/o) and women of reproductive age (15–49 y/o) were tested for P. falciparum using an RDT (SD BIOLINE Malaria Ag P.f. HRPII, Standard Diagnostics, Yonginsi, Republic of Korea).
Mosquitoes were sampled from 5pm to 7am using Suna traps (Biogents AG, Regensburg, Germany) with MB5 blend plus $C{O}_{2}$ to mimic human odour (Hiscox et al., 2014; Mburu et al., 2019). For a selected household in a surveillance round, the trap was set for one night indoors and one night outdoors, with the order of indoor/outdoor determined randomly. For households where the residents were sleeping in more than one building, a trap was set at each building. Trapped female anophelines were preserved using a desiccant and identified using standard morphological and molecular techniques (Gillies and Coetzee, 1987; Koekemoer et al., 2002; Scott et al., 1993). Female anophelines were further tested for the presence of P. falciparum in their head and thorax, after removing the abdomen, using quantitative polymerase chain reaction (qPCR ) (Bass et al., 2008; Perandin et al., 2004). Specimens with a Ct value below 37.0 were considered positive for P. falciparum.
Environmental and climatic factors
Request a detailed protocolEnvironmental and climatic factors affect the abundance and suitability of water bodies that support the development of immature mosquitoes (Madder et al., 1983; Loetti et al., 2011), the duration of mosquito development (Ciota et al., 2014; Loetti et al., 2011; Craig et al., 1999), mosquito hostseeking and biting behaviour, and the development rate of malaria parasites in mosquitoes (Rumisha et al., 2014; Amek et al., 2011).
Using hourly measurements of temperature and relative humidity (RH) from a weather station in each focal area, we computed the average temperature and RH for different ranges of days before the day of data collection (Appendix 1 – Procedure for building the HBR, PfSR and PfPR models).
Spectral indices, namely normalised difference vegetation index (NDVI) and enhanced vegetation index (EVI), were computed using remotely sensed multispectral imagery from the Landsat 8 satellite. These data are freely available from the United States Geological Survey (USGS) Earth Explorer (earthexplorer.usgs.gov) as raster files at a spatial resolution of 30 × 30 m for every 16 days. For our analysis, we averaged each spectral index over 5 years, from April 2013 to April 2018, while omitting scenes that were dominated by cloud artefacts.
We extracted raster data of surface elevation from the global digital elevation model (DEM) generated using measurements from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) (Tachikawa et al., 2011). These data are also freely available for download from the USGS Earth Explorer. Using a flow accumulation map derived from the DEM, a river network map was generated and used to calculate and store as raster images the distance to small rivers and large rivers (henceforth, DSR and DLR, respectively).
Geostatistical analysis
The number of mosquitoes trapped by Suna traps can be used to estimate HBR, as these traps primarily target hostseeking mosquitoes. Hence, we first estimated HBR and the P. falciparum sporozoite rate (PfSR). We then estimated PfEIR as the product of these two quantities.
We carried out separate analyses for A. arabiensis and A. funestus s.s., using explanatory variables and random effects structures that we found to be suitable for each species. Details of the variable selection process and the final sets of explanatory variables for each of the models later described in this section are given in Appendix 1 – Procedure for building the HBR, PfSR, and PfPR models. The correlation structures adopted for the geostatistical models were informed by the variogrambased algorithm described in Giorgi et al., 2018. The geostatistical models for the HBR and PfPR data described below were fitted using PrevMap (Giorgi and Diggle, 2017), freely available from the Comprehensive R Archive Network (CRAN, www.rproject.org). The PfSR models were fitted using the glmm package, also available on CRAN.
Human biting rate
Request a detailed protocolLet $Y({x}_{i},{t}_{i}),i=1,\mathrm{\dots},M$, where $M=2432$ is the total number of houses, denote counts of mosquitoes trapped at location x_{i} in month ${t}_{i}\in \{1,\mathrm{\dots},38\}$, where ${t}_{i}=1$ denotes April 2015. We modelled the $Y({x}_{i},{t}_{i})$ using Poisson mixed models expressed by the following linear predictor
where: $d({x}_{i},{t}_{i})$ is a vector of spatiotemporal explanatory variables with associated regression coefficients β; the $f({t}_{i};\alpha )$ is a linear combination of several functions of time, including sines, cosines and splines, with an associated vector of regression parameters α, accounting for trends and seasonal patterns; the ${Z}_{i}$ are independent and identically distributed Gaussian random variables with variance ${\tau}^{2};\text{}S(x)$ is a zeromean stationary and isotropic Gaussian process with variance ${\sigma}^{2}$ and exponential correlation function $\rho (u)=\mathrm{exp}(u/\varphi )$, where $\varphi $ regulates the pace at which the spatial correlation decays for increasing distance $u$ between any two locations. We allow the explanatory variables $d({x}_{i},{t}_{i})$ and $f({t}_{i};\alpha )$ to differ between different mosquito species since different species may respond differently to environmental changes. We point out that the stationarity of the process $S(x)$ implies that all of its properties, including the variance (${\sigma}^{2}$) and scale of the spatial correlation ($\varphi $), are constant over space. The estimation of the model parameters is then carried out using Monte Carlo Maximum Likelihood (Christensen, 2004).
Plasmodium falciparum sporozoite rate
Request a detailed protocolLet ${Y}^{*}({x}_{i},{t}_{i})$ be the number of mosquitoes that tested positive for the presence of P. falciparum sporozoites. We assumed that the ${Y}^{*}({x}_{i},{t}_{i})$ follow a Binomial mixed model with number of trials ${N}^{*}({x}_{i},{t}_{i})$, that is the total number of successfully tested mosquitoes, and probability of testing positive $PfSR({x}_{i},{t}_{i})$. We model the latter as a logitlinear regression given by
where each term in (2) has an analogous interpretation to those of (1). A spatial process $S(x)$ was not included in the sporozoite rate model because we found no evidence of residual spatial correlation in the sporozoite rate data (Appendix 1—figure 1).
Estimating the Plasmodium falciparum entomological inoculation rate
Request a detailed protocolLet $PfEI{R}_{f}(x,t)$ and $PfEI{R}_{a}(x,t)$ denote the PfEIR for A. funestus s.s. and A. arabiensis at a given location $x$ and month $t$. We estimated each of these two as
where $l(t)$ is the number of days in month $t$. Finally, we estimated the overall PfEIR as
We then mapped PfEIR as in (3) over a 30 × 30 m regular grid covering the whole of the study area for each month across 38 months.
To map PfEIR for each month, we first simulate at each prediction location (i.e. the centroid of each grid cell) 10,000 samples from the conditional distribution of the random effects (corresponding to $S(x)+Z$ in the case of the PfHBR and $Z$ in the case of PfSR) given the data. We then transform these to obtain 10,000 predicted surfaces for PfHBR and PfSR, and by applying (3) we obtain 10,000 predictive samples for PfEIR. The predicted PfEIR at each prediction location is taken to be the median of the 10,000 samples at that location. The associated 95% predictive interval is the 2.5th to 97.5th percentile of the 10,000 predictive samples.
In this procedure, all the parameters corresponding to the regression coefficients, the scale and variance of the spatial process, and variance of Gaussian noise are fixed at their MCML estimates.
Plasmodium falciparum prevalence
Request a detailed protocolWe mapped PfPR in women and in children by fitting a geostatistical model to each group. More specifically, let $I({x}_{i},{t}_{i})$ denote the number of RDT positives out of ${N}_{it}$ sampled individuals at location x_{i} in month t_{i}. We then assumed that the $I({x}_{i},{t}_{i})$ follow a Binomial mixed model with probability of a positive RDT result $p({x}_{i},{t}_{i})$, such that
where $T({x}_{i})$ is a stationary and isotropic Gaussian process with exponential correlation function and ${U}_{i}$ are Gaussian noise, $g({t}_{i},\varrho )$ is a linear combination of splines, and sine and cosine functions of time accounting for trends and seasonality, and φ and $\varrho $ are vectors of regression parameters to be estimated.
Hotspot detection using PfEIR and PfPR
Request a detailed protocolWe demarcated hotspots for PfEIR and PfPR using an exceedance probability approach. Using the resulting 10,000 predictive samples for PfEIR and PfPR, as described above, we then obtained the exceedance probability for each outcome at each spacetime location by computing the proportion of the 10,000 predictive samples that exceeded the respective, predefined thresholds, which were set at 0.1 ib/person/month for EIR, 31% for PfPR in children, and 17% for PfPR in women. Finally, we mapped these exceedance probabilities and demarcated hotspots as areas where these probabilities were at least 0.9. The PfPR thresholds were defined to correspond to the PfEIR threshold based on the best of six functional relationships between PfEIR and PfPR as described in the next section.
Modelling the relationship between PfEIR and PfPR
Because PfEIR may have a delayed effect on PfPR, possibly due to the time taken for P. falciparum to develop in the human host, we considered that current PfPR may depend on PfEIR $l$ months prior. In particular, we considered $l=0,1,2$. We then assumed that the number of RDT positive individuals, $I({x}_{i},{t}_{i})$, follow independent Binomial distributions such that
where $h(\cdot )$ is a function depending on a vector of parameters θ that governs the relationship between PfPR and PfEIR, and $\widehat{PfEIR}({x}_{i},{t}_{i}l)$ is the estimated PfEIR as in Equation (3). We considered six models, each of which provided a different specification for $h(\cdot )$.
We now describe the six models for $h(\cdot )$. Models 1 to 4 make explicit assumptions on the underlying mechanism of transmission, whereas models 5 and 6 describe the functional relationship between PfEIR and PfPR through regression methods.
Model 1: The susceptibleinfectedsusceptible (SIS) model
Request a detailed protocolLet $b$ be the probability that an infectious mosquito bite results in an infection, referred to as the transmission efficiency. Then, infections at $({x}_{i},{t}_{i}l)$ are assessed to occur at a rate of $b\times PfEIR({x}_{i},{t}_{i}l)$. We assumed that each infection cleared independently over a duration $1/r$ so that the ratio $\gamma =b/r$ is the time taken to clear infection per infectious bite. We assumed that the relationship between PfEIR and PfPR holds throughout the study region. If $PfEIR(x,tl)$ is constant, the relationship between $PfEIR(x,tl)$ and $PfPR(x,t)$ is described by Ross, 1911
We obtained our first model as the nonzero equilibrium solution of (6), given by
Model 2: The SIS model with different infection/recovery rates (D.I/R)
Request a detailed protocolModel 1 assumes that women and children get infected and recover at the same rate. However, the transmission and recovery rates in children may differ from those in women. We, therefore, modified Model 1 by allowing different values of $b$ and $r$ for each category of people. Let ${\xi}_{1,it}$ and ${\xi}_{2,it}$ respectively be the proportion of children and women sampled at $({x}_{i},{t}_{i})$ and ${\gamma}_{k}={b}_{k}/{r}_{k}$, where $k=1$ denotes children and $k=2$ denotes women. The resulting Model 2 is
Model 3: The SIS model with superinfection (S.I.)
Request a detailed protocolIf individuals are superinfected with P. falciparum, then the rate at which infections clear depends on the infection rate, with clearance being faster when infection rate is low, and slower when infection rate is high. To capture this feature, we modelled infection clearance rate as $g(\vartheta ,r)=\vartheta /({e}^{\vartheta /r}1)$, where $\vartheta =b\times PfEIR$(Smith et al., 2005; Walton, 1947; Dietz et al., 1974; Aron and May, 1982). The resulting model for $PfPR(x,t)$ is
Model 4: The SIS model with S.I and D.I/R
Request a detailed protocolCombining the assumptions of heterogeneous infection/recovery rates, as in Model 2 and superinfection, as in Model 3, we obtain Model 4,
Model 5: The Beier model
Request a detailed protocolBeier et al., 1999 assumed that the log of PfEIR is linearly related to PfPR, and fitted the regression model
the so called ‘loglinear model’.
Model 6: The logitlinear model
Request a detailed protocolThe Beier model has the limitation that PfPR approaches $\mathrm{\infty}$ as PfEIR goes to 0 and approaches $\mathrm{\infty}$ as PfEIR goes to $\mathrm{\infty}$. To constrain PfPR to lie between 0 and 1, we applied the logitlink function to PfPR to give Model 6,
Parameter estimation of the PfEIRPfPR relationship models
Request a detailed protocolWe estimated the parameters of each of the six models by maximising the loglikelihood function
To fit each model, we first obtained 10,000 bootstrapped data sets of predicted PfEIR as in (3) at the set of all spacetime locations sampled for the rMIS. We did this for two reasons: to obtain PfEIR data at locations $({x}_{i},{t}_{i})$ that were sampled for rMIS but not for the entomological surveillance; and to account for the uncertainty in PfEIR. The predicted PfEIR values were then paired with respective empirical PfPR values at $({x}_{i},{t}_{i})$. By fitting each model to each of the 10,000 datasets, we then obtained 10,000 bootstrapped samples $\{{\widehat{\theta}}_{1},\mathrm{\dots},{\widehat{\theta}}_{10000}\}$ for the vector of parameter estimates $\widehat{\theta}$ of each the six candidate models. We then summarised these samples by their mean and central 95% probability interval. We repeated this process for $l=0,1,2$.
We compared the fit of the six models based on their predictive ability as measured by the bias and rootmeansquare error when each model is used to predict prevalence at all the observed spacetime locations.
Results
rMIS and mosquito sampling
From April 2015 to May 2018, a total of 6870 traps (3439 indoors; 3431 outdoors) were placed at 2432 houses over 17 rounds of sampling (Figure 2), resulting in the collection of 657 female Anopheles mosquitoes (Table 1). Following PCR of the 478 A. gambiae s.l. collected, 92% were identified as A. arabiensis, 2% as A. gambiae s.s., 1% as A. quadriannulatus, and 5% could not be identified further. From the 179 A. funestus s.l. collected, 95% were identified as A. funestus s.s. by PCR, while the remaining 5% could not be identified further. The observed vector composition is therefore 71%, 27%, and 2% for A. arabiensis, A. funestus s.s., and A. gambiae s.s., respectively.
Despite the relatively low abundance of A. funestus s.s. compared to A. arabiensis, the higher sporozoite rate of the former made the contribution of A. funestus s.s. to the total PfEIR almost equivalent to that of A. arabiensis (Table 1). The total PfEIR for the 38 months was 8.24 ib/person, equivalent to an average 2.60 ib/person/year.
Over the same 38month period, 5685 individual P. falciparum RDT tests were conducted across 3096 household visits (Figure 2). Among the 2401 tests conducted on children aged 6 to 59 months, 25.5% were positive, while 14.3% of the 3284 tests conducted on women aged 15–49 y/o were positive.
Spatiotemporal patterns of PfEIR and PfPR
We observed clear spatiotemporal heterogeneities in PfEIR, PfPR in children, and PfPR in women when mapped across the study region at a fine spatial resolution (30 x 30 m) and 1 month intervals. For convenient visualisation of the main features of the spatiotemporal maps, we have developed an interactive webbased application to show the maps at http://chicas.lancasteruniversity.uk/projects/malaria_in_malawi/pfpr/. We show selected predictive maps of PfEIR and PfPR in Figure 3 and exceedance probability of PfEIR and PfPR in Figure 4 for June 2015, August 2016 and November 2017, which are representative of high, medium and low transmission months, relative to the full study period. Spatially, there were differences both within and between the three focal areas. Focal Area A generally showed the lowest PfEIR and PfPR, while Focal Areas B and C showed similar, higher levels of PfEIR. Within each focal area, the spatial patterns changed from month to month, with hotspots of both PfEIR and PfPR proceeding through seasonal cycles of expansion and retraction over time. Over the 3year study period, hotspots of PfEIR and PfPR generally disappeared during the low transmission seasons, except for residual hotspots of PfPR that persisted throughout the study period.
When summarised over the whole study region, each of PfPR and PfEIR exhibited seasonal patterns with a single annual peak. The monthly predicted PfEIRs and PfPRs were similar to the observed values (Figure 5). PfEIR increased from November to a peak in May and decreased to a trough in November. PfPR started increasing from December to a peak around July, after which it decreased to a trough between November/December.
Three observations are clear from both the spatiotemporal maps and the monthly summarised data (Figure 5). First, children aged 6–59 months consistently had a higher level of PfPR than women throughout the study period. Second, PfPR in both groups generally decreased from the start of the study in April 2015 to December 2016, after which there was a general trend of increasing PfPR in both age groups. Finally, PfEIR was steady in the first 2 years of the study, followed by a general decrease after May 2016. Strikingly, the observed PfEIR was 0 between June 2017 and the end of the study, while the PfPR increased in both children and women between November 2017 and May 2018.
The relationship between PfEIR and PfPR
Temporally, the seasonal patterns of PfEIR and PfPR within each year were nearly concurrent, with the estimated peak in PfEIR preceding that of PfPR by one month (Figure 5).
Spatially, PfEIR and PfPR showed broadly similar patterns. When comparing the hotspots of PfEIR and PfPR using spatiotemporal maps of exceedance probabilities, the hotspots of PfEIR and PfPR partly overlapped during the high transmission season (http://chicas.lancasteruniversity.uk/projects/malaria_in_malawi/pfpr/). However, there were hotspots of PfEIR that were not necessarily hotspots of PfPR and vice versa (Figure 6).
Scatter plots of the logit of PfPR against the log of PfEIR show an approximately direct linear relationship (Figure 7).
For each of the six classes of model, those with a 1month laggedeffect were found to be the best based on the rootmeansquareerror (RMSE) and bias indices of predictive performance. Among these, the empirical models (i.e. logitlinear and Beier) yielded lower values for the bias and RMSE values than the mechanistic models (Appendix 1 – Procedure for building the HBR, PfSR, and PfPR models, Appendix 1—table 6). The logitlinear model, albeit by a small margin, outperformed all the models (Appendix 1 – Procedure for building the HBR, PfSR, and PfPR models, Appendix 1—table 6).
The fitted logitlinear model (Figure 8) shows that PfPR rises quickly with increasing PfEIR at very low PfEIR, followed by a flattening off or saturation. From the estimated relationship for women and children combined (Figure 8a), a decrease in PfEIR from 1 ib/person/month to 0.85 ib/person/month (i.e. a 0.15 decrease in PfEIR) is associated with a reduction in PfPR from 27.17% to 26.85% on average (i.e. a 1.18% decrease in PfPR). However, at the lower ranges of EIR, the same decrease of 0.15 ib/person/month from 0.2 ib/person/month to 0.05 ib/person/month is associated with a reduction in PfPR from 24.10% to 21.66% on average (i.e. a 10.13% decrease in PfPR). The resulting nonlinear reltionship emphasizes two aspect of the PfPRPfEIR relationship. First, for large values in PfEIR, reductions in PfEIR are associated with smaller reductions in PfPR, whilst if PfEIR is low, reductions in PfEIR are associated with greater reductions in PfPR. Second, even when transmission, as measured by PfEIR, has been driven close to zero, PfPR can still remain substantial.
An indication of possible differences in the PfEIR–PfPR relationship between children and women lies in the logitlinear model fitted to children and women separately (Figure 8b). The average trajectories of PfPR and corresponding 95% confidence intervals with varying PfEIR are distinct for women and children. PfPR in children tends to show a steeper rise with increasing PfEIR than in women. From the estimated relationship for children, a decrease in PfEIR from 0.1 ib/person/month to 0.01 ib/person/month is associated with a reduction in PfPR from 31.07% to 25.52% on average (i.e. a 17.86% decrease in PfPR). From the estimated relationship for women, the same decrease in PfEIR is associated with a reduction in PfPR from 16.84% to 14.33% (i.e. a 14.90% decrease in PfPR) on average. We make two observations. (1) With decreasing PfEIR, the percentage reduction in PfPR achieved in children tends to be higher than in women. (2) When transmission has been driven almost to zero, PfPR remains consistently high in children.
Discussion
Using data from 38 months of repeated crosssectional surveys, we have mapped the finescale spatiotemporal dynamics of PfEIR and PfPR in a region of Malawi with moderately intense, seasonally variable malaria parasite transmission. We found evident spatial heterogeneity in both PfEIR and PfPR, with areas of higher PfEIR and PfPR expanding and contracting over time. We also found that hotspots of PfEIR and hotspots of PfPR overlapped at times, but the amount of overlap varied over time. Finally, we showed that monthtomonth variations in PfEIR over the study period are strongly associated with changes in PfPR. These findings highlight the dynamic nature of malaria parasite transmission and underscore the value of monitoring both PfEIR and PfPR at fine spatial and temporal resolutions.
In this study, we used modelbased geostatistics (Diggle and Giorgi, 2019) to analyse repeated crosssectional data with a unique sampling framework (RocaFeltrer et al., 2012; Kabaghe et al., 2017) and estimate the finescale spatial patterns of PfEIR and PfPR across our study site at monthly intervals. This approach was essential for identifying hotspots of PfEIR and PfPR through the use of exceedance probabilities, because it allowed us to acknowledge the effects of unmeasured malaria risk factors on each metric through the inclusion of spatially structured random variations. Consequently, we were able to directly compare the spatial patterns of PfEIR and PfPR over an extended period of time in a single geographical region. As expected, there were hotspots identified by each of the two metrics of malaria transmission risk, which expanded and retracted over time. However, the hotspots of PfEIR and PfPR only partially overlapped, with the most substantial amount of overlap observed during the high transmission seasons. Withinvillage and betweenvillage spatial heterogeneities of malaria parasite transmission are well documented across many sites (Greenwood, 1989; Thompson et al., 1997; Amek et al., 2012; Mwandagalirwa et al., 2017), but few previous studies have compared the spatial distributions of different transmission metrics in the same site (but see Stresman et al., 2017). Given that all available metrics of malaria parasite transmission have their own strengths and weaknesses (Tusting et al., 2014), our findings suggest that monitoring multiple transmission metrics, each aligned with widely separated steps of the parasite transmission cycle, provides a more complete understanding of the underlying spatial heterogeneity in malaria parasite transmission (Cohen et al., 2017). Furthermore, monitoring multiple metrics may provide an opportunity to optimise the impact of control activities by targeting different control activities to different locations based on differences in the metrics (Cohen et al., 2017). For example, areas with higher transmission risk according to an entomological metric (e.g. PfEIR) than a measure of the potential transmission reservoir (e.g. PfPR) may indicate a need for increased vector control, whereas areas with lower PfEIR and higher PfPR may indicate a need for increased treatment of malaria parasite infections.
Our geostatistical modelling approach also provided a principled framework for estimating PfEIR in our study, allowing us to robustly investigate the PfEIRPfPR relationship despite the low mosquito densities observed in this region. Previous studies Beier et al., 1999; Smith et al., 2005 have estimated the relationship between PfEIR and PfPR using paired estimates of these metrics from several sites throughout Africa, characterised by a wide range of transmission intensities (PfEIR <1 to >500 ib/person/year). These studies demonstrated that small changes in PfEIR are associated with large changes in PfPR when PfEIR is low, for example below about 15 ib/person/year (Beier et al., 1999). However, estimating PfEIR in settings with low mosquito density is challenging because the accuracy and precision of PfEIR depend on the accuracy and precision of the human biting rates and sprorozite rates used to calculate PfEIR (Tusting et al., 2014). We overcame these challenges by using modelbased geostatistics to improve the precision of our PfEIR estimates and, just as importantly, used bootstrap procedures to propagate the uncertainty from each modelling step to the next. Altogether, these methods allowed us to investigate the PfEIRPfPR relationship with a focus on much lower ranges of PfEIR than previous studies.
An additional advantage of using geostatistical models in this study was the prediction of entomological data at unsampled geographical locations. For a number of logistical reasons (e.g. mosquito sampling was conducted over 2 consecutive nights at each sampled location), we sampled for mosquitoes at roughly 75% of the locations visited for parasitaemia sampling in each round. For the geographical locations where empirical parasitaemia data were available but entomological data were not, our geostatistical model predictions of PfEIR were combined with the empirical PfPR data. This rigorous statistical solution enriched the data used in our assessment of the PfEIRPfPR relationship.
Prior to our study, the most recent assessment of PfEIR in this district of Malawi was from 2002, with an estimated annual PfEIR of 183 ib/person/year (Mzilahowa et al., 2012). The drastic reduction in annual PfEIR since then to an estimated 2.60 ib/person/year in our study is likely due, at least in part, to an increase in the use of ITNs and ACTs as observed elsewhere (Bayoh et al., 2010; Mwangangi et al., 2013). Nationwide, use of ITNs by children under 5 years old in Malawi has increased from nil in 2000 and 14.8% in 2004 (Mathanga et al., 2012) to 67.8% in 2014 (Malawi National Malaria Control Programme and ICF International, 2014). Treatment for malaria in Malawi switched from sulfadoxine–pyrimethamine to ACT with artemether–lumefantrine in 2007 (Mathanga et al., 2012), and by 2014, 39.3% of children under five reporting a fever had taken ACT (Malawi National Malaria Control Programme and ICF International, 2014). In addition to these longterm, nationwide trends, changes in malaria intervention coverage over time also likely impacted malaria parasite transmission in the more specific context of our study. Mass distributions of ITNs took place across Malawi in 2012 (World Health Organization, 2013), that is about 3 years prior to our study, and again in April 2016, that is, 1 year into our 38month study. Additionally, randomly selected villages in our study site implemented communityled larval source management, house improvement, or both as part of a randomised trial between May 2016 and May 2018 (McCann et al., 2017b; van den Berg et al., 2018). Although no differences in PfEIR or PfPR were observed between the trial arms (McCann et al., 2021), the entire study site, including the trial’s control arm, was included in the 2016 mass ITN distribution, as well as other National Malaria Control Programme interventions and a community engagement programme as part of the Majete Malaria Project (McCann et al., 2017b; van den Berg et al., 2018). Therefore, the changes in PfEIR and PfPR over time observed in this study reflect the combined effects of seasonal weather patterns, yeartoyear climatic variation and all malaria control activities.
Within this context of observing monthtomonth changes in both PfEIR and PfPR in a single geographical region, we have demonstrated that fluctuations in PfEIR over short periods are associated with predictable changes in PfPR in the same region. We found that a logitlinear model explained the PfEIRPfPR relationship better than any of the other five model classes examined, and our data better supported a onemonth delayed effect of PfEIR on PfPR than no delayed effect or a 2month delayed effect. The onemonth delayed effect is likely due to the incubation period of the parasite (Ruan et al., 2008) and the duration of infections (Felger et al., 2012). As shown in previous studies (Beier et al., 1999; Smith et al., 2005), we observed that small changes in PfEIR led to relatively large changes in PfPR at lower ranges of PfEIR, while PfPR saturated rather than changing at a constant rate at higher ranges of PfEIR. These previous studies were based on estimates of PfEIR and PfPR from 31 locations (Beier et al., 1999) and more than 90 locations (Smith et al., 2005) representing a wide range of PfEIR and PfPR in Africa. By assuming that the PfEIRPfPR relationship is constant across space on a continental scale, results from these previous studies suggested that variation between PfEIR and PfPR across geographical location is representative of variation between PfEIR and PfPR over time. In our study, we explicitly confirmed this association between PfEIR and PfPR within a single location. Whilst we also assumed that the PfEIRPfPR relationship does not change dynamically across space and time, we believe this to be a realistic assumption for the restricted geographical setting of our study.
The saturation in PfPR with increasing PfEIR may be explained by several factors, which are not mutually exclusive. One set of factors relates to people being heterogeneously exposed to vectors (Guelbéogo et al., 2018) because of differences in attractiveness (Knols et al., 1995; Qiu et al., 2006), behaviour (SherrardSmith et al., 2019; Finda et al., 2019), access to ITNs (Bhatt et al., 2015b), housing design (Tusting et al., 2015; Tusting et al., 2017), or the spatial distribution of vector habitat (McCann et al., 2017a), so that as PfEIR increases, it is more likely that infectious vectors are biting already infected individuals (Smith et al., 2007b; Smith et al., 2010). The second set of factors relates to interindividual variation in acquired immunity, which in some individuals may prevent vectorinoculated sporozoites from progressing to bloodstage infection (John et al., 2005; Offeddu et al., 2017), keep bloodstage infections at densities lower than the level of detection (Doolan et al., 2009) (about 50–200 parasites/µl for RDTs as used in our study), or increase the rate at which bloodstage infections are cleared (Hviid et al., 2015).
Regardless of the underlying factors driving the PfEIRPfPR relationship, our results have practical implications for the selection and interpretation of malaria parasite transmission metrics. In settings with higher ranges of PfEIR, moderate changes in PfEIR will not be associated with appreciable changes in PfPR (Beier et al., 1999; Smith et al., 2005; Churcher et al., 2015). Framed in terms of a public health goal to decrease PfPR in these settings with high baseline PfEIR, relatively large reductions in PfEIR would be required to achieve appreciable reductions in PfPR. In terms of selecting an appropriate metric for monitoring changes in transmission in these same highPfEIR settings, PfPR may only be suitable for measuring very large changes in transmission.
Conversely, in settings with a lower range of PfEIR, our results show that PfPR is sensitive to smaller, shortterm changes in malaria parasite transmission. This finding highlights the importance of sustaining vector control efforts in settings with relatively low PfEIR, because a small increase in the rate of infectious bites (PfEIR) could result in a rapid increase in the proportion of people infected (PfPR). This sensitivity of PfPR to shortterm changes in parasite transmission, when PfEIR is low, also confirms that PfPR can be used for monitoring changes in the intensity of parasite transmission linked to either environmental conditions or the effects of malaria interventions. However, this sensitivity of PfPR in these settings suggests that annual crosssectional surveys aiming for a transmission peak are more likely to miss the actual peak than in settings with higher parasite transmission intensity, as shown previously by others (Kang et al., 2018). National malaria control programs and others planning malaria indicator surveys to measure yeartoyear variation in PfPR should therefore consider approaches to identify and account for any potential bias in PfPR estimates from a single time point, for example incorporating continuous or ‘rolling’ surveys at sentinel sites (RocaFeltrer et al., 2012; Kabaghe et al., 2017), monitoring ‘easyaccess groups’ (Sesay et al., 2017), or modelling subannual trends based on health facility data (Sturrock et al., 2014; Awine et al., 2018). These considerations likely apply to settings where increasing coverage of ITNs (Bhatt et al., 2015a) and ACTs (Bennett et al., 2017) has reduced PfPR from ≥50% (i.e. holo and hyperendemic [Hay et al., 2008]) to between 10–50% (i.e. mesoendemic), which have become increasingly common over the last 20 years (Weiss et al., 2019).
The monthly PfEIR in our study region was 0 ib/person/month in multiple months. This likely indicates that the number of infectious bites received per person during these months was below the level of detection, rather than an actual interruption of transmission during those months, especially in the first 2 years of the study when these periods only lasted 2–3 months. Our finding that a monthly PfEIR near or equal to zero is associated with substantial PfPR is, therefore, unsurprising given that previous studies have had similar findings when comparing annual PfEIR to PfPR (Kabiru, 1994; Mbogo et al., 1995; Beier et al., 1999; Smith et al., 2005). On the other hand, we observed an increase in PfPR from about November 2017 to May 2018 while PfEIR remained at zero. It remains unclear whether this rise in PfPR was due to new infectious bites that nevertheless remained below the level of detection or to previously infected individuals with parasite densities that increased to detectable levels (Drakeley et al., 2018). Either way, this result shows that a rise in PfPR may be observed even when PfEIR cannot be detected by current methods, and, therefore, both interventions and monitoring need to continue for some time after PfEIR has not been detected. Our results also highlight the importance of monitoring additional metrics of parasite transmission (in addition to PfEIR) when the annual PfEIR is <10 ib/person/year, especially when expecting a reduction in transmission as in the case of testing malaria interventions. Nonetheless, when PfEIR is above the level of detection, it provides information about the vector species involved in transmission, which is critical because different mosquito species may respond differently to vector control interventions (Ferguson et al., 2010; Wilson et al., 2020).
We observed a consistently higher PfPR in children 0.5–5 y/o than in women 15–49 y/o throughout the study region and study period, as expected. The extent of difference in PfPR between children and adults for a given region generally increases with parasite transmission intensity. However, even in mesoendemic settings (PfPR between 10–50%), it is common for PfPR in children to be appreciably higher than in adults (Smith et al., 2007a). This pattern is due to increasing acquired immunity with increased exposure to malaria parasites over time (Baird, 1995), which may decrease transmission efficiency and time to clear a P. falciparum infection in adults compared to children (see Appendix 1 – Procedure for building the HBR, PfSR, and PfPR models, Appendix 1—table 6 and Smith et al., 2005). Moreover, the higher PfPR in children than adults, even at the lowest levels of transmission, suggests that children may play a more significant role in transmission, consistent with other studies (Walldorf et al., 2015; Ouédraogo et al., 2016).
Although the functional form of our best fitting PfEIRPfPR model matches that of previous studies (Beier et al., 1999; Smith et al., 2005), the estimated values of PfPR as a function of PfEIR show nonnegligible differences. For example, based on our best model, a 0.15 ib/person/month decrease in PfEIR from 0.2 to 0.05 ib/person/month leads to a decrease in PfPR from 24.10% to 21.66% on average, a reduction of about 10.13%. In Beier et al., 1999, instead, the same decrease in PfEIR corresponds to a decrease in PfPR from 33.88% to 19.31%, a reduction of about 43.00%, whilst in the case of the best model of Smith et al., 2005, that yields a decrease in PfPR from 33.46% to 16.51%, a reduction of about 50.66%. One possible reason for these differences is that our study focuses on a geographically limited area where lower values in PfEIR are observed. Secondly, both previous studies excluded data from sites reporting mosquito control activities whilst our study site included multiple interventions. Finally, another important difference with our study is that, our focus was on a relatively small subnational area of Malawi, whereas Beier et al., 1999 and Smith et al., 2005 used data from across Africa and implicitly assumed that the properties of the PfEIRPfPR relationship do not vary over such continental scale.
In this study, the empirical PfPR was used as the response variables of six different statistical models, while the modelled PfEIR from a geostatistical model was used as the independent variable in each of the six models. As shown in this paper, this approach has two main advantages: (1) it allowed us to develop a bootstrap procedure for propagating the uncertainty arising from the estimates of PfEIR into the PfPRPfEIR relationship; (2) it allowed us to avoid the generation of spurious correlations in the estimation of the PfPRPfEIR relationship. The risk of spurious correlation may in fact occur when using approaches that are based on the estimates of both PfPR and PfEIR, which are obtained from statistical models informed by the same set of covariates.
One limitation of this study was that the six PfEIRPfPR models do not allow for overdispersion in the estimation of the PfEIRPfPR relationship. However, the use of standard Binomial likelihoods still delivers unbiased estimates of the functional relationship between PfEIR and PfPR, even in the presence of overdispersion (Godambe and Kale, 1991). Furthermore, given that the uncertainty around the PfPRPfEIR relationship is mainly driven by the predictive distribution of PfEIR, which we account for through our bootstrap procedure, we do not expect overdispersion to have nonnegigible influence on the parameter estimates. Finally, the development of models that allow for overdispersion may be achieved in several different ways, for example, by modelling the parameters that modulate the PfEIRPfPR relationship as stochastic processes. However, these approaches would require a larger amount of data than those available in this study and should be the subject of future research.
A second limitation was the use of RDTs to estimate PfPR. RDTs can show false positives after antimalarial drug treatment due to persistence of the antigens detected by RDTs (Dalrymple et al., 2018). Also, the limit of detection (usually 50–200 parasites/µl) is higher than that of expert microscopy or PCR (Chiodini, 2014). In modelling the relationship between PfEIR and PfPR, we did not account for the sensitivity and specificity of the RDT used to detect P falciparum infection. If the sensitivity α and specificity β were known, we could account for them by setting $PfPR(x,t)$ as used in our analysis to $\alpha (PfPR(x,t)\beta 1)/(\alpha +\beta 1)$. Thus, strictly, what we have called PfPR should be interpreted as the probability of testing positive for P. falciparum using RDT. However, the use of RDTs as a diagnostic test for the detection of malaria infection provides PfPR estimates that are comparable to national malaria indicator surveys.
An additional limitation of our study was the unidirectional relationship implicitly assumed in our models of PfEIRPfPR. In reality, PfPR and PfEIR are causally linked by the malaria parasite transmission cycle, which alternates between the human host and the mosquito vector. A higher rate of infectious bites received per person (i.e. PfEIR) increases the probability of the person becoming infected when bitten. Therefore, any factor that reduces mosquito populations, biting rates or humantomosquito parasite transmission (e.g. effective vector control interventions) will reduce PfEIR and consequently translate to reductions in PfPR. Similarly, a higher rate of parasite infection in people (i.e. PfPR) increases the probability of a mosquito becoming infected after any given blood meal. Therefore, factors that directly reduce PfPR (e.g. treatment of infections with effective drugs) will consequently reduce PfEIR. The impact of interventions may therefore affect both PfEIR and PfPR in such a way that a cyclic relationship may better describe the association between these metrics. Future modelling efforts may thus be improved by taking into account the cyclic aspect of the PfEIRPfPR relationship.
Conclusion
Measuring PfEIR and PfPR using the rolling MIS sampling framework and a geostatistical modelling approach allowed us to assess the finescale spatial and temporal distributions of malaria parasite transmission over 38 months in a mesoendemic setting. The relationship between PfEIR and PfPR estimated here shows that at low levels of PfEIR, changes in PfEIR are associated with rapid changes in PfPR, while at higher levels of PfEIR, changes in PfEIR are not associated with appreciable changes in PfPR. Comparing hotspots of PfEIR and PfPR revealed that each metric could identify potential transmission hotspots that the other fails to capture. Our results emphasise that PfEIR and PfPR are essential, complementary metrics for monitoring short term changes in P. falciparum transmission intensity in mesoendemic settings, which have become increasingly common as many regions reduce transmission and shift from the highest malaria endemicity levels. Our study emphasises the need to couple vector control with identifying and treating infected individuals to drive malaria to elimination levels and to monitor both entomological and parasitaemia indices in malaria surveillance.
Appendix 1
Procedure for building the HBR, PfSR and PfPR models
Let $Avg(\text{Temp}({x}_{i},{t}_{i}),{s}_{1},{s}_{2})$ and $Avg(\text{RH}({x}_{i},{t}_{i}),{s}_{1},{s}_{2})$ respectively denote the average temperature and relative humidity taken over s_{1} to s_{2} days prior to the data collection. Procedure for building the HBR, PfSR and PfPR models Appendix 1—table 1 shows the s_{1} and s_{2} values over which average temperature and relative humidity were computed. A set of these variables were selected as the best predictors each of the outcome variables based on the procedure in the next section.
We selected the best combination of fixed and random effects that best explain HBR, PfSR and PfPR using the following procedure.
We first built a generalized linear model in which temperature and RH are considered together with time trends and sine and cosine functions for seasonality. For $Avg(\text{Temp}({x}_{i},{t}_{i}),{s}_{1},{s}_{2})$, $Avg(\text{RH}({x}_{i},{t}_{i}),{s}_{1},{s}_{2})$, the choice of s_{1} and s_{2}, as illustrated by Procedure for building the HBR, PfSR and PfPR models Appendix 1—table 1, was based on the deviance profile of the variable involved, that is, either temperature or RH. Piecewiselinear transformations of temperature and RH were considered based on visual inspection and epidemiological knowledge.
Potential confounding between seasonal sinusoids, temperature and RH were checked. Covariates that did not improve the model fit as judged by the AIC were excluded. Sincosine terms were always considered together as if they were one covariate.
With the current model as a basic model we include other available explanatory variables based on forward selection, and check for interactions.
When no more explanatory variables significantly improve the model fit, we fit a generalized linear mixed model with a random effect for each unique spacetime location.
We then check for the presence of residual spatial, temporal, and spatiotemporal correlations using the algorithm described in Giorgi et al., 2018, and then include the random effect terms that improve the model fit.
The selected fixed effects for the HBR, PfSR and PfPR models
We specify the set of fixed effects we selected to be in the final model for the A. arabiensis HBR, A. funestus s.s. HBR, PfSR, and the PfPR models. Detailed description of the terms involved in the fixed effects and the estimates of all the parameters of each model are given in S1 Appendix 1—table 2–5.
A. arabiensis human biting rate
$\begin{array}{ll}{\displaystyle d({x}_{i},{t}_{i}{)}^{\mathrm{\top}}\beta +f({t}_{i};\alpha )=}& {\displaystyle {\beta}_{1}\mathbf{1}({x}_{i}\in \mathcal{\mathcal{A}})+{\beta}_{2}\mathbf{1}({x}_{i}\in \mathcal{\mathcal{B}})+{\beta}_{3}\mathbf{1}({x}_{i}\in \mathcal{\mathcal{C}})+{\beta}_{4}\mathbf{1}(\text{Indoor})+}\\ & {\displaystyle {\beta}_{5}\text{DSR}({x}_{i})+{\beta}_{6}Avg(\text{RH}({x}_{i},{t}_{i}),14,35)+}\\ & {\displaystyle {\beta}_{7}min\{Avg(\text{Temp}({x}_{i},{t}_{i}),7,14),22.9\}}\\ & {\displaystyle +{\beta}_{8}max\{Avg(\text{Temp}({x}_{i},{t}_{i}),7,14)22.9,0\}+}\\ & {\displaystyle {\alpha}_{1}\mathrm{sin}(2\pi {t}_{i}/12)/t+{\alpha}_{2}\mathrm{cos}(2\pi {t}_{i}/12)/t}\end{array}$A. funestus s.s. human biting rate
$\begin{array}{ll}{\displaystyle d({x}_{i},{t}_{i}{)}^{\mathrm{\top}}\beta +f({t}_{i};\alpha )=}& {\displaystyle {\beta}_{0}+\text{Elevation}({x}_{i})+{\beta}_{1}\text{DSR}({x}_{i})+{\beta}_{2}\text{NDVI}({x}_{i})+}\\ & {\displaystyle {\beta}_{3}Avg(\text{Temp}({x}_{i},{t}_{i}),0,7)+{\beta}_{4}Avg(\text{Temp}({x}_{i},{t}_{i}),7,14)+}\\ & {\displaystyle +{\beta}_{5}Avg(\text{RH}({x}_{i},{t}_{i}),14,21)}\\ & {\displaystyle +{\alpha}_{1}\mathrm{sin}(2\pi {t}_{i}/12)+{\alpha}_{2}\mathrm{cos}(2\pi {t}_{i}/12)+}\\ & {\displaystyle +{\alpha}_{3}min\{{t}_{i},12\}+{\alpha}_{4}max\{{t}_{i}12,0\}}\end{array}$A. arabiensis sporozoite rate
$\begin{array}{ll}{\displaystyle d({x}_{i},{t}_{i}{)}^{\mathrm{\top}}{\beta}^{\ast}+{f}^{\ast}({t}_{i};{\alpha}^{\ast})=}& {\displaystyle {\beta}_{0}^{\ast}+{\beta}_{1}^{\ast}\text{DLR}({x}_{i})+{\beta}_{2}^{\ast}\text{DSR}({x}_{i})+{\beta}_{3}^{\ast}\text{Elevation}({x}_{i})+}\\ & {\displaystyle {\beta}_{4}^{\ast}\text{EVI}({x}_{i})+{\alpha}_{1}^{\ast}\mathrm{sin}(2\pi {t}_{i}/12)+{\alpha}_{2}^{\ast}\mathrm{cos}(2\pi {t}_{i}/12)+}\\ & {\displaystyle {\alpha}_{3}^{\ast}min\{{t}_{i},12\}+{\alpha}_{4}^{\ast}max\{{t}_{i}12,0\}}\end{array}$A. funestus s.s. sporozoite rate
$\begin{array}{ll}{\displaystyle d({x}_{i},{t}_{i}{)}^{\mathrm{\top}}{\beta}^{\ast}+{f}^{\ast}({t}_{i};{\alpha}^{\ast})=}& {\displaystyle {\beta}_{0}^{\ast}+{\alpha}_{1}^{\ast}\mathrm{sin}(2\pi {t}_{i}/12)+{\alpha}_{2}^{\ast}\mathrm{cos}(2\pi {t}_{i}/12)+}\\ & {\displaystyle {\alpha}_{3}^{\ast}min\{{t}_{i},12\}+{\alpha}_{4}^{\ast}max\{{t}_{i}12,0\}}\end{array}$P. faciparum prevalence
$\begin{array}{ll}{\displaystyle d({x}_{i},{t}_{i}{)}^{\mathrm{\top}}\phi +g({t}_{i};\varrho )=}& {\displaystyle {\phi}_{1}\mathbf{1}({x}_{i}\in \mathcal{\mathcal{A}})+{\phi}_{2}\mathbf{1}({x}_{i}\in \mathcal{\mathcal{B}})+{\phi}_{3}\mathbf{1}({x}_{i}\in \mathcal{\mathcal{C}})+}\\ & {\displaystyle {\phi}_{4}\text{Elevation}({x}_{i})+{\phi}_{5}\text{DLR}({x}_{i})+}\\ & {\displaystyle {\phi}_{6}Avg(\text{Temp}({x}_{i},{t}_{i}),14,42)+{\phi}_{7}\text{NDVI}({x}_{i})+}\\ & {\displaystyle {\phi}_{8}\text{Wealth}({x}_{i})+{\varrho}_{1}min\{{t}_{i},21\}+{\varrho}_{2}max\{{t}_{i}21,0\}+}\\ & {\displaystyle +{\varrho}_{3}\mathrm{cos}(2\pi {t}_{i}/12)+{\varrho}_{4}\mathrm{sin}(2\pi {t}_{i}/12)}\end{array}$
Data availability
The data are available at Lancaster University Library online (https://doi.org/10.17635/lancaster/researchdata/479). Geolocations of individuals have been removed to anonymise the data. However, geolocation can be assessed by sending proposals to Professor Michèle van Vugt (m.vanvugt@amsterdamumc.nl). Data requestors will need to apply for ethical clearance and sign a data access agreement.

Lancaster University Research DataEntomology and parasitaemia data.https://doi.org/10.17635/lancaster/researchdata/479
References

Spatiotemporal modeling of sparse geostatistical malaria sporozoite rate data using a zero inflated binomial modelSpatial and SpatioTemporal Epidemiology 2:283–290.https://doi.org/10.1016/j.sste.2011.08.001

BookThe population dynamics of malariaIn: Aron J. L, editors. The Population Dynamics of Infectious Diseases: Theory and Applications. Springer. pp. 139–179.https://doi.org/10.1007/9781489929013_5

Short report: entomologic inoculation rates and Plasmodium falciparum malaria prevalence in africaThe American Journal of Tropical Medicine and Hygiene 61:109–113.https://doi.org/10.4269/ajtmh.1999.61.109

Asymptomatic malaria infections: detectability, transmissibility and public health relevanceNature Reviews Microbiology 12:833–840.https://doi.org/10.1038/nrmicro3364

Spatial targeting of interventions against malariaBulletin of the World Health Organization 78:1401–1411.

Malaria diagnostics: now and the futureParasitology 141:1873–1879.https://doi.org/10.1017/S0031182014001371

Adaptive geostatistical design and analysis for prevalence surveysSpatial Statistics 15:70–84.https://doi.org/10.1016/j.spasta.2015.12.004

Monte carlo maximum likelihood in ModelBased geostatisticsJournal of Computational and Graphical Statistics 13:702–718.https://doi.org/10.1198/106186004X2525

Humantomosquito transmission efficiency increases as malaria is controlledNature Communications 6:6054.https://doi.org/10.1038/ncomms7054

The effect of temperature on life history traits of culex mosquitoesJournal of Medical Entomology 51:55–62.https://doi.org/10.1603/me13003

A malaria model tested in the african savannaBulletin of the World Health Organization 50:347–357.

Acquired immunity to malariaClinical Microbiology Reviews 22:13–36.https://doi.org/10.1128/CMR.0002508

BookUnderstanding the Importance of Asymptomatic and Low Density Infections for Malaria EliminationIn: Manguin S, Dev V, editors. Towards Malaria Eliminationa Leap Forward. InTech. pp. 1–8.https://doi.org/10.5772/intechopen.77293

BookA Supplement to the Anophelinae of Africa South of the SaharaPublications of the South African Institute for Medical Research.

Geostatistical methods for disease mapping and visualisation using data from Spatiotemporally referenced prevalence surveysInternational Statistical Review 86:571–597.https://doi.org/10.1111/insr.12268

Prevmap: an r package for prevalence mappingJournal of Statistical Software 78:1–29.https://doi.org/10.18637/jss.v078.i08

BookEstimating functions: An overviewIn: Godambe V. P, editors. Estimating Functions. New York: Oxford Science Publications. pp. 3–20.

The microepidemiology of malaria and its importance to malaria controlTransactions of the Royal Society of Tropical Medicine and Hygiene 83 Suppl:25–29.https://doi.org/10.1016/00359203(89)905993

Measuring malaria endemicity from intense to interrupted transmissionThe Lancet Infectious Diseases 8:369–378.https://doi.org/10.1016/S14733099(08)700690

Trying to remember: immunological B cell memory to malariaTrends in Parasitology 31:89–94.https://doi.org/10.1016/j.pt.2014.12.009

Correlation of high levels of antibodies to multiple preerythrocytic Plasmodium falciparum antigens and protection from infectionThe American Journal of Tropical Medicine and Hygiene 73:222–228.https://doi.org/10.4269/ajtmh.2005.73.222

Climate change in semiarid Malawi: perceptions, adaptation strategies and water governanceJàmbá: Journal of Disaster Risk Studies 8:255.https://doi.org/10.4102/jamba.v8i3.255

BookSporozoite challenge and transmission patterns as determinants of occurrence of severe malaria in residents of kilifi district, kenyaNairobi: University of Nairobi.

A simplified model for predicting malaria entomologic inoculation rates based on entomologic and parasitologic parameters relevant to controlThe American Journal of Tropical Medicine and Hygiene 62:535–544.https://doi.org/10.4269/ajtmh.2000.62.535

Differential attractiveness of isolated humans to mosquitoes in TanzaniaTransactions of the Royal Society of Tropical Medicine and Hygiene 89:604–606.https://doi.org/10.1016/00359203(95)904069

A cocktail polymerase chain reaction assay to identify members of the anopheles funestus (Diptera: culicidae) groupThe American Journal of Tropical Medicine and Hygiene 66:804–811.https://doi.org/10.4269/ajtmh.2002.66.804

Relationships between Plasmodium falciparum transmission by vector populations and the incidence of severe disease at nine sites on the kenyan coastThe American Journal of Tropical Medicine and Hygiene 52:201–206.https://doi.org/10.4269/ajtmh.1995.52.201

Assessing the efficacy of candidate mosquito repellents against the background of an attractive source that mimics a human hostMedical and Veterinary Entomology 28:407–413.https://doi.org/10.1111/mve.12061

A novel synthetic odorant blend for trapping of malaria and other african mosquito speciesJournal of Chemical Ecology 38:235–244.https://doi.org/10.1007/s1088601200888

Indicators for the forecasting of malaria epidemicsBulletin of the World Health Organization 58:91–98.

HRP2: transforming malaria diagnosis, but with caveatsTrends in Parasitology 36:112–126.https://doi.org/10.1016/j.pt.2019.12.004

Interindividual variation in the attractiveness of human odours to the malaria mosquito anopheles gambiae s. sMedical and Veterinary Entomology 20:280–287.https://doi.org/10.1111/j.13652915.2006.00627.x

Rolling malaria indicator surveys (rMIS): a potential districtlevel malaria monitoring and evaluation (M&E) tool for program managersThe American Journal of Tropical Medicine and Hygiene 86:96–98.https://doi.org/10.4269/ajtmh.2012.110397

On the delayed RossMacdonald model for malaria transmissionBulletin of Mathematical Biology 70:1098–1114.https://doi.org/10.1007/s115380079292z

Identification of single specimens of the anopheles gambiae complex by the polymerase chain reactionThe American Journal of Tropical Medicine and Hygiene 49:520–529.https://doi.org/10.4269/ajtmh.1993.49.520

ConferenceIEEE international symposium on geoscience and remote sensing (IGARSS)2011 IEEE International Geoscience and Remote Sensing Symposium.https://doi.org/10.1109/IGARSS18752.2011

The matola malaria project: a temporal and spatial study of malaria transmission and disease in a suburban area of Maputo, MozambiqueThe American Journal of Tropical Medicine and Hygiene 57:550–559.https://doi.org/10.4269/ajtmh.1997.57.550

On the control of malaria in Freetown, sierra leone; Plasmodium falciparum and anopheles gambiae in relation to malaria occurring in infantsAnnals of Tropical Medicine and Parasitology 41:380–407.https://doi.org/10.1080/00034983.1947.11685341

The importance of vector control for the control and elimination of vectorborne diseasesPLOS Neglected Tropical Diseases 14:e0007831.https://doi.org/10.1371/journal.pntd.0007831
Decision letter

Elena A LevashinaReviewing Editor; Max Planck Institute for Infection Biology, Germany

Dominique SoldatiFavreSenior Editor; University of Geneva, Switzerland

Chihhao HsiehReviewer; National Taiwan University, Taiwan

Patrick BrownReviewer; University of Toronto, Canada
Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.
Acceptance summary:
Using an impressive dataset and a suit of advanced statistical models, the authors illustrated the spatialtemporal patterns of parasite prevalence (PR) and the entomological inoculation rate (EIR) in a rural area in Malawi. This research approach highlights the importance of integrated control malaria strategies based on both entomological and epidemiological surveys.
Decision letter after peer review:
Thank you for submitting your article "Identifying Plasmodium falciparum transmission patterns through parasite prevalence and entomological inoculation rate" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Dominique SoldatiFavre as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Chihhao Hsieh (Reviewer #1); Patrick Brown (Reviewer #2).
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential revisions:
1. The total number of households under survey is M=2432. The total number of survey months = 38. However, if I understand correctly, not every household is surveyed in each month. However, according to the explanation in the Data section, I cannot understand how many households were surveyed for each month. The data structure needs to be illustrated better so that the readers can get better sense out of the Geostatistical Model section. One simple solution is to make animation that shows the sampling sites for each month. The authors can include the animation as Supplementary Information (can be included in the web portal they provided).
2. How to deal with the uncertainty of PfEIR and PfPR in later analyses? Particularly, error propagates through HBRf (x, t)*PfSRf(x, t)l(t) and PfEIR(x,t) = PfEIRf(x,t) + PfEIRa(x,t). Several estimations have this error propagation issue. I am sure that authors can do these properly using bootstrap. Please explain the procedure clearly and earlier in their Geostatistical analysis section. The bootstrapped approach described in the later section is a good solution.
3. Please provide a clear demonstration of the extra value brought by the spatiotemporal model in terms of: (i) space – how are the spatial of the aspects of the association better characterized; and (ii) time – how does the relationship between PfEIR and PfPR vary dynamically over time.
4. On page 6, it would be worth mentioning that frequentist inference was used to estimate model parameters, many readers might associate geostatistical models as Bayesian.
5. The formula on p.16 for calculating PfPR while accounting for sensitivity and specificity seems to be wrong. Shouldn't it be (PfPR + β – 1)/(α + β – 1)? If the formula in the manuscript is correct, please provide a reference.
6. Figure 1: Attribution is required for the map figure (are the raster tiles from Google, OpenStreetMap, or other sources?)
7. Figure 4: is the confidence region calculated with error propagation of estimates of PfSR and PfEIR? Likewise, the uncertainty associated with estimates needs to be carefully addressed in all estimates and clearly depicted in the figures.
8. While the online visualisation is good, the manuscript could include static map figures of the results (PfEIR and PfPR). For example, select certain time points with combination of high/low PfEIR and high/low PfPR. This is useful for archive purposes, and also because only the maps in the manuscript can be said to have been peerreviewed.
Reviewer #1 (Recommendations for the authors):
1. The total number of households under survey is M=2432. The total number of survey months = 38. However, if I understand correctly, not every household is surveyed in each month. However, according to the explanation in the Data section, I cannot understand how many households were surveyed for each month. The data structure needs to be illustrated better so that the readers can get better sense out of the Geostatistical Model section. One simple solution is to make animation that shows the sampling sites for each month. The authors can include the animation as Supplementary Information (can be included in the web portal they provided).
2. How to deal with the uncertainty of PfEIR and PfPR in later analyses? Particularly, error propagates through HBRf (x, t)*PfSRf(x, t)l(t) and PfEIR(x,t) = PfEIRf(x,t) + PfEIRa(x,t). Several estimations have this error propagation issue. I am sure that authors can do these properly using bootstrap. Please explain the procedure clearly and earlier in their Geostatistical analysis section. The bootstrapped approach described in the later section is a good solution.
3. Figure 4: is the confidence region calculated with error propagation of estimates of PfSR and PfEIR? Likewise, the uncertainty associated with estimates needs to be carefully addressed in all estimates and clearly depicted in the figures.
Reviewer #2 (Recommendations for the authors):
On page 6, it would be worth mentioning that frequentist inference was used to estimate model parameters, many readers might associate geostatistical models as Bayesian.
The formula on p.16 for calculating PfPR while accounting for sensitivity and specificity seems to be wrong. Shouldn't it be (PfPR + β – 1)/(α + β – 1)? If the formula in the manuscript is correct, please provide a reference.
Figures and Tables:
Figure 1: Attribution is required for the map figure (are the raster tiles from Google, OpenStreetMap, or other sources?)
While the online visualisation is good, the manuscript could include static map figures of the results (PfEIR and PfPR). For example, select certain time points with combination of high/low PfEIR and high/low PfPR. This is useful for archive purposes, and also because only the maps in the manuscript can be said to have been peerreviewed.
Reviewer #3 (Recommendations for the authors):
To strengthen this work, I would want a clear demonstration of the extra value brought by the spatiotemporal model in terms of: (i) space – how are the spatial of the aspects of the association better characterized; and (ii) time – how does the relationship between PfEIR and PfPR vary dynamically over time.
[Editors' note: further revisions were suggested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Identifying Plasmodium falciparum transmission patterns through parasite prevalence and entomological inoculation rate" for further consideration by eLife. Your revised article has been evaluated by Dominique SoldatiFavre (Senior Editor) and a Reviewing Editor.
The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:
1. Highlight and clarify in the main text the issue related to potential spurious correlation (statistical artifact) between the modelled PfPR and PfEIR.
2. Address in the main text how the inclusion of an over dispersion term will change the results of modelling the relationship between PfPR and PfEIR in the six models.
3. Introduce into the main text comparison of the estimation PfEIR vs PfPR relationship with previously published models from the literature and explain the observed differences beyond the differences between the modelling approaches.
https://doi.org/10.7554/eLife.65682.sa1Author response
Essential revisions:
1. The total number of households under survey is M=2432. The total number of survey months = 38. However, if I understand correctly, not every household is surveyed in each month. However, according to the explanation in the Data section, I cannot understand how many households were surveyed for each month. The data structure needs to be illustrated better so that the readers can get better sense out of the Geostatistical Model section. One simple solution is to make animation that shows the sampling sites for each month. The authors can include the animation as Supplementary Information (can be included in the web portal they provided).
We have included two plots as Figure 2 in the Results section (rMIS and mosquito sampling subsection), which indicate how the number of rMIS participants and houses for mosquito samplings were distributed over the course of the study. We have also added a sentence to the Data subsection of Methods as further clarification.
2. How to deal with the uncertainty of PfEIR and PfPR in later analyses? Particularly, error propagates through HBRf (x, t)*PfSRf(x, t)l(t) and PfEIR(x,t) = PfEIRf(x,t) + PfEIRa(x,t). Several estimations have this error propagation issue. I am sure that authors can do these properly using bootstrap. Please explain the procedure clearly and earlier in their Geostatistical analysis section. The bootstrapped approach described in the later section is a good solution.
We used bootstrapping as the reviewer suggests. We have added text to the Estimating the Plasmodium falciparum entomological inoculation rate section to provide this further clarification and explanation.
3. Please provide a clear demonstration of the extra value brought by the spatiotemporal model in terms of: (i) space – how are the spatial of the aspects of the association better characterized; and (ii) time – how does the relationship between PfEIR and PfPR vary dynamically over time.
Thanks to the reviewer for raising this point, which we hope will clarify our assumptions about the PfEIRPfPR relationship and highlight the contributions of this paper. We have made a number of significant changes to the Discussion section to address these points, most notably paragraphs number 24 and 6.
4. On page 6, it would be worth mentioning that frequentist inference was used to estimate model parameters, many readers might associate geostatistical models as Bayesian.
We have clarified that parameters were estimated using Monte Carlo Maximum Likelihood and added relevant references.
5. The formula on p.16 for calculating PfPR while accounting for sensitivity and specificity seems to be wrong. Shouldn't it be (PfPR + β – 1)/(α + β – 1)? If the formula in the manuscript is correct, please provide a reference.
The reviewer is correct. This was a mistake on our part, and we have corrected this as suggested.
6. Figure 1: Attribution is required for the map figure (are the raster tiles from Google, OpenStreetMap, or other sources?)
Done.
7. Figure 4: is the confidence region calculated with error propagation of estimates of PfSR and PfEIR? Likewise, the uncertainty associated with estimates needs to be carefully addressed in all estimates and clearly depicted in the figures.
All confidence regions shown in all plots were obtained using bootstrapping. We have now clarified this.
8. While the online visualisation is good, the manuscript could include static map figures of the results (PfEIR and PfPR). For example, select certain time points with combination of high/low PfEIR and high/low PfPR. This is useful for archive purposes, and also because only the maps in the manuscript can be said to have been peerreviewed.
We have selected three time points representing high, low and medium transmission and have included static maps thereof in the paper. We have highlighted the selected time points in the main text.[Editors' note: further revisions were suggested prior to acceptance, as described below.]
The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:
1. Highlight and clarify in the main text the issue related to potential spurious correlation (statistical artifact) between the modelled PfPR and PfEIR.
We have included a paragraph discussing this issue in the Discussion section.
2. Address in the main text how the inclusion of an over dispersion term will change the results of modelling the relationship between PfPR and PfEIR in the six models.
We have included a paragraph discussing this issue in the Discussion section.
3. Introduce into the main text comparison of the estimation PfEIR vs PfPR relationship with previously published models from the literature and explain the observed differences beyond the differences between the modelling approaches.
We have introduced a paragraph in the Discussion to make this comparison.
https://doi.org/10.7554/eLife.65682.sa2Article and author information
Author details
Funding
Stichting Dioraphte (13050800)
 Willem Takken
National Institutes of Health (T32AI007524)
 Robert S McCann
National Institutes of Health (K01TW011770)
 Robert S McCann
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This study was generously supported by Dioraphte Foundation, The Netherlands. RSM received additional support from the NIH (T32AI007524 and K01TW011770). The content is solely the responsibility of the authors and does not necessarily represent the official views of the funders. We thank African Parks and The Hunger Project for their significant and practical contributions in facilitating the study. We are grateful to the entire Majete Malaria Project team for their tireless efforts in carrying out the study. The population of the study area is thanked for their cooperation with the study. We also thank: Alexandra Hiscox for advice on mosquito sampling and study design; Jeroen Spitzen for logistical assistance; and Martin Donnelly, Karl Seydel, and their respective laboratory teams for assistance in molecular identification of malaria parasites and anopheline mosquitoes.
Ethics
Human subjects: Written informed consent was sought from individual participants (and their parents or guardians when appropriate) during household surveys. The College of Medicine Research and Ethics Committee in Malawi approved this study (proposal number P.05/15/1731).
Senior Editor
 Dominique SoldatiFavre, University of Geneva, Switzerland
Reviewing Editor
 Elena A Levashina, Max Planck Institute for Infection Biology, Germany
Reviewers
 Chihhao Hsieh, National Taiwan University, Taiwan
 Patrick Brown, University of Toronto, Canada
Publication history
 Received: December 11, 2020
 Preprint posted: January 16, 2021 (view preprint)
 Accepted: September 30, 2021
 Version of Record published: October 21, 2021 (version 1)
Copyright
© 2021, Amoah et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,071
 Page views

 148
 Downloads

 6
 Citations
Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Epidemiology and Global Health
Background:
Shortterm forecasts of infectious disease burden can contribute to situational awareness and aid capacity planning. Based on best practice in other fields and recent insights in infectious disease epidemiology, one can maximise the predictive performance of such forecasts if multiple models are combined into an ensemble. Here, we report on the performance of ensembles in predicting COVID19 cases and deaths across Europe between 08 March 2021 and 07 March 2022.
Methods:
We used opensource tools to develop a public European COVID19 Forecast Hub. We invited groups globally to contribute weekly forecasts for COVID19 cases and deaths reported by a standardised source for 32 countries over the next 1–4 weeks. Teams submitted forecasts from March 2021 using standardised quantiles of the predictive distribution. Each week we created an ensemble forecast, where each predictive quantile was calculated as the equallyweighted average (initially the mean and then from 26th July the median) of all individual models’ predictive quantiles. We measured the performance of each model using the relative Weighted Interval Score (WIS), comparing models’ forecast accuracy relative to all other models. We retrospectively explored alternative methods for ensemble forecasts, including weighted averages based on models’ past predictive performance.
Results:
Over 52 weeks, we collected forecasts from 48 unique models. We evaluated 29 models’ forecast scores in comparison to the ensemble model. We found a weekly ensemble had a consistently strong performance across countries over time. Across all horizons and locations, the ensemble performed better on relative WIS than 83% of participating models’ forecasts of incident cases (with a total N=886 predictions from 23 unique models), and 91% of participating models’ forecasts of deaths (N=763 predictions from 20 models). Across a 1–4 week time horizon, ensemble performance declined with longer forecast periods when forecasting cases, but remained stable over 4 weeks for incident death forecasts. In every forecast across 32 countries, the ensemble outperformed most contributing models when forecasting either cases or deaths, frequently outperforming all of its individual component models. Among several choices of ensemble methods we found that the most influential and best choice was to use a median average of models instead of using the mean, regardless of methods of weighting component forecast models.
Conclusions:
Our results support the use of combining forecasts from individual models into an ensemble in order to improve predictive performance across epidemiological targets and populations during infectious disease epidemics. Our findings further suggest that median ensemble methods yield better predictive performance more than ones based on means. Our findings also highlight that forecast consumers should place more weight on incident death forecasts than incident case forecasts at forecast horizons greater than 2 weeks.
Funding:
AA, BH, BL, LWa, MMa, PP, SV funded by National Institutes of Health (NIH) Grant 1R01GM109718, NSF BIG DATA Grant IIS1633028, NSF Grant No.: OAC1916805, NSF Expeditions in Computing Grant CCF1918656, CCF1917819, NSF RAPID CNS2028004, NSF RAPID OAC2027541, US Centers for Disease Control and Prevention 75D30119C05935, a grant from Google, University of Virginia Strategic Investment Fund award number SIF160, Defense Threat Reduction Agency (DTRA) under Contract No. HDTRA119D0007, and respectively Virginia Dept of Health Grant VDH215010141, VDH215010143, VDH215010147, VDH215010145, VDH215010146, VDH215010142, VDH215010148. AF, AMa, GL funded by SMIGE  Modelli statistici inferenziali per governare l'epidemia, FISR 2020Covid19 I Fase, FISR2020IP00156, Codice Progetto: PRJ0695. AM, BK, FD, FR, JK, JN, JZ, KN, MG, MR, MS, RB funded by Ministry of Science and Higher Education of Poland with grant 28/WFSN/2021 to the University of Warsaw. BRe, CPe, JLAz funded by Ministerio de Sanidad/ISCIII. BT, PG funded by PERISCOPE European H2020 project, contract number 101016233. CP, DL, EA, MC, SA funded by European Commission  DirectorateGeneral for Communications Networks, Content and Technology through the contract LC01485746, and Ministerio de Ciencia, Innovacion y Universidades and FEDER, with the project PGC2018095456BI00. DE., MGu funded by Spanish Ministry of Health / REACTUE (FEDER). DO, GF, IMi, LC funded by Laboratory Directed Research and Development program of Los Alamos National Laboratory (LANL) under project number 20200700ER. DS, ELR, GG, NGR, NW, YW funded by National Institutes of General Medical Sciences (R35GM119582; the content is solely the responsibility of the authors and does not necessarily represent the official views of NIGMS or the National Institutes of Health). FB, FP funded by InPresa, Lombardy Region, Italy. HG, KS funded by European Centre for Disease Prevention and Control. IV funded by Agencia de Qualitat i Avaluacio Sanitaries de Catalunya (AQuAS) through contract 2021021OE. JDe, SMo, VP funded by Netzwerk Universitatsmedizin (NUM) project egePan (01KX2021). JPB, SH, TH funded by Federal Ministry of Education and Research (BMBF; grant 05M18SIA). KH, MSc, YKh funded by Project SaxoCOV, funded by the German Free State of Saxony. Presentation of data, model results and simulations also funded by the NFDI4Health Task Force COVID19 (https://www.nfdi4health.de/taskforcecovid192) within the framework of a DFGproject (LO342/171). LP, VE funded by Mathematical and Statistical modelling project (MUNI/A/1615/2020), Online platform for realtime monitoring, analysis and management of epidemic situations (MUNI/11/02202001/2020); VE also supported by RECETOX research infrastructure (Ministry of Education, Youth and Sports of the Czech Republic: LM2018121), the CETOCOEN EXCELLENCE (CZ.02.1.01/0.0/0.0/17043/0009632), RECETOX RI project (CZ.02.1.01/0.0/0.0/16013/0001761). NIB funded by Health Protection Research Unit (grant code NIHR200908). SAb, SF funded by Wellcome Trust (210758/Z/18/Z).

 Epidemiology and Global Health
Background:
Affectionate touch, which is vital for mental and physical health, was restricted during the Covid19 pandemic. This study investigated the association between momentary affectionate touch and subjective wellbeing, as well as salivary oxytocin and cortisol in everyday life during the pandemic.
Methods:
In the first step, we measured anxiety and depression symptoms, loneliness and attitudes toward social touch in a large crosssectional online survey (N = 1050). From this sample, N = 247 participants completed ecological momentary assessments over 2 days with six daily assessments by answering smartphonebased questions on affectionate touch and momentary mental state, and providing concomitant saliva samples for cortisol and oxytocin assessment.
Results:
Multilevel models showed that on a withinperson level, affectionate touch was associated with decreased selfreported anxiety, general burden, stress, and increased oxytocin levels. On a betweenperson level, affectionate touch was associated with decreased cortisol levels and higher happiness. Moreover, individuals with a positive attitude toward social touch experiencing loneliness reported more mental health problems.
Conclusions:
Our results suggest that affectionate touch is linked to higher endogenous oxytocin in times of pandemic and lockdown and might buffer stress on a subjective and hormonal level. These findings might have implications for preventing mental burden during social contact restrictions.
Funding:
The study was funded by the German Research Foundation, the German Psychological Society, and German Academic Exchange Service.