Seasonal Dynamics of Anopheles stephensi and its Implications for Mosquito Detection and Emergent Malaria Control in the Horn of Africa

Systematic Review of Anopheles stephensi Surveys

We collated references from published systematic reviews of literature relating to An. stephensi^16,35, and updated these previous searches by searching Web of Science and PubMed from Jan.2017 to Sep.2020. We included all records containing temporally disaggregated adult mosquito catch data with monthly (or finer) temporal resolution spanning at least 10 months, that had not been conducted as part of vector control intervention trials and where at least 25 An. stephensi mosquitoes had been caught over the study period. A total of 36 references were collated containing 65 time-series with monthly catch data (no study presented data at a finer temporal resolution) from surveys carried out across Afghanistan, Djibouti, India, Iran, Myanmar and Pakistan. See Supplementary Information for further details and references therein.

Clustering of Similar Time-Series & Random Forest Prediction of Cluster Membership

Following methodologies developed in previous work³⁵, we fitted a Gaussian Process-based model to smooth these mosquito count time-series, using a Negative Binomial likelihood to account for overdispersion and a periodic kernel function to capture the repeating patterns often observed seasonally in mosquito populations. Model fitting was carried out within a Bayesian framework, using the probabilistic programming language STAN³⁶. We then calculated summary statistics for each smoothed time-series to characterise their temporal properties (Supplementary Information), generating a set of parameters for each time-series that summarises their temporal properties. We then scaled and normalised each summary statistic to give a mean 0 and unit variance –a process necessary for the principal components analysis (PCA) we apply to identify a lower-dimensional representation of the structure present in the data amenable to visualisation.

Using k-means clustering, we identified clusters of time-series with similar temporal properties – the output of this process is a label for each time-series indicating which cluster (of time-series with similar temporal properties) each specific time-series was assigned to. For each study location, we extracted a suite of satellite-derived environmental variables (Supplementary Table 2) and used these variables alongside empirically calculated rainfall seasonality and average monthly catch as covariates within a random-forest based classification framework to predict cluster membership of each time-series. These models were fitted using the R package Ranger³⁷ with 6-fold cross-validation utilised to optimise hyperparameters. Results are based on averaging the results of 25 iterations of cross-validation and model fitting and predictions made using out-of-bag estimates. There were significant imbalances in class size across clusters and so we carried out upsampling using the SMOTE³⁸ algorithm. For results without upsampling see Supplementary Information.

Probability of Detecting Anopheles stephensi With Different Surveillance Strategies

We explore the implications of seasonal variation in An. stephensi abundance on the probability of detecting the vector in entomological surveys using a theoretical sampling method with a defined amount of effort (such as a human landing catch conducted by a single volunteer for one night). We use a statistical framework (described further in Supplementary Information) that calculates the cumulative probability of detection from: i) an overall assumed An. stephensi annual biting rate (ABR, arbitrarily set to 20 for illustrative purposes here); ii) changes in vector density over the course of the year (from our collated time-series); and iii) various factors relating to timing of, and effort expended in, the entomological survey. Specifically, for each time-series, we identified the month with the highest rainfall, and the month in which vector density was highest (noting that these months were very rarely the same month). We then calculated the cumulative probability of An. stephensi detection using the framework, under a range of different surveillance strategies. Specifically, three strategies were simulated:

• Vector-Peak Timed: Starting the survey at the month with peak vector density (noting that in the absence of pre-existing detailed entomological information this is a hypothetical quantity designed to illustrate an approximate upper bound on the detection probability that could be achieved).

• Rainfall-Peak Timed: Starting the survey at the month with peak rainfall.

• Random Month Timed: The expected probability of detection achieved if the survey was started during a random month (calculated by simulating survey starting in each of the year’s 12 months and calculating the average cumulative probability).

In addition to varying the survey’s starting time, we also varied the amount of sampling effort (number of days sampled within each month) and overall duration of the survey (consecutive months sampled given a defined number of nights sampling per month). Note that the aim here is not to describe the exact probability of missing An. stephensi in any given entomological survey, as this will depend on a wide array of other, poorly defined and heterogeneous factors (e.g. type of catch methodology used, location etc). We also assume that the collection method is unbiased (i.e. not biased towards catching mosquitoes with particular resting or biting properties) which is also highly unlikely. Instead, the aim is to highlight how variation in seasonal dynamics can influence the nature of surveillance required to successfully detect a single An. stephensi (i.e. successfully establish presence), and the probability of detection should be viewed as a relative measure (i.e. viewed in relation to other sampling efforts and survey timings possible for surveys) and not an absolute value. Note that this framework assumes no seasonal variation in factors other than mosquito abundance (such as the ability of the sampling method to accurately record ABR) that might influence the probability of An. stephensi being caught by our theoretical sampling method.

Modelling Anopheles stephensi-Driven Malaria Dynamics and Control

We integrated these vector abundance time-series into a published population-level model of Plasmodium falciparum malaria transmission and disease dynamics^39–41 to explore the implications of An. stephensi seasonality on malaria control in settings in the Horn of Africa where malaria is currently largely absent (see Supplementary Information for full description of the modelling framework). We use the modelling framework to understand how variation in mosquito seasonality might influence the impact of IRS, a key vector control intervention. As the dynamics of An. stephensi’s establishment and influence on temporal trends in malaria transmission during its establishment remain unclear, we focused on the time-period immediately following establishment (when the disease is at equilibrium) and provide an illustrative example of how seasonality of An. stephensi driven malaria transmission could influence the effectiveness of IRS in a site with no pre-existing history of malaria control. For simplicity we assume that all transmission is due to An. stephensi and that IRS efficacy against this species is consistent with that observed against other species across the continent⁴². We simulate the impact of a single, illustrative IRS spraying campaign in a setting with an annual EIR of 1.5 (average malaria prevalence of 8-9%), timed either for optimal impact, randomly or alongside peak rainfall, and assume that 80% of the vector’s resting sites are successfully sprayed (noting the vector is thought to also rest in animal houses which are not typically sprayed in public-health campaigns). For further details, see Supplementary Information.

Diversity in Temporal Dynamics Across the Collated Anopheles stephensi Time-Series

A total of 65 time-series from across Afghanistan, Djibouti, India, Iran, Myanmar and Pakistan were identified (Fig.1A, Supp Fig.1). Substantial variation in the degree and timing of vector seasonality was observed, with the maximum percentage of annual vector density in any consecutive 4-month period (a proxy for degree of seasonality) ranging from 35-99% across the collated studies (average=62%). This contrasted with rainfall seasonality, where highly seasonal rainfall patterns were consistently observed across the locations the surveys had been carried out in (maximum percentage of annual rainfall in any consecutive 4-month period, mean=82, range 47-99%). We also observed a diverse range of temporal patterns, ranging from highly seasonal dynamics with a single seasonal peak, to bimodal population dynamics with two peaks within a single year, or more perennial patterns of abundance (Fig.1B).

• Download figure
• Open in new tab

Figure 1: Sources and Locations of Anopheles stephensi Time-Series Data and Examples for Each Country.

(A) Map of the geographical range over which time-series entomological collections have been carried out. Countries with studies are highlighted in light grey, and the locations of individual studies indicated by the individual points, coloured according to country (Afghanistan=red, Djibouti=yellow, India=green, Iran=turquoise, Myanmar=blue and Pakistan=pink). (B) Example An. stephensi time-series from each country, with the empirical monthly mosquito catch (black points), fitted gaussian process curves (mean=coloured line, ribbon=95% Bayesian Credible Interval) and monthly rainfall (matching sampling location and year of sampling) for each (light blue bars). The x-axis indicates the month of sampling, the y-axis either the monthly rainfall (left hand side y-axis) or number of vectors caught in each month (right hand side y-axis; note that the absolute number of mosquitoes caught between time-series are not comparable due to variable sampling effort). n indicates the number of time-series in each country.

Statistical Characterisation and Clustering of Temporal Properties Highlights Distinct Archetypes

Summary statistics were calculated for each time-series to characterise their temporal properties (Supp Fig.2), followed by k-means clustering of the results to cluster the time-series into groups with similar temporal patterns. Our results highlight two distinct clusters of time-series, each characterised by distinct temporal patterns (Fig.2B). Cluster 1 time-series had single seasonal peaks and were more seasonal (average of 68% of annual vector density in the consecutive 4-month period with highest density) than Cluster 2 time-series, which had more perennial patterns of annual abundance (average 44% of annual vector density in the consecutive 4-month period with highest density) and contained several time-series with two peaks across the year. Despite differing significantly in mean vector abundance seasonality (Fig.2C and 2D, p<0.001), there was no significant difference between Clusters in rainfall seasonality (Fig.2D, p=0.59). Seasonality of rainfall (defined as the highest proportion of total annual rainfall occurring in any consecutive 4 month period) across sampled locations was high (average 82% and 84% for Clusters 1 and 2 respectively) despite wide variation in vector abundance seasonality. Timing of peak rainfall relative to peak vector density significantly differed between clusters (Supp Fig.3), with peak rainfall and vector abundance separated by <1 month on average for Cluster 1 compared to 2.2 months for Cluster 2. There was, however, considerable within-cluster variation in timing–within Cluster 1, timing of peak vector density relative to rainfall ranged from -5.8 months to +5.3 months (with 6 months the maximum gap that can occur within an annually repeating 12 month time-series, highlighting that the peaks in vector density relative to rainfall were found across the entirety of the year). We also explored varying the number of clusters specified in the k-means algorithm. Specifying 4 clusters resulted in further disaggregation of the 49 time-series in Cluster 1 into 3 separate clusters, each characterised by a single seasonal peak, but which differed in the timing of peak vector-density relative to peak rainfall (Supp Fig.4).

• Download figure
• Open in new tab

Figure 2: Characterisation and Clustering to Identify Time-Series with Similar Temporal Properties.

(A) Results of principal components analysis (PCA) and k-means clustering for 2 clusters. Points on main figure indicate individual time-series, with point colour indicating cluster membership. Ellipsoids demarcate the 75^th quantile of the density associated with each cluster. Principal components 1 and 2 are plotted, together explaining 69% of the total variation in temporal properties across the time-series. (B) Time-series belonging to each cluster. Pale lines represent individual time-series, brighter line the mean of all the time-series belonging to that cluster – in all cases vector density is normalised to sum to 1 over the course of the year, and time-standardised so that the highest vector density for each time-series is arbitrarily set to occur at month 7. (C) Plot comparing the percentage of annual total mosquito catch and percentage of annual total rainfall occurring in any consecutive 4-month period for each time-series, coloured by cluster membership. (D) Boxplots show the percentage of annual total mosquito catch (left) and annual total rainfall (right) series occurring any in consecutive 4-month periods for each time-series. Rainfall data comes from the CHIRPS dataset³⁷ and is specific to study location and time-period. Each point indicates an individual time-series. (E) Variable importance plot for each of the covariates included in the random forest model used to predict cluster membership– bar height indicates the mean variable importance across the 25 individual iterations of random forest fitting, with error bars representing the 95% confidence interval. Inset plots are the partial dependence plots for the top 5 most important variables in the model showing how the average prediction for Cluster 2 (y-axis, with higher values indicating an increased probability of Cluster 2 membership) varies with (normalised) variable value (x-axis).

Random-Forest Modelling of Seasonal Dynamics Highlights Urbanicity as a Key Factor

We fitted a random forest-based classification framework to predict cluster membership (Cluster 1 or Cluster 2, as defined in Fig.2A). Due to the significant class size imbalance between Cluster 1 (n=49) and Cluster 2 (n=16), we up-sampled Cluster 2 data to generate balanced classes. Across 25 iterations of random forest model fitting, mean AUC was 0.89 (indicating good predictive performance) and the model was able to correctly classify Cluster 1 and Cluster 2 time-series equally well (83% and 85% accuracy respectively).

We calculated the relative importance of each variable to the model’s predictive ability (Fig.2E). Patterns of land-use were strongly associated with different clusters–time-series from surveys in locations with lower population density (a proxy for rurality) more likely to belong to Cluster 2 (less seasonal), as were areas with a high proportion of land occupied by irrigated cropland. By contrast, a high proportion of land occupied by rainfed cropland was associated with Cluster 1 (more seasonal) dynamics. We also observed strong associations with temperature covariates, including the mean temperature of the driest quarter (where a high temperature was associated with Cluster 2), temperature seasonality (where a non-monotonically increasing relationship was observed, see Fig.2E inset panels and Supp Fig.5 for all covariate response plots) and whether the study had been conducted in Iran (indicating potential spatial confounding). By contrast, rainfall seasonality was not an important predictor of temporal dynamics and was in the least 5 important variables. Examining the association between cluster membership and rurality/urbanicity (defined by study authors), there was indication of an association (chi-squared test, p=0.07), though this was not statistically significant at the 5% level. 88% (n=22/25) time-series from urban settings were assigned to Cluster 1, and only 12% (n=3/25) assigned to Cluster 2. 65% (n=24/37) of time-series from rural settings were assigned to Cluster 1, and 35% (n=13/37) to Cluster 2.

Model predictive performance and variable importance rankings were similar when no up-sampling was applied (AUC=0.81, Supp Fig.6), though predictive accuracy was highly unbalanced (Cluster 2 accuracy=50%, Cluster 1 accuracy=94%). Model performance and variable importance ordering remained similar when fitting the model and explicitly holding out a subset of the data to subsequently evaluate model performance (n=7 time-series, Supp Fig.7). Predictive power for seasonality (percentage of vector catch in any 4-month period) was more modest, although estimates were positively correlated (r=0.43, Supp Fig.8).

Implications of Seasonal Dynamics for Entomological Surveillance of Anopheles stephensi across the Horn of Africa

We collated the same covariates for countries across the Horn of Africa and used the random forest model to predict cluster membership and potential temporal dynamics of An. stephensi across the region (Fig.3A). Our results highlight distinct geographical areas considered more likely to belong to Cluster 1 (more seasonal) and Cluster 2 (less seasonal), as well as areas of significant uncertainty. We next asked what consequences this seasonality might have on entomological surveillance of the vector, with a focus on how these seasonal dynamics might interact with features of surveillance programmes such as the timing and duration of entomological surveys. Across the collated temporal profiles, in a setting with an ABR of 20, surveys consisting of 3 months sampling and 3 sampling days per month that were timed to start at periods of peak An. stephensi density were on average 64% more likely to detect the vector compared to starting the survey at a random month of the year; and 57% more likely to successfully detect the vector compared to starting the survey in the month of peak rainfall (Fig.3B). Timing entomological surveys to coincide with peaks in rainfall did not lead to a significant increase in the probability of successfully detecting An. stephensi (average 4% increase), suggesting that the timing of peak rainfall may be a poor measure for guiding entomological surveys searching for the vector. We next stratified these results by temporal Cluster (Fig.3C). For Cluster 1 (and a survey lasting 3 months, with 3 days sampling per month) we observed differences in the cumulative probability of detection when comparing strategies which start surveys at the location’s rainfall peak, compared to starting them at peak An. stephensi abundance – on average, the latter strategy increased the cumulative probability of detection by 62% compared to a randomly timed survey, compared to only a 40% increase over random timing for Cluster 2.

• Download figure
• Open in new tab

Figure 3: Possible Seasonal Dynamics of Anopheles stephensi Across the Horn of Africa, and Consequences for Entomological Surveillance and Monitoring.

(A) Environmental covariates were collated across countries in the Horn of Africa where An. stephensi has been found, and the random forest classification model from Fig 2E used to predict potential temporal dynamics. Map shows the probability of temporal dynamics belonging to Cluster 1 (more seasonal), with pink corresponding to Cluster 1 dynamics being more likely and black indicating Cluster 2 dynamics (more perennial) are more likely, with white indicating both are equally likely. (B) For a setting with an annual biting rate of 20, the average probability (across all 65 collated An. stephensi temporal profiles) of detecting An. stephensi (where detection is defined as catching ≥1 mosquito) for a range of different sampling efforts (number of consecutive months sampled and number of sampling days in each month) in a setting with an annual biting rate of 20 bites per person. These results were generated for 3 different sampling strategies: i) with sampling starting at a random month in the year (left hand panel, subsequently averaged over all possible sampling start months in the year); ii) with sampling starting in the month of peak rainfall (centre panel); or iii) with sampling starting in the month of peak vector density (right hand panel). (C) For setting with an annual biting rate of 20 and a sampling effort of 3 days per month, the cumulative probability of An. stephensi detection as a function of the number of consecutive months sampled for each individual time-series, stratified by sampling strategy (starting at peak vector abundance=orange, at peak rainfall=blue) and Cluster. In both panels, pale, thin lines indicate the cumulative probability curve for a specific temporal profile, and thicker lines indicate the average for the specific sampling strategy.

Modelling the Impact of Anopheles stephensi Seasonality On Vector Control Measures

Integrating the temporal profiles of An. stephensi abundance with a malaria transmission model, we explored how variation in temporal dynamics influences the impact of IRS (with two different insecticides, Fig.4A). Across the An. stephensi temporal profiles, optimal timing of IRS delivery resulted in an average of 47.6% reduction in annual malaria incidence in the 12 months following spraying for pirimithos methyl, and 28.9% for bendiocarb (Fig.4B). These results represent 1.12x and 1.41x increases over the average impact achieved if the campaign is timed to a random month of the year. The extent to which optimal timing provided greater impact than random timing was dependent on the degree of seasonality and insecticide however–it increased with the degree of seasonality and was consistently larger for bendiocarb than pirimithos methyl (due to the latter’s longer duration and retention of residual activity following spraying). Timing the IRS campaign to occur when rainfall peaks did not significantly increase impact compared to timing the campaign to a random month (with less than a 2% average increase in impact for both pirimithos methyl (Fig.4C) and bendiocarb (Fig.4D)) and had significantly lower impact than optimally timed campaigns (39% and 15% lower impact for pirimithos methyl and bendiocarb respectively).

• Download figure
• Open in new tab

Figure 4: Modelling the Public-Health Impact of Indoor Residual Spraying and the Influence of Anopheles stephensi Seasonality.

(A) Probability of mosquitoes dying upon exposure to each IRS compound in the time-period following spraying⁴² – pink indicates bendiocarb, turquoise indicates pirimiphos methyl. (B) For each temporal profile, the public health impact of annual IRS campaign with each insecticide according to the timing of the campaign. Points in grey correspond to the average reduction in incidence occurring from picking a random month to conduct the IRS campaign for each An. stephensi temporal profile, coloured points indicate the reduction in incidence arising from optimally timing the IRS campaign relative to the vector density peak for each An. stephensi temporal profile (and coloured arrows indicate the difference). The arrows and coloured points are coloured according to i) the insecticide used and ii) the degree of An. stephensi seasonality in each temporal profile (defined as the proportion of total annual abundance in any consecutive 4-month period). (C) For pirimithos methyl and for each An. stephensi temporal profile (coloured points) the percentage reduction in malaria incidence achieved if the IRS campaign is timed randomly, timed to start when rainfall is at its peak, or optimally timed based on peak vector density. Points correspond to specific An. stephensi temporal profiles and are coloured according to their degree of seasonality. Boxplot shows the minimum, first quartile, third quartile and maximum for the different individual projections. (D) As for (C), but for the insecticide bendiocarb.