Hierarchical models for West Nile virus risk mapping in Connecticut, USA

doi:10.21203/rs.3.rs-3293255/v1

Download PDF

Research Article

Hierarchical models for West Nile virus risk mapping in Connecticut, USA

https://doi.org/10.21203/rs.3.rs-3293255/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background: Predicting risk of West Nile virus (WNV) to humans in spaces without mosquito surveillance data is a key limitation of many WNV surveillance programs. To address this knowledge gap, we analyzed 20 years (2001 – 2020) of statewide, point-level mosquito and WNV surveillance data from Connecticut (CT), USA, using boosted regression trees (BRT) and generalized linear models (GLMs) to determine the most influential climate variables, land cover classes, and seasonality factors (such as Month of collection) associated with Culex pipiens abundance and WNV presence/absence in tested mosquito pools.

Methods: Candidate Cx. pipiens collection models were assessed based on explained deviance and root mean square error then optimized using a backward selection process. We then used predicted Cx. pipiens abundance in models of WNV presence/absence to predict WNV detection probabilities throughout CT. We validated these WNV predictions by testing the association between predicted WNV detection probabilities in mosquitoes and observed WNV incidence in mosquitoes from 2021 – 2022 and in humans from 2001 – 2022 using binomial-error generalized linear mixed effects models.

Results: Predicted mosquito WNV detection probabilities in unsampled spaces were significantly associated with the odds of a mosquito pool testing positive as well as a human case occurring within the geopolitical boundaries of a town.

Conclusion: This methodology has broad utility in the US and abroad to inform the public of risk of WNV quickly and easily in mosquitoes using only a few online and easily accessible data sources.

machine learning

West Nile virus

surveillance

Culex pipiens

risk mapping

Vector-borne disease (VBD) transmission is an inherently focal process that is influenced by a myriad of biological, climatic, environmental, and sociodemographic factors (1). Despite the ultra-fine scale nature of VBD transmission, patterns can emerge through the aggregation of data at increasing spatial and temporal levels (2, 3). These patterns may identify signals of public health risk which can be predicted (i.e., estimate future risk or risk unsampled space) across various time frames and geopolitical units. Such outcomes are central to public health messaging for VBDs. The diseases, datasets and methods used to understand VBD risks to humans vary widely, though all methods aggregate and extrapolate information from sampled to un(der)-sampled spaces. The West Nile virus (WNV) transmission system in the United States is an excellent case study for understanding and predicting spatiotemporal risk. Since the WNV was first detected in New York, New York in 1999, over 50,000 clinical cases have been reported in the United States (4, 5). West Nile virus is maintained by Culex spp. mosquitoes and birds in a zoonotic cycle, and humans are considered dead end hosts (i.e., development of viremia cannot be transferred to a blood feeding mosquito); human exposure to WNV occurs when wildlife transmission cycles spillover into human habitats and populations. While clinical infections of WNV are rare, it is estimated that only 1% of cases are reported (4, 6). Due to this large though uncertain burden of disease, there have been tremendous efforts to predict and project the spatiotemporal dynamics of WNV epidemics in the US and abroad.

Given the issues of scale at which mosquito and human data is collected and reported (respectively, very fine vs. very coarse spatiotemporal scales), it is difficult to predict risk of WNV in humans based on observed or predicted incidence of WNV in mosquitoes (7). This is driven by two primary difficulties: 1) human WNV cases are generally detected/reported after observing infections in mosquitoes, and 2) human cases can occur in the absence of collected/reported mosquito data. Researchers must aggregate mosquito data to the scale of human incidence (typically US counties for US studies) to determine the link between mosquito and human infections (8, 9). Predictive models of human WNV cases tend to analyze patterns of incidence at regional to country scales, often in the absence of mosquito data as an important predictor. At these scales, multiple research objectives have driven research into predicting WNV risk to humans, including determining climate factors that increase seasonal risk (10, 11), exploring range expansions under climate change scenarios (12), forecasting incidence at a county-level in real-time (13, 14), and developing methodologies that treat WNV clinical incidence as a probabilistic rather than deterministic process (15, 16). Predictive studies of mosquito WNV activity tend to be local, meaning they focus on examining the factors that influence mosquito infection rates at specific sites in an area typically no larger than a city or county (17–20). These studies have shown that hotter, drier, and more urban spaces (for Culex pipiens/quinquefasciatus) and suburban to rural spaces (for Culex tarsalis) strongly predict intensity of WNV mosquito infections (21–25). Studies of human risk also identify hot, dry summers as important risk factors for clinical incidence of WNV; how transmission intensity in zoonotic cycles links to incidence of human WNV cases at shorter temporal and smaller spatial extents remains an important knowledge gap.

To address how to predict WNV enzootic activity in unsampled spaces and how those predictions may relate to human incidence of WNV, we analyzed longitudinal mosquito, arbovirus, and human case data from the State of Connecticut (CT), USA using a series of nested boosted regression tree models. The Connecticut Agricultural Experiment Station’s (CAES) statewide mosquito and arbovirus surveillance system, established for eastern equine encephalitis virus (EEEV) surveillance in 1997, currently monitors mosquito populations and arbovirus activity at 108 sites in CT (26, 27): major expansions of the surveillance system occurred in 2001 and 2019. The CAES surveillance program uses standardized sampling methodologies across years and sites, and all captured female mosquitoes are identified to species and are tested for any of 9 arboviruses. Here, we used boosted regression tree models (28) informed by the point level surveillance data in the CAES program to generate monthly WNV detection probabilities across the entirety of Connecticut. Our goal was to develop a WNV risk projection pipeline to forecast the probability of WNV being detected in areas not covered by the surveillance system. To this end, we developed a series of nested hierarchical models to predict WNV detection probabilities in Culex pipiens mosquitoes. We then validated these models by examining the relationships between predicted town-level WNV detection probabilities and observed incidence of WNV in mosquitoes and humans.

All analyses took place in R V4.2.2 (29) using a variety of packages listed in each subsection. All spatial data used the WGS84 coordinate projection system. Conversions to this CRS, when needed, were performed using the corresponding functions available in the “raster”, “sp”, and “sf” packages. All displayed maps were created using a combination of the sf and ggplot2 packages (30, 31).

Mosquito and Arbovirus Surveillance Data: Mosquito and arbovirus data were obtained using hay-lactalbumin infusion-baited gravid traps at 87 surveillance sites in CT from 2001 – 2020. These traps are the most efficient surveillance tool for detecting WNV-infected Cx. pipiens mosquitoes (32), which is the primary enzootic and epidemic vector of WNV in the US northeast (33, 34). Briefly, sites were sampled by operating traps overnight on a 10-day rotation each summer from June through October, and all female mosquitoes were identified to species using a dichotomous key (35) and tested for nine arboviruses using viral cell culture and RT-PCR techniques (36). Collection sites were sampled more frequently (1-2 times weekly) for the remainder of a season if mosquitoes tested positive for an arbovirus of primary concern (WNV or EEEV). The average distance between analyzed sites was 61.3 km (+/- 0.42km).

Climate Data: Monthly temperature and precipitation records from 2001 to 2020 as well as monthly climate 30-year normals were obtained from PRISM (https://prism.oregonstate.edu/) (37) by extracting data to the specific surveillance site coordinates using the Data Explorer tool (https://prism.oregonstate.edu/explorer/). PRISM uses a default 4km grid system, meaning climate values were assigned to points based on the grid in which they are located. Since the CT surveillance program only operates traps from June to October, we extracted the current and 1-month lag values for each site that aligned with the mosquito collection period. We also calculated the absolute difference in the observed and normal values of these time intervals. We chose a 1-month lag in climate variables to allow for the serial interval of WNV infections in mosquitoes.

Landcover Data: Land cover data was obtained from UCONN CLEAR (https://clear.uconn.edu/projects/landscape/index.htm). Briefly, this site hosts data on twelve land cover classes for Long Island Sound from 1985 to 2015 (at intervals of 5 years). To build a suite of land cover variables for model training, shapefiles of the 2015 land cover estimate were processed in ArcMap (ESRI©) using a combination of clipping and spatial joins to produce a .csv of the area of each land cover type within ~500m of a site; percent land cover per trap buffer was then calculated in R. Published research has shown that overall land cover patterns have changed little in CT (38), so we did not include any other estimates of land cover concurrent with the surveillance period. Shapefiles for the political boundaries of CT were downloaded from CT.gov’s geospatial data hub (https://ct-deep-gis-open-data-website-ctdeep.hub.arcgis.com/search?groupIds=71c5c4a9c6ea4ea8ab54d1bf1faaeed8).

Population Density Data: Preliminary analyses with a training data set revealed one important aspect of the data: algorithms that did not directly include trapping effort (i.e., the total number of times a site was sampled in a month) as a predictor term performed poorly by explained deviance and root mean square error. Furthermore, algorithms that included trapping effort directly as a dependent variable often concluded trapping effort was the most important explanatory variable. This created a need to develop a prediction pipeline for trapping effort in unsampled spaces. To do this, we investigated the relationship between trapping effort and population density at any given surveillance site; we chose to analyze this relationship since CAES tends to choose surveillance sites near a variety of human-dominated habitats. We first downloaded the gridded global population density raster image for 2020 from Columbia University’s Socioeconomic Data and Applications Center (https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-adjusted-to-2015-unwpp-country-totals-rev11/data-download), cropped the image to match the spatial extent of CT, resampled the image to match the 4x4km resolution of the PRISM data, then extracted the average population density estimates within 5km of each surveillance site. We chose a 5km buffer since prior analyses with the CAES surveillance data set identified spatial patterns of synchronous WNV detection in mosquitoes up to 5km (39). We then used a Poisson-error generalized linear mixed effects model to confirm the strength of the association between trapping effort and population density while controlling for month of the year.

Hierarchical models: We used two types of modeling approaches. Our primary modeling approach was the use of gradient boosting machine, specifically boosted regression trees (BRT), in the “dismo” and “gbm” R packages (28).All processed surveillance, climate, land cover, and population density data were randomly assigned as training and test data using a 70/30 training/test ratio. Full models of Cx. pipiens trap night collections included all land cover and climate variables, Month of the year, and population density within 5 km of a surveillance site; full models of WNV detection (1 – one or more positive trap nights in a month, 0 – no positive trap nights in a month) included the same variables with the addition of the observed Cx. pipiens trap night collection. For the Cx. pipiens collections and WNV detection models, we first compared full models containing a subset of pairwise combinations of the climate variables and trapping location coordinates. In these subsets, only one climate form (average or difference from normal) per month or per climate variable (temperature or precipitation) was used; full models were also assessed with and without the x y coordinate of each surveillance site. We simplified these models through sequential elimination of the variables with lowest influence; the gbm package allows for model simplification up to X-2 variables. The best performing model for each variable was assessed from the vantage points of model deviance, root mean square error (RMSE, assessed using the test data), and simplicity (i.e., fewest variables). Once all BRT models were assessed, we compared their performance to an all pairwise interaction general linear model (GLM) using the same final terms as in the BRT models; we then simplified these GLMs using a backward selection procedure based on AIC (a similar function to model simplification for BRTs). The use of BRT and GLM models was done to compare, respectively, their ability (or lack of ability) to capture non-linear relationships. Final models utilized for validation were chosen based on lowest RMSE of the simplified BRT and GLMs.

Projecting Predictions: To generate monthly predicted WNV detection probabilities in Cx. pipiens mosquitoes, we first defined a 4 km x 4 km raster grid at the spatial extent of the state of CT. We downloaded the temperature and precipitation raster files from PRISM that matched the forms of the final models for the Years 2001 to 2022, re-projected the PRISM data to match the specified raster grid, and cropped the images to match the political boundaries of CT using a combination of functions available in the “sp”, “sf”, and “raster” packages (30, 40); we used the same methodology to generate a raster image of population density. To generate a percent land cover per raster cell for the entirety of CT, we downloaded a rasterized CT land cover image, converted all pixels to points, and listed the location and landcover value of each point in a data frame. Then for each land cover class, we counted the number of points within a 4km x 4km raster cell, and the percentage landcover within each cell was determined as the number of pixels per cell of each landcover type divided by the total number of pixels within the cell. For Cx. pipiens predictions, we created a RasterStack of all variables’ layers in the final models and generated predictions using the built in prediction function in the “raster” package for both BRT- and GLM-based models. The predicted vector abundance raster layer was then added to the RasterStack containing the variables in the final WNV model, and WNV detection probabilities were generated.

Validating predictions: We performed two different model validations. First, we compared the predicted probability of WNV detection in mosquitoes from our final candidate models to the observed presence/absence of WNV detections in mosquitoes for all sites in the CAES surveillance network in 2021 and 2022. To do this, we generated a data frame of WNV detection predictions for the summer months in towns with at least one trapping location. Because multiple raster cells could occur within a town’s political boundary, we assigned the maximum WNV detection probability to the town, rather than the average. Our rationale for this was that it communicated the maximum level of risk experienced by any resident within a town; using the maximum prediction may also better represent the actual occurrence of transmission within a town’s boundary. We then used binomial-error GLMMs to test the association between our predicted WNV detections and the observed WNV detections. Fixed effects included the current and prior month’s predicted WNV mosquito detection probability, year and month of collection, the town’s population size (log-transformed), and town name as a random effect (to control for repeated observations (41).

Our second validation approach was to compare predicted WNV detection probabilities in mosquitoes to observed detections of human cases from 2001 to 2022. A prime issue with human spillover of WNV is that cases can occur in spaces where little to no surveillance takes place. Because our primary objective was to determine if predicted WNV risk in unsampled spaces was a reliable indicator of WNV activity, testing the association between predicted mosquito detections and incidence of human cases was of great interest. We first generated a data frame of all town by month by year predictions of WNV detection in mosquitoes following a similar procedure for assigning raster values to geopolitical units as described for validating against mosquito data. We then used a GLMM to test the association between observed WNV human cases (modeled as detected 1, not detected 0) and our predicted WNV mosquito detection probabilities; we chose occurrence rather than absolute number of cases since there were rarely more than 1 human case per town for any year-month combination. Fixed effects included the current and prior month’s predicted WNV detection probability in mosquitoes, month of the year, the town’s population size (log-transformed), and town name as a random effect (to control for repeated observations).

To select the best model for both validation approaches, we first randomly split the predicted monthly WNV mosquito detection probability datasets into training and test data using a 70/30 ratio ensuring that data points were not resampled between the training and test data sets; for mosquito validation models, this included all data from 2021 – 2022 and for human validation models this included all data from 2001 - 2022. We then constructed a null model containing any seasonality terms (year and/or month), town population size, and the random effect terms. We then assessed improvements in AIC, BIC, and RMSE (assessed using the test data) due to the inclusion of single monthly detection terms (current and prior month’s predicted WNV detection). If both terms improved AIC, we then evaluated additive and interaction terms. If only one term improved AIC, we assessed whether interactions with the null model terms further improved model fit. GLMMs were built and assessed using a combination of the glmmTMB and sjPlot packages (42, 43).

From 2001 to 2020, CAES captured and tested 261,092 Cx. pipiens females in gravid traps, which resulted in 14,871 total mosquito pools of which 1,417 were WNV positive. During that time, 163 human cases were detected in CT: an additional 274 WNV positive Cx. pipiens pools and 13 human cases were detected between 2021 and 2022. Culex pipiens collections were generally higher in early summer, peaking in July, while peak incidence of WNV in mosquitoes occurred in August as did peak incidence of WNV in humans (Figure 1). West Nile virus detections in mosquitoes declined sharply from August to October while human incidence was similar between August and September (Figure 1). All three variables – Cx. pipens collections, WNV detection in mosquitoes, and WNV incidence in humans – were higher in the southwestern and central portions of CT (Figure 1).

Raster images of population density, percent landcover, and example temperature and climate variables are available as Supporting Information (SI) I: Hierarchical Model Data, as are other supporting images and tables referenced throughout the results section (SI II: Supplemental Results, Figures & SI III: Supplemental Results, Tables).

Relationship between population density and trapping effort: Because the CAES surveillance program tends to overall sample urban spaces (see (33)) and will increase trapping frequency at a site after detection of WNV or other arbovirus positive mosquitoes, we examined the relationship between population density and trapping frequency. Results from a Poisson-error GLMM confirmed a positive relationship between population density within 5 km of a surveillance site and the number of monthly trapping events (Figure 2, S. Table 1). Population density was thus used as a proxy for trapping effort in all subsequent models.

Culex pipiens trap night models: While site coordinate improved explained deviance in all Cx. pipiens trap night collection models, simplified models typically only kept the latitudinal coordinate as a predictor term (with trap night collections increasing along a North to South gradient); RMSE also tended to be higher for models when coordinate was included. The best performing BRT model by RMSE (0.29) included only two terms: population density (variable importance, VI, 51.0%) and the current month’s average temperature (VI 49.0%) (Figure 3, S. Figure 1). All GLMs of the simplified BRT models further reduced RMSE (Table 1); however, explained deviance in GLMs were lower compared to simplified BRTs. The best fitting GLM-based Cx. pipiens collection model was the same as the simplified BRT (Table 1), though the rank order of variables was different: average monthly temperature (VI 63%), population density (VI 22%), and the interaction between the two terms (VI 15%) (S. Figure 1). The observed predicted fit of the best fitting simplified BRT and GLM models was graphically similar (Figure 3 & S. Figure 2).

Table 1: Mosquito training data results. Boosted regression tree and general linear model (GLM) performance for the abundance of Culex pipiens females in gravid traps set at 87 sample sites in Connecticut, US from 2001 – 2020. Empty cells indicate the GLMs did not further simplify beyond the initial model. AVG – average observed monthly temperature (°C) or precipitation (mm); NRM – difference from observed 30-year normal, temperature (°C) or precipitation (mm); RMSE – root mean squared error

Climate Forms				Location Included	Full BRT Model (N variables = 20 when coordinate included)		Reduced Model
Climate Forms					Full BRT Model (N variables = 20 when coordinate included)		BRT			GLM comparison
Current Month		Prior Month			Estimated Deviance explained (%)	RMSE	Number variables	Estimated Deviance explained (%)	RMSE	Simplified BRT		Backwards Selection
Temperature	Precipitation	Temperature	Precipitation		Estimated Deviance explained (%)	RMSE	Number variables	Estimated Deviance explained (%)	RMSE	Estimated Deviance Explained (%)	RMSE	Estimated Deviance Explained (%)	RMSE
AVG	AVG	AVG	AVG	YES	53.5	0.88	11	54.0	1.06	51.9	0.40	49.26	0.13
AVG	AVG	AVG	AVG	NO	50.8	0.79	9	52.2	0.95	48.0	0.39	44.66	0.0053
NRM	NRM	AVG	AVG	YES	54.7	1.15	5	57.8	0.91	42.8	0.07	42.44	0.12
NRM	NRM	AVG	AVG	NO	52.3	1.07	4	56.7	0.76	40.8	0.12	---	---
NRM	AVG	NRM	AVG	YES	54.3	1.15	6	58.3	1.19	43.4	0.17	43.03	0.10
NRM	AVG	NRM	AVG	NO	51.6	1.24	5	54.6	0.89	40.8	0.13	---	---
AVG	NRM	AVG	NRM	YES	52.8	1.04	11	54.0	1.33	50.1	0.42	49.73	0.36
AVG	NRM	AVG	NRM	NO	49.4	0.96	2	50.8	0.29	42.3	0.001	---	---
AVG	AVG	NRM	NRM	YES	52.3	0.84	10	54.2	1.43	48.2	0.19	47.85	0.073
AVG	AVG	NRM	NRM	NO	47.8	1.29	3	49.3	0.49	42.3	0.01	---	---
NRM	NRM	NRM	NRM	YES	55.0	1.20	6	57.7	1.30	56.5	0.09	43.08	0.11
NRM	NRM	NRM	NRM	NO	52.2	1.33	6	54.6	1.28	41.7	0.14	41.55	0.061

West Nile virus detection models: There was no clear influence of trap coordinate on model performance, and explained deviance and RMSE were generally similar for BRT models of WNV detections in mosquitoes with and without coordinate. The best performing BRT-based WNV model by RMSE included five terms: Cx. pipiens trap night collections (VI 44%), month of the year (VI 19.5%), population density (VI 17.5%), average current monthly precipitation (VI 10%), and current difference from monthly normal temperatures (VI 9%) (Figure 4). This model also ranked two influential two-way interaction terms: one with month of the year and population density and one with Cx. pipiens trap night collections and the monthly temperature difference from normal (Figure 4). GLMs of the simplified BRT models did reduce RMSE (Table 2) yet explained deviance in GLMs was lower than the simplified BRTs. Furthermore, there was a noticeable difference in explained deviance and RMSE between the simplified BRT and the corresponding GLMs: explained deviance tended to increase and RMSE tended to decrease with increasing number of final variables in the GLM versions (Table 2). The best fitting simplified GLM by RMSE included eight terms: Cx. pipiens trap night collections, month of the year, population density, temperature difference from normal for the current and prior month, average current monthly precipitation, latitudinal trap coordinate, and the percent land classified as developed surrounding a trap site (Figure 5). Where terms were equal or of a similar nature (such as average verses different from normal), prediction trends were mostly similar between the best performing BRT and GLM (Figures 4 & 5).

Table 2: West Nile virus training data results. Boosted regression tree and general linear model (GLM) performance for the probability of West Nile virus detection in Culex pipiens females in gravid traps set at 87 sample sites in Connecticut, US from 2001 – 2020. When a full model containing trap coordinates simplified to a model without coordinates, the full model versions without coordinates were not analyzed. This is represented by the dashed lines. Cells containing a backslash indicate that the reduced model simplified to a version already discovered in a different variable combination run. AVG – average observed monthly temperature (°C) or precipitation (mm); NRM – difference from observed 30-year normal, temperature (°C) or precipitation (mm); RMSE – root mean squared error

Climate Forms				Location Included	Full BRT Model (N variables = 21 when coordinate included)		Reduced Model
Climate Forms					Full BRT Model (N variables = 21 when coordinate included)		Simplified BRT			GLM comparison
Current Month		Prior Month			Estimated Deviance explained (%)	RMSE	Number variables	Estimated Deviance explained (%)	RMSE	Simplified BRT		Backwards selection
Temperature	Precipitation	Temperature	Precipitation		Estimated Deviance explained (%)	RMSE	Number variables	Estimated Deviance explained (%)	RMSE	Estimated Deviance Explained (%)	RMSE	Estimated Deviance Explained (%)	RMSE
AVG	AVG	AVG	AVG	YES	44.58	0.0030	5	45.31	0.0017	41.67	0.00098	41.67	0.00098
AVG	AVG	AVG	AVG	NO	---	---	---	---	---	---	---	---	---
NRM	NRM	AVG	AVG	YES	45.28	0.0019	5	45.75	0.00081	40.19	0.00088	39.62	0.0010
NRM	NRM	AVG	AVG	NO	---	---	---	---	---	---	---	---	---
NRM	AVG	NRM	AVG	YES	44.82	0.0023	8	46.15	0.0012	43.02	0.00013	42.21	0.00022
NRM	AVG	NRM	AVG	NO	44.92	0.0019	5	45.47	0.00073	40.21	0.00095	39.84	0.00092
AVG	NRM	AVG	NRM	YES	45.09	0.0027	5	45.47	0.0014	41.62	0.00089	41.37	0.00091
AVG	NRM	AVG	NRM	NO	---	---	---	---	---	---	---	---	---
AVG	AVG	NRM	NRM	YES	45.05	0.0032	5	/	/	/	/	/	/
AVG	AVG	NRM	NRM	NO	---	---	---	---	---	---	---	---	---
NRM	NRM	NRM	NRM	YES	45.23	0.002	17	45.12	0.0017	49.76	0.00096	Did not converge
NRM	NRM	NRM	NRM	NO	45.13	0.0016	5	/	/	/	/	/	/

Projecting WNV detections into unsampled space: Because of the similar predictive profiles yet minor differences in explained deviance and RMSE among simplified BRTs and GLMs, we projected WNV detections using all possible pair-wise combinations of the best fitting BRT and GLM training models by RMSE for Cx. pipiens collections and WNV detections.

Visual assessments of the WNV detection probabilities confirmed known geographies and ecologies of the pathogen in CT. The highest probabilities of detecting WNV in Cx. pipiens mosquitoes concentrated in coastal regions in the southwest portion of CT (including towns in Fairfield and New Haven counties) as well as in central CT (predominately the Hartford metropolitan area) (Figure 6 & 7). The form of the WNV detection model had the greatest influence on visual patterns of predicted WNV presence/absence in mosquitoes. Overall, GLM-based WNV models generated a greater range of detection probabilities yet restricted maximum detections to a few core areas (Greenwich/Stamford, Bridgeport, New Haven, and Hartford) (Figure 6 & 7). BRT-based WNV models generated an overall smaller range of detection rates yet also identified a greater prevalence of detection rates above the average in more towns.

West Nile virus detection model validation: Validation approaches confirmed the utility of predicted WNV detection probabilities in mosquitoes for forecasting risk of WNV activity in into unsampled spaces (Figures 8 & 9). For WNV presence/absence in mosquito validations, all but two candidate models were either equivalent or ΔAIC < 2 to the null model containing population size, month and year, and town as a random effect (S. Table 2). The two models that did improve ΔAIC > 2 contained a single term for the predicted detection of WNV in mosquitoes in the current month: both candidate models generated WNV predictions from a GLM-based model while Cx. pipiens trap night collections were generated from either a GLM- or BRT-based model. Predictive fits from the model displaying the lowest AIC are shown in Figure 8.

For presence/absence of a human case within a town, six candidate models improved ΔAIC > 2 compared to the null model containing town population size and month as a seasonal term (S Table 3). Four of these candidate models contained only a single term – the hierarchical model’s predicted prior month’s WNV presence in mosquitoes. Predictive fits from the single term model displaying the lowest AIC are shown in Figure 9. The best fitting model by AIC and RMSE contained the nested BRT-CXP and BRT-WNV training model’s prior month’s WNV detection rate in mosquitoes and an interaction between this term and town population size (S. Figure 3). Null incidence among towns, based on an assessment of random effect intercepts, were equivalent between the two candidate models of human incidence (S. Figure 4).

Inability to estimate risk of VBDs in un(der)-sampled spaces is a critical limitation of many VBD surveillance systems. To enhance the utility of a current public health surveillance system to project risk of disease within the entire political boundary of a US state, we developed methodologies that utilized nested models to first predict trap collections of Cx. pipiens mosquitoes throughout CT; these predictions were then used to predict the probability of WNV being detected in mosquitoes. Our results were successfully validated against patterns of WNV detection in mosquitoes collected from data not included as training data. Our results also showed a small, though significant, association between our predicted presence of WNV in mosquitoes and the occurrence of human WNV cases at the level of a town. Thus, the surveillance data pipeline developed in this report provides a validated tool to quantify and communicate the risks of WNV to humans based on information obtained from mosquito samples that is then interpolated into unsampled spaces.

To determine if WNV is circulating at any given point in space and time, surveillance systems must maximize their collections of the primary mosquito vector, in our case Cx. pipiens. Thus, our risk projection pipeline was sensitive to our Cx. pipiens modeling methodologies. To approximate trapping effort, we examined the relationship between the number of monthly trapping events and the density of humans within 5km of a surveillance site. We chose this radius since prior work with the CAES dataset identified 5km as the upper limit of spatiotemporal correlations in WNV detections (39). We confirmed a strong positive relationship between trapping effort and population density, and all resulting best performing models included population density as an important predictor of Cx. pipiens collections and incidence of WNV in mosquitoes and humans. This important association between population density, Cx. pipiens, and WNV supports the overwhelming evidence that, in the eastern half of the US, urban environments are the most at risk for WNV. This result also compliments prior ecological investigations of Cx. pipiens and WNV in CT, which identified non-forested habitats as the best predictors of Cx. pipiens collections and WNV detections; the highest risk sites are also in the State’s primary urban centers (33, 44). Despite these known associations, few to no specific land cover classes were not included in any final candidate training model of Cx. pipiens collections or WNV detections. Simplified models that did contain land cover variables found increased collections in sites with a greater prevalence of developed and open water (i.e., urban coastal) or decreased collections in sites with a greater prevalence of forest land cover classes. For the Cx. pipiens training models, it could be that land cover is highly correlated with population density and temperature, such that these two variables were sufficient to explain the observed variation in collections. Our methods may have also been more sensitive to other measures of urban sites which provided greater predictive precision. Given that climate data, in general, is more freely accessible online and captured at greater frequencies than land cover estimates, more detailed analyses of climatic data sources for mosquito species distributions may prove a more powerful approach to projecting vector-borne disease risk

The predicted functional relationships in the final Cx. pipiens and WNV models support prior studies showing elevated Cx. pipiens collections and WNV detections during the hotter summer months (21, 22, 45). More specifically, the models predicted Cx. pipiens collections were highest when average monthly temperatures approached or exceeded ~ 25°C. A similar trend was observed for WNV detections in mosquitoes, however from the perspective of temperatures different from normal. These results support by prior work on the thermal biology of vector-borne diseases, which has identified temperatures between 23 and 25°C as optimal temperatures for WNV activity (17, 46, 47). The additional variables included in the WNV training models, and their corresponding rank order, likely has to do with the form of risk modeled (presence/absence of detection) and the objectives of our research. Seasonality and population density were the second and third ranking order variables in the final BRT-based WNV model, and each variable carried more influence than the current month’s average precipitation and the temperature difference from normal. This could be because we modeled the presence/absence of WNV rather than the magnitude of WNV infections. In any given year, a greater number of WNV positive pools are detected in a greater number of sites in August, which would reflect the rank importance of month in the BRT-based model. Modeling presence/absence of WNV in mosquitoes could further explain why elevated precipitation levels were associated with increased WNV detections: if WNV detections are sensitive to the number of captured and collected Cx. pipiens and Cx. pipiens populations in general increase with increasing rainfall, then there is a greater chance of detecting WNV in non-traditional sites.

Rank importance of the variables in the GLM-based WNV training model were far more complicated than the BRT-based models. Because the GLMs allowed for all pair-wise combinations of variables, variable importance scores were much lower across all variables and many interaction terms were ranked more influential than single, additive terms. Prediction ranges for Cx. pipiens collections and WNV detections from GLM-based models were also generally larger compared to the BRT-based models, possibly because GLMs allow for only linear relationships which may have over or under-predicted collections and detections at the limits of a variables range. Despite these differences, GLM-based WNV training models were more often associated with improved fit when validated against 2021–2022 mosquito data and 2001–2022 human data. This again could be because GLM-based training models allowed for a greater upper limit of WNV detection probabilities, which then provided a slightly better fit to the validation data sets.

It is important to note that many of the best performing training models for Cx. pipiens collection and WNV detections did not include any climate lags as important predictor terms. Other studies have identified lags in extreme events or used measures of climatic variability as important predictors of Cx. pipiens collections (12, 48, 49). Many of these studies’ objectives focused on seasonal variation in WNV transmission intensity and relied on more sophisticated methodologies, which applied in this investigation may have achieved a better fit to the data. While extreme events are difficult to forecast, they can have greater proportional impacts in urban spaces, especially when it comes to extreme heat. Such results would also not change the primary message of our study, which is that the probability of detecting WNV-infected Cx. pipiens mosquitoes is highest in urbanized spaces during the hottest months of the year. Another component missing from our risk estimation pipeline concerns secondary vectors, and prior work with the CAES surveillance data set has identified Culex restuans, Culex salinarius, and Culiseta melanura as important secondary vectors of WNV in CT (33). We suspect that methods that project risk into Cx. restuans and Cx. salinarius would mostly identify similar risk factors as for Cx. pipiens. Analyses of Cs. melanura could identify expansions of risk later in the summer into non-urbanized spaces. This is based on prior knowledge of Cs. melanura in CT (50, 51), though future studies of WNV incidence in all three species, as well studies into the relationships between vector community composition and WNV transmission, are needed to determine the exact associations between collections of multiple species and WNV detections.

We validated our WNV risk estimation pipeline against observed monthly incidence WNV in mosquitoes and humans at the level of a town, which, for human-based studies, is at a much smaller spatial scale and finer temporal resolution than many prior studies (15). For mosquitoes, only two candidate models containing the current month’s predicted probability of detecting WNV mosquitoes improved fit beyond a null model that accounted for town population size and month of the year. The lack of improved fits from all candidate prediction models could be due to analyzing the validation data at the scale of a town; some towns contain multiple trapping locations, and validating at the scale of a surveillance site may have improved fit. This could also be due to the overall dearth of validation data compared to the training data. We subset our initial training data set into training and test data, and all RMSE values were calculated against the test data set. Thus, we do evidence that our training models fit the test data at the scale of a surveillance site. For human clinical data, overall incidence of WNV clinical cases is extremely low, such that our predicted WNV detections in mosquitoes were only minimally associated with presence/absence of a human case. Despite the small increase in overall odds of detecting a human case based on WNV infections in mosquitoes, these associations remained significant after accounting for town population size and month of the year. Other metrics of human populations, such as demography, social inequity, and access to healthcare are other important metrics to consider when predicting incidence of WNV in humans (52–54); these variables can also influence patterns of mosquito community composition and collections (55, 56). Future studies that examine these relationships may refine the predicted incidence of WNV in mosquitoes and humans.

Together, the primary risk factors for detecting a human case of WNV at any given town in CT include month of the year (highest incidence is in August and September), population size (more opportunities for a case to occur as well as an indication of urbanization), and a high probability of WNV activity in mosquitoes in the prior month. Mosquito-based surveillance remains the gold standard for assessing the risk of WNV activity in mosquitoes; surveillance provides early evidence of virus activity as well as information on the abundance, distribution, identity, and infection rates of potential mosquito vectors. Nevertheless, mosquito surveillance is labor intensive and is limited by practical constraints on the number of locations that can be sampled. Our presented methods provide a means to extend the utility of information obtained from point-source surveillance information by predicting WNV risk to the public in regions that are not monitored by traps. Our approach uses only a handful of freely available data sets on climate and socioeconomics, and these data sets and methods could be used to identify high risk areas of WNV in the northeast that are not covered by an existing surveillance network.

BRT: boosted regression tree model

CT: State of Connecticut

GLM: generalize linear model

GLMM: generalized linear mixed effects model

WNV: West Nile virus

Human data: Any visualized and presented human data is freely available from online sources. There are additionally no identifiers provided.

Author contributions: JRM, LFC and PMA conceived of the research; JRM and LFC performed the analyses; JRM wrote the manuscript first draft; all authors contributed to manuscript writing and approved the submitted version; JRM and PMA acquired funding for the project.

Acknowledgements: This work would not be possible without the contributions of numerous field and laboratory seasonals employed by CAES during each surveillance season. Special mention is given to John Shepard, Tanya Petroff, Angela Bransfield, Michael Misencik, and Theodore Andreadis for their commitment to the past and present operation of the CAES surveillance network.

Data availability: All data is archived in Mendeley Data and available for review upon request. Data will be released with a citable DOI at time of acceptance.

Conflicts of interest: None

Funding: This publication was supported in part by the American Mosquito Control Association 2022 Research Fund Number 2022-02 and Cooperative Agreement Number U01CK000509, funded by the Centers for Disease Control and Prevention. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of the American Mosquito Control Association, the Centers for Disease Control and Prevention or the Department of Health and Human Services. LFC conducted this work while a Scholar in the NIH Climate and Health Scholars Program with additional support from Indiana University. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Reisen WK. Landscape epidemiology of vector-borne diseases. Annu Rev Entomol. 2010;55:461-83.
Chaves LF, Harrington LC, Keogh CL, Nguyen AM, Kitron UD. Blood feeding patterns of mosquitoes: random or structured? Front Zool. 2010;7(1):3.
Kitron U, Clennon JA, Cecere MC, Gurtler RE, King CH, Vazquez-Prokopec G. Upscale or downscale: applications of fine scale remotely sensed data to Chagas disease in Argentina and schistosomiasis in Kenya. Geospat Health. 2006;1(1):49-58.
Soto RA, Hughes ML, Staples JE, Lindsey NP. West Nile Virus and Other Domestic Nationally Notifiable Arboviral Diseases - United States, 2020. MMWR Morb Mortal Wkly Rep. 2022;71(18):628-32.
McDonald E, Mathis S, Martin SW, Staples JE, Fischer M, Lindsey NP. Surveillance for West Nile Virus Disease - United States, 2009-2018. MMWR Surveill Summ. 2021;70(1):1-15.
Petersen LR, Carson PJ, Biggerstaff BJ, Custer B, Borchardt SM, Busch MP. Estimated cumulative incidence of West Nile virus infection in US adults, 1999-2010. Epidemiol Infect. 2013;141(3):591-5.
Uelmen JA, Irwin P, Bartlett D, Brown W, Karki S, Ruiz MO, et al. Effects of Scale on Modeling West Nile Virus Disease Risk. Am J Trop. Med. Hyg. 2021;104(1):151-65.
Keyel AC, Raghavendra A, Ciota AT, Elison Timm O. West Nile virus is predicted to be more geographically widespread in New York State and Connecticut under future climate change. Glob Change Biol.n/a(n/a).
Keyel AC, Elison Timm O, Backenson PB, Prussing C, Quinones S, McDonough KA, et al. Seasonal temperatures and hydrological conditions improve the prediction of West Nile virus infection rates in Culex mosquitoes and human case counts in New York and Connecticut. PloS one. 2019;14(6):e0217854.
DeGroote JP, Sugumaran R, Ecker M. Landscape, demographic and climatic associations with human West Nile virus occurrence regionally in 2012 in the United States of America. Geospat Health. 2014;9(1):153-68.
Wimberly MC, Lamsal A, Giacomo P, Chuang TW. Regional variation of climatic influences on West Nile virus outbreaks in the United States. Am J Trop Med Hyg. 2014;91(4):677-84.
Farooq Z, Rocklöv J, Wallin J, Abiri N, Sewe MO, Sjödin H, et al. Artificial intelligence to predict West Nile virus outbreaks with eco-climatic drivers. Lancet. 2022;17:100370.
DeFelice NB, Little E, Campbell SR, Shaman J. Ensemble forecast of human West Nile virus cases and mosquito infection rates. Nat Commun. 2017;8(1):14592.
DeFelice NB, Birger R, DeFelice N, Gagner A, Campbell SR, Romano C, et al. Modeling and Surveillance of Reporting Delays of Mosquitoes and Humans Infected With West Nile Virus and Associations With Accuracy of West Nile Virus Forecasts. JAMA Netw Open. 2019;2(4):e193175.
Keyel AC, Gorris ME, Rochlin I, Uelmen JA, Chaves LF, Hamer GL, et al. A proposed framework for the development and qualitative evaluation of West Nile virus models and their application to local public health decision-making. PLoS Neg Trop Dis. 2021;15(9):e0009653.
Holcomb KM, Mathis S, Staples JE, Fischer M, Barker CM, Beard CB, et al. Evaluation of an open forecasting challenge to assess skill of West Nile virus neuroinvasive disease prediction. Parasit Vector. 2023;16(1):1-13.
Skaff NK, Cheng Q, Clemesha RES, Collender PA, Gershunov A, Head JR, et al. Thermal thresholds heighten sensitivity of West Nile virus transmission to changing temperatures in coastal California. Proc Biol Sci. 2020;287(1932):20201065.
Little E, Campbell SR, Shaman J. Development and validation of a climate-based ensemble prediction model for West Nile Virus infection rates in Culex mosquitoes, Suffolk County, New York. Parasit Vector. 2016;9(1):1-13.
Shand L, Brown WM, Chaves LF, Goldberg TL, Hamer GL, Haramis L, et al. Predicting West Nile Virus Infection Risk From the Synergistic Effects of Rainfall and Temperature. J Med Entomol. 2016;53(4):935-44.
Poh KC, Chaves LF, Reyna-Nava M, Roberts CM, Fredregill C, Bueno Jr R, et al. The influence of weather and weather variability on mosquito abundance and infection with West Nile virus in Harris County, Texas, USA. SciTotal Environ. 2019;675:260-72.
Ruiz MO, Chaves LF, Hamer GL, Sun T, Brown WM, Walker ED, et al. Local impact of temperature and precipitation on West Nile virus infection in Culex species mosquitoes in northeast Illinois, USA. Parasit Vectors. 2010;3(1):19.
Chaves LF, Hamer GL, Walker ED, Brown WM, Ruiz MO, Kitron UD. Climatic variability and landscape heterogeneity impact urban mosquito diversity and vector abundance and infection. Ecosphere. 2011;2(6).
Reisen WK, Carroll BD, Takahashi R, Fang Y, Garcia S, Martinez VM, et al. Repeated West Nile virus epidemic transmission in Kern County, California, 2004-2007. J Med Entomol. 2009;46(1):139-57.
Rochlin I, Faraji A, Healy K, Andreadis TG. West Nile Virus Mosquito Vectors in North America. J Med Entomol. 2019;56(6):1475-90.
Holcomb KM, Reiner RC, Barker CM. Spatio-temporal impacts of aerial adulticide applications on populations of West Nile virus vector mosquitoes. Parasit Vector. 2021;14(1):120.
Brown SC, Cormier J, Tuan J, Lier AJ, McGuone D, Armstrong PM, et al. Four Human Cases of Eastern Equine Encephalitis in Connecticut, USA, during a Larger Regional Outbreak, 2019. Emerg Infect Dis. 2021;27(8).
Morens DM, Folkers GK, Fauci AS. Eastern Equine Encephalitis Virus - Another Emergent Arbovirus in the United States. N Engl J Med. 2019;381(21):1989-92.
Elith J, Leathwick JR, Hastie T. A working guide to boosted regression trees. J Anim Ecol. 2008;77(4):802-13.
R Development Core Team R. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2008.
Pebesma EJ. Simple features for R: standardized support for spatial vector data. R J. 2018;10(1):439.
Wickham H, Chang W, Wickham MH. Package ‘ggplot2’. Create elegant data visualisations using the grammar of graphics Version. 2016;2(1):1-189.
Reiter P. A portable battery-powered trap for collecting gravid Culex mosquitoes. Mosq News. 1983;43(4):496-8.
McMillan JR, Armstrong PM, Andreadis TG. Patterns of mosquito and arbovirus community composition and ecological indexes of arboviral risk in the northeast United States. PLoS Neg Trop Dis. 2020;14(2):e0008066.
Andreadis TG. The contribution of Culex pipiens complex mosquitoes to transmission and persistence of West Nile virus in North America. J Am Mosq Control Assoc. 2012;28(4 Suppl):137-51.
Darsie R, Ward R. Identification and geographical distribution of the mosquitoes of North America, north of Mexico: University Press of Florida; 1981.
Armstrong PM, Andreadis TG, Finan SL, Shepard JJ, Thomas MC. Detection of infectious virus from field-collected mosquitoes by vero cell culture assay. J Vis Exp. 2011(52).
Daly C, Halbleib M, Smith JI, Gibson WP, Doggett MK, Taylor GH, et al. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Ont. J. Climatol. 2008;28(15):2031-64.
Arnold C, Wilson E, Hurd J, Civco D. 30 years of land cover change in Connecticut, USA: A case study of long-term research, dissemination of results, and their use in land use planning and natural resource conservation. Land. 2020;9(8):255.
McMillan JR, Chaves LF, Armstrong PM. Ecological predictors of mosquito population and arbovirus transmission synchrony estimates. J Med Entomol. 2023;60(3):564-74.
Hijmans RJ, Van Etten J, Mattiuzzi M, Sumner M, Greenberg J, Lamigueiro O, et al. Raster package in R. Version https://mirro rs sjtug sjtu edu cn/cran/web/packa ges/rast e r/raster pdf. 2013.
Chaves LF. An entomologist guide to demystify pseudoreplication: data analysis of field studies with design constraints. J Med Entomol. 2010;47(3):291-8.
Brooks ME, Kristensen K, Van Benthem KJ, Magnusson A, Berg CW, Nielsen A, et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 2017;9(2):378-400.
Ludecke D. sjPlot. Data visualization for statistics in social science. R package version 241. 2018.
Diuk-Wasser MA, Brown HE, Andreadis TG, Fish D. Modeling the spatial distribution of mosquito vectors for West Nile virus in Connecticut, USA. Vector Borne Zoonotic Dis. 2006;6(3):283-95.
Bisanzio D, Giacobini M, Bertolotti L, Mosca A, Balbo L, Kitron U, et al. Spatio-temporal patterns of distribution of West Nile virus vectors in eastern Piedmont Region, Italy. Parasit Vector. 2011;4.
Mordecai EA, Caldwell JM, Grossman MK, Lippi CA, Johnson LR, Neira M, et al. Thermal biology of mosquito-borne disease. Ecol Lett. 2019;22(10):1690-708.
McMillan JR, Blakney RA, Mead DG, Koval WT, Coker SM, Waller LA, et al. Linking the vectorial capacity of multiple vectors to observed patterns of West Nile virus transmission. J Appl Ecol. 2019;56(4):956-65.
Lim AY, Cheong HK, Chung Y, Sim K, Kim JH. Mosquito abundance in relation to extremely high temperatures in urban and rural areas of Incheon Metropolitan City, South Korea from 2015 to 2020: an observational study. Parasit Vectors. 2021;14(1):559.
Johnson BJ, Sukhdeo MV. Drought-induced amplification of local and regional West Nile virus infection rates in New Jersey. J Med Entomol. 2013;50(1):195-204.
Armstrong PM, Andreadis TG. Ecology and Epidemiology of Eastern Equine Encephalitis Virus in the Northeastern United States: An Historical Perspective. J Med Entomol. 2022;59(1):1-13.
Skaff NK, Armstrong PM, Andreadis TG, Cheruvelil KS. Wetland characteristics linked to broad-scale patterns in Culiseta melanura abundance and eastern equine encephalitis virus infection. Parasit Vectors. 2017;10(1):501.
Chaves LF, Cohen JM, Pascual M, Wilson ML. Social Exclusion Modifies Climate and Deforestation Impacts on a Vector-Borne Disease. PLoS Negl Trop Dis. 2008;2(2):e176.
Bisanzio D, Martello E, Izenour K, Stevens K, Kaur R, McKenzie BA, et al. Arboviral diseases and poverty in Alabama, 2007-2017. PLoS Negl Trop Dis. 2021;15(7):e0009535.
Ruiz MO, Tedesco C, McTighe TJ, Austin C, Kitron U. Environmental and social determinants of human risk during a West Nile virus outbreak in the greater Chicago area, 2002. Int J Health Geogr. 2004;3(1):8.
Katz G, Leisnham PT, LaDeau SL. Aedes albopictus Body Size Differs Across Neighborhoods With Varying Infrastructural Abandonment. J Med Entomol. 2020;57(2):615-9.
Rothman SE, Jones JA, LaDeau SL, Leisnham PT. Higher West Nile Virus Infection in Aedes albopictus (Diptera: Culicidae) and Culex (Diptera: Culicidae) Mosquitoes From Lower Income Neighborhoods in Urban Baltimore, MD. J Med Entomol. 2021;58(3):1424-8.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Hierarchical models for West Nile virus risk mapping in Connecticut, USA

Status:

Version 1

Abstract

Figures

Background

Methods

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1