The importance of human population characteristics in modeling Aedes aegypti distributions and assessing risk of mosquito-borne infectious diseases
© The Author(s) 2017
Received: 19 July 2017
Accepted: 30 October 2017
Published: 15 November 2017
The mosquito Aedes aegypti has long been a vector for human illness in the Southeastern United States. In the past, it has been responsible for outbreaks of dengue, chikungunya, and yellow fever and, very recently, the Zika virus that has been introduced to the region. Multiple studies have modeled the geographic distribution of Ae. aegypti as a function of climate factors; however, this ignores the importance of humans to the anthropophilic biter. Furthermore, Ae. aegypti thrives in areas where humans have created standing water sites, such as water storage containers and trash. As models are developed to examine the potential impact of climate change, it becomes increasingly important to include the most comprehensive set of predictors possible.
This study uses Maxent, a species distribution model, to evaluate the effects of adding poverty and population density to climate-only models. Performance was evaluated through model fit statistics, such as AUC, omission, and commission, as well as individual variable contributions and response curves. Models which included both population density and poverty exhibited better predictive power and produced more precise distribution maps. Furthermore, the two human population characteristics accounted for much of the model contribution—more so than climate variables.
Modeling mosquito distributions without accounting for their dependence on local human populations may miss factors that are very important to niche realization and subsequent risk of infection for humans. Further research is needed to determine if additional human characteristics should be evaluated for model inclusion.
As the 1963 training manual on “Mosquitoes of Public Health Importance and Their Control” puts it, “Mosquitoes have probably had a greater influence on human health and well-being throughout the world than any other insects” . This is particularly true in the USA for the mosquito Aedes aegypti, the vector for Zika virus, dengue, chikungunya, and yellow fever. It has long been, and continues to be, the focus of much of the local concern for human health [2–4]. Furthermore, the Southeastern region has historically borne the brunt of these infectious diseases when they appear within the conterminous USA. The only outbreak of chikungunya in the USA occurred in Florida in 2014 and, before that, dengue was spread locally in Texas in 2005 [5, 6]. Historically, yellow fever caused several major outbreaks in the Southeast, particularly along the Mississippi River , with the last major outbreak occurring in New Orleans in 1905 . And, now the USA is bracing itself for the spread of Zika virus from Central and South America into the Southeast .
Mosquitoes have specific habitat requirements that vary based on each particular species. For Ae. aegypti, the preferred habitat is in tropical and subtropical climates, but a recent study estimated that the current distribution does extend into some temperate regions . While Ae. aegypti is sensitive to temperature, other climatic factors such as precipitation and altitude-influenced proxies are important as well [4, 11]. However, a close examination of the life cycle of this mosquito species also reveals that anthropogenic factors may strongly influence their geographic distribution. Ae. aegypti is a species that rests indoors and takes its blood meals primarily from humans . It also lays its eggs in containers of standing water, which are often found near the humans who use those containers to store their water . Poverty may be an indicator of areas with increased amounts of standing water, due to a lack of local sanitary services. Poverty may also portend a lack of air conditioning, which would raise the need for opening windows, thereby increasing mosquito entry. As these examples show, many of the predictors of mosquito-human interactions are affected by socioeconomic status within an area [4, 14, 15].
Previous modeling efforts focused primarily on temperature and other climatic variables for niche identification in Ae. aegypti distribution models. Campbell et al.  modeled Ae. aegypti as a vector for chikungunya and dengue globally and their models restricted parameters entirely to climatic variables, and Brady et al.  only examined temperature as a predictor of vector competence for dengue. In modeling the distribution of Ae. aegypti globally using boosted tree regression models, Kraemer et al.  explored the contributions of enhanced vegetation indices, climatic variables, and urbanization. To model urbanization, they generated a categorical variable from the Global Rural Urban Mapping Project (GRUMP), population density metrics, and night-time satellite imagery displaying land usage. The categories created consisted of urban, peri-urban, and rural classifications. They found a relatively small effect of urban designation on habitat suitability but hypothesized that satellite data showing human habitats may be more informative . However, each of these modeling attempts identified roughly the same geographic regions within the Southeastern United States as suitable for Ae. aegypti, but while it is informative to know where the mosquito can live, this does not help identify where individuals are at the highest risk of interactions with mosquitoes.
To understand the full picture of niche suitability and appropriately model the risk of mosquito-human interaction, it is necessary to model possible niches as a function of both bioclimatic and human population characteristics. As Eisen and Moore  indicated, understanding how temperature and precipitation affect suitability does not account for factors that may confound the relationship between niche suitability and niche realization. Restricting models to climate variables, while convenient, oversimplifies the actual niche dynamics of the species . This problem is magnified when models are projected for climate change. Human population is changing, both due to population expansion and the impacts of climate change and conflict-related migration [18, 19]. It is unreasonable to assume this will not affect the interactions between humans and mosquitoes.
To address this divide, this analysis seeks to examine how human population characteristics impact species distribution modeling for Ae. aegypti. This analysis will use bioclimatic variables identified in previous modeling efforts  for an initial model. Then, additional variables that account for population density, as a measure of availability of feeding opportunities, and poverty, as a surrogate for likelihood of human-mosquito interaction, will be added and tested for model performance and potential improvement.
While Ae. aegypti is present in a variety of global locations, modeling a small and well-characterized region can increase the predictive capabilities of a model , which should be more appropriate here, in testing variable contribution to model fit. Since the Southeastern United States is highly vulnerable to Ae. aegypti and has well-catalogued climate, population density, and poverty data, this region will be used to test the relative effects of these variables on the model. The results of this model will inform future Ae. aegypti models with respect to the importance of human population characteristics.
Occurrence data for Ae. aegypti was taken from the dataset available through the “Global Compendium of Aedes aegypti and Aedes albopictus Occurrence” . This compendium was based primarily on a review of published occurrence datasets, but also included surveillance system data, national mosquito surveys, and records from governmental health agencies, where available. All occurrences were recorded between 1960 and 2014 and contained latitude and longitude for georeferencing . Since the dataset contains occurrences with georeferencing and metadata information but no information on sampling or absence, this is a presence-only dataset. To avoid issues common to “opportunistic” sampling (i.e., oversampling of more accessible and higher risk areas), Kraemer et al.  applied methods of spatial and temporal standardization. For example, a single occurrence was defined as an observation within one calendar year at a given unique location (5 km2). This resulted in 1112 occurrences of Ae. aegypti being removed from the global compendium dataset .
Publicly available data from WorldClim.org were used for current bioclimatic variables. This dataset contains 19 variables known as the “bioclim” set, which includes characteristics designed to model annual trends, such as annual mean temperature, isothermality, mean temperature in warmest and coldest quarters, and mean precipitation in wettest and driest quarters. The bioclim variables are calculated from other variables within the set, meaning that multicollinearity is a concern . To determine which variables should be included, previous modeling efforts were consulted. Bio8 (mean temperature of wettest quarter), bio9 (mean temperature of driest quarter), bio18 (precipitation of warmest quarter), and bio19 (precipitation of coldest quarter) were omitted from the analysis due to artifacts in the data found in similar research . The WorldClim climate surfaces were calculated based on average observations from field sites taken between 1960 and 1990. Data were interpolated across regions with few observation sites. However, data coverage in the USA, and particularly the Southeastern portion, were extensive and needed only minimal interpolation, therefore reducing associated surface error .
Population density was accessed through the Socioeconomic Data and Applications Center (SEDAC), which makes global socioeconomic data publicly available in formats that can be used in ArcGIS. In this case, the Population Density Grid v.3 data were downloaded in ascii format. Data represented the population density of the USA in the year 2000 to match occurrence and climate data and were displayed as population count/land area with a unit of persons per square kilometer. The data were calculated and displayed at 2.5 arc minutes . Because Ae. aegypti is highly anthropophilic, population density data are very important to include in model development and comparisons.
To model poverty, a factor that strongly influences an area’s environmental conditions, a raster surface from the 2000 US Census was added. Data were displayed at 30 arc seconds (~ 1 km2) but resampled at the resolution of the population density grids. The poverty variable was “Proportion of population living below the poverty level” for the continental USA in 2000 . Poverty levels for individuals were defined as annual incomes at or below $8350 and $17,050 for a family of four. Also, taken from the US Census was a regional shapefile for the Southeast to use as the extent and a state boundaries shapefile for visual reference [23, 24].
The Maxent model was chosen for this analysis because it is designed for presence-only data and known to outperform other common presence-only models . The method uses presence-only data to identify cells in which the species has been found. It then samples the user-provided predictors, in this case climate and human population characteristics, to identify other areas that match the presence areas of these predictors. Since occurrences that are geographically close are likely to exhibit spatial autocorrelation, Maxent requires rarefication (reducing the number of occurrences to no more than one per cell) prior to model implementation and the resulting output is a relative probability surface. It is important to remember that Maxent output should be interpreted as the probability of a location containing a suitable niche, not the probability of species presence [26, 27]. Finally, while global modeling efforts are common, Maxent may perform better when smaller geographical regions are modeled , which is particularly relevant in this paper.
The Maxent software version 3.3.3k was used for all modeling efforts. This is a GUI that requires the user to select the species training file and the environmental layers. Users are given several options for model output, and we chose to create response curves and perform jackknife analysis to test variable importance. Additionally, the file for the 20% testing dataset was specified and the model threshold was set at 10 percentile training presence, meaning that models which excluded more than 10% of training observations were rejected . Additional model settings comprised the inclusion of occurrence data within the background sampling process. This process replaces the uniform background data that is randomly sampled and instead creates a background dataset drawn from the distribution of the occurrence data . Maximum entropy distribution is then selected relative to the provided background, effectively eliminating selection bias. In addition to response curves, variable jackknifing, and variable gain values, Maxent outputs a probability surface ranging from 0.0 to as high as 1.0.
The Maxent model was run in four iterations. The first was climate only, as a reference for models with additional parameters. Then, models were run with climate and poverty, and again with climate and population density to assess how the addition of each predictor individually affected the model. The final model included the original climate variables and both poverty and population density. Model evaluation depended on the AUC, omission error, and commission. Omission error was calculated based on omission of test data (20% excluded from the model) at the 10-percentile training presence level, so once a model was produced based on the inclusion of at least 90% of training data, the test data were overlaid and all data that were found in areas of 0.5 logistic probability or less (< 50% likelihood of suitability) were deemed omitted. Commission was calculated as the percent of total study area predicted to have a probability greater than 0.5, also at the 10-percentile training presence level. The selection of a threshold is a topic of much debate (e.g., [30, 31]), and for our study, we selected a threshold (0.5) that maximizes the percent of points correctly classified (> 90%). This is a modification of the lowest presence threshold method . Models were ranked by order of AUC, so long as AUC was higher than 0.6 (0.5 means a model’s predictions are effectively random) . Individual components were also evaluated under jackknife and response curve results to determine how much each predictor contributed to the model. Additionally, AUC scores were compared between models to test for significant increases in improvement. To test for statistical significance, standard errors were obtained for each model (based on number of samples and model standard deviation), then a two-sample z-statistic was calculated using the difference of AUC values and standard errors between each model pairing.
Model fit statistics for all four models
Model 2: climate + poverty
Model 3: climate + population density
Full model: climate + population density + poverty
Population density was also the single largest contributor in the full model at 43.7%. Poverty was fourth at 11%, behind bio1 (mean annual temperature) at 15.4%, and bio6 (minimum temperature of the coldest month) at 14.6%. The jackknife analysis for the training set of occurrences showed that population density contributed the most useful information to the model, followed by the poverty variable.
Tests for significant improvement between each model
Mean AUC differences
Climate only and climate + poverty
Climate only and climate + pop. density
Climate only and full model
Climate + poverty and climate + pop. density
Climate + poverty and full model
Climate + pop. density and full model
The model fit indices for these models clearly demonstrate that the addition of poverty and population density improves the model’s predictive power. Additionally, while each predictor performs well when added independently, the best fitting and most accurate model is the full model, with both human-interaction predictors included. Investigating variable contributions and response curves supports this finding. The variable that contributes the most information varies depending on which model fit index is being discussed, meaning that both make valuable contributions to the overall model. While human-interaction predictors proved to be highly influential for the models and AUC scores showed significant differences, AUC scores alone were fairly high for all models. This indicates that AUC scores alone provide an incomplete picture, as the modeled surfaces appeared to be quite different between the climate-only model and the human-interaction models.
These findings support the assertion that models which only explore climate variables, such as temperature and precipitation, are incomplete and may be missing a significant source of information on habitat suitability. This is particularly highlighted in Maxent when considering that population density contributed almost three times more information than the next predictor, mean annual temperature. Since Ae. aegypti is highly anthropophilic in its feeding preferences and relies on standing water for breeding, it is logical that human population characteristics have proven to be important in modeling possible habitats and high-risk areas. While much of the Southeast United States may present a suitable habitat for Ae. aegypti, this study identified preferred habitat for the mosquito, and thus areas where Ae. aegypti is more likely to spread various diseases, including Zika. Variables used for this study, especially poverty and population density, may infer an observation/sampling bias for mosquitoes, but the mosquito locality database  applied both spatial and temporal methods of standardization to reduce sample bias, and our own process of rarefication further reduced any sampling bias that may be more related to abundance (i.e., number of observations in densely populated areas). Additionally, the union of occurrence data with background samples eliminated bias associated with random background sampling, resulting in a Gibbs distribution that accounted for bias in the final predictions .
Interestingly, the addition of human population characteristics illustrates that probability surfaces become more a probability of human-mosquito interaction that leads to niche realization than a probability of suitability or presence. This is an important distinction in understanding the risk mosquitoes pose to human health, particularly when seeking to implement mosquito control programs in regions that are at the highest risk of mosquito-borne illnesses. A high probability of suitability in these maps still does not guarantee species presence, but it does show that there is a high likelihood that the combination of ideal climatic conditions and human population characteristics provides suitable conditions for feeding and breeding. In many of the northern states within the study region (e.g., Kentucky, Tennessee, Oklahoma), occurrences of Ae. aegypti have thus far been minimal, but models indicate that many elements of a suitable habitat, and consequently a higher risk of mosquito-borne illnesses, are present, especially in the larger urban environments. Long-term changes in climate in these regions may result in an increased risk to mosquito-borne illnesses as mosquitoes establish populations at higher latitudes .
This research also indicates that, as Eisen and Moore  feared, there is the possibility of human population characteristics confounding the relationship between climate and mosquito presence. This is important considering the on-going focus on modeling both current species presence and changes in habitats due to climate change. As climate change alters the environmental landscape, it will also affect the human landscape . Additionally, human population dynamics will change on their own, further altering the manmade landscape, in terms of population characteristics and built environment . Failing to account for these changes may likely result in only a partial understanding of future species expansion, especially since so much of the current model’s explanatory power is due to population density.
Exploring human characteristics in climate projections currently faces some problems. For example, WorldClim projects climate change to the year 2050 but SEDAC only projects population density to 2030. Another issue is that poverty data were obtained from the US Census and not every country has data that is so complete or reliable. Therefore, adding human population characteristics to global maps or projections is problematic and requires data that may be difficult to obtain. Additional environmental parameters may be useful and could be compared to models that used population density (e.g., an enhanced vegetation index (EVI)). The EVI, or similar satellite-derived variables, may serve as acceptable proxies to population density.
An important limitation to consider in assessing these models is the lack of literature on which human population characteristics are most important. Poverty was included as a surrogate for likelihood of being in an environment that is suitable for mosquitoes (i.e., has untreated standing water) and population density was included since Ae. aegypti is highly anthropophilic. However, there may be more suitable surrogates for these characteristics or, for areas where it is available, measuring the presence of standing water may be a much better predictor. It is clear from these results that significantly more research is needed to better understand which human characteristics increase niche suitability (or preference) for vector mosquitoes and how to model those interactions.
Mosquito distributions have been modeled and predicted globally, regionally, and locally using many different methods and datasets; however, applying the concept of “human risk” to infection has not been well studied. In many cases, mosquito populations depend on local human populations and our inclusion of human population density in distribution models resulted in risk probability surfaces, not species distribution probabilities. This is an important distinction in understanding the models and is crucial for niche realization and subsequent risk of infection for humans. While the models (Figs. 4, 5, and 6) not surprisingly follow population density patterns, many cities, especially those in the northern extent of the study area (e.g., Nashville and Oklahoma City), show a higher risk probability when compared to the climate-only model (Fig. 3). Further research is needed to determine if additional human characteristics should be evaluated for model inclusion.
We would like to thank Drs. Megan Quinn and Ying Li for their editing and advice during the study design and writing phases.
Availability of data and materials
The datasets generated during and/or analyzed during the current study are available in the Dryad Digital Repository, http://datadryad.org/resource/doi:10.5061/dryad.47v3c.
JFO processed the data, developed the models, and wrote or co-wrote all sections of the manuscript; TAJ aided in the project design, data processing, model development, writing, and editing for all sections; JBH contributed to the project development and overall writing/editing. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Pratt HD, Moore CG. Mosquitoes of public health importance and their control, Insect control series. Atlanta: Communicable Disease Center (U.S.). Self-study Course 3013-G; 1960.Google Scholar
- Amaku M, Coutinho FA, Massad E. Why dengue and yellow fever coexist in some areas of the world and not in others? Biosystems. 2011;106:111–20.View ArticlePubMedGoogle Scholar
- Higgs S. Zika virus: emergence and emergency. Vector Borne Zoonotic Dis. 2016;16:75–6.View ArticlePubMedGoogle Scholar
- Kraemer M, Sinka M, Duda K, Mylne A, Shearer F, Barker C, et al. The global distribution of the arbovirus vectors Aedes aegypti and Ae. albopictus. elife. 2015;4:e08347.View ArticlePubMedPubMed CentralGoogle Scholar
- CDC. Dengue background. 2010. http://www.cdc.gov/dengue/epidemiology/index.html. Accessed 26 Apr 2016.Google Scholar
- CDC. Chikungunya virus: 2014 final data for the United States. 2015. https://www.cdc.gov/chikungunya/geo/united-states-2014.html. Accessed 28 Apr 2016.Google Scholar
- Runge JS. Lost trust: a yellow fever patient response. Yale J Biol Med. 2013;86:571–82.PubMedPubMed CentralGoogle Scholar
- PBS. Major American epidemics of yellow fever 1793–1905: American experience. 2006. http://www.pbs.org/wgbh/amex/fever/map/. Accessed 26 Apr 2016.Google Scholar
- CDC. Areas with Zika. 2016. https://www.cdc.gov/zika/geo/. Accessed 15 Oct 2016.Google Scholar
- Kraemer M, Sinka ME, Duda KA, Mylne A, Shearer FM, Brady OJ, et al. The global compendium of Aedes aegypti and Ae. albopictus occurrence. Sci Data. 2015;2:150035.View ArticlePubMedPubMed CentralGoogle Scholar
- Lozano-Fuentes S, Hayden MH, Welsh-Rodriguez C, Ochoa-Martinez C, Tapia-Santos B, Kobylinski KC, et al. The dengue virus mosquito vector Aedes aegypti at high elevation in Mexico. Am J Trop Med Hyg. 2012;87:902–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Christophers SR. Aedes aegypti (L.) the yellow fever mosquito: its life history, bionomics and structure. Cambridge: Cambridge University Press; 1960.Google Scholar
- Garcia-Rejon JE, Lorono-Pino MA, Farfan-Ale JA, Flores-Flores LF, Lopez-Uribe MP, Najera-Vazquez MR, et al. Mosquito infestation and dengue virus infection in Aedes aegypti females in schools in Merida, Mexico. Am J Trop Med Hyg. 2011;84:489–96.View ArticlePubMedPubMed CentralGoogle Scholar
- Eisen L, Moore CG. Aedes (Stegomyia) aegypti in the continental United States: a vector at the cool margin of its geographic range. J Med Entomol. 2013;50:467–78.View ArticlePubMedGoogle Scholar
- Jansen CC, Beebe NW. The dengue vector Aedes aegypti: what comes next. Microbes Infect. 2010;12:272–9.View ArticlePubMedGoogle Scholar
- Campbell LP, Luther C, Moo-Llanes D, Ramsey JM, Danis-Lozano R, Peterson AT. Climate change influences on global distributions of dengue and chikungunya virus vectors. Philosophical transactions of the Royal Society of London B Biol Sci. 2015;370.Google Scholar
- Brady OJ, Golding N, Pigott DM, Kraemer M, Messina JP, Reiner Jr RC, et al. Global temperature constraints on Aedes aegypti and Ae. albopictus persistence and competence for dengue virus transmission. Parasit Vectors. 2014;7:1–17.View ArticleGoogle Scholar
- Dawson IG, Johnson JE. Growing pains: how risk perception and risk communication research can help to manage the challenges of global population growth. Risk Anal. 2014;34:1378–90.View ArticlePubMedGoogle Scholar
- Reinmann AB, Hutyra LR, Trlica A, Olofsson P. Assessing the global warming potential of human settlement expansion in a mesic temperate landscape from 2005 to 2050. Sci Total Environ. 2016;545-546:512–24.View ArticlePubMedGoogle Scholar
- Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25:1965–78.View ArticleGoogle Scholar
- CIESN. Gridded Population of the World, version 3 (GPWv3): Population Density Grid. 2005. Palisades, NY: NASSA Socioeconomic Data and Applications Center (SEDAC). https://doi.org/10.7927/H4XK8CG2. Accessed 29 Jan 2017.
- Seirup L, Yetman G. U.S. Census Grids (Summary File 3), 2000. Palisades, NY: NASA Socioeconomic Data and Applications Center (SEDAC); 2006. https://doi.org/10.7927/H42R3PMN. Accessed 29 Jan 2017.Google Scholar
- US Census Bureau. Cartographic boundary shapefiles-regions. 2014. https://www.census.gov/geo/maps-data/data/cbf/cbf_region.html. Accessed 15 Jan 2017.
- US Census Bureau. Cartographic boundary shapefiles-states. 2014. https://www.census.gov/geo/maps-data/data/cbf/cbf_state.html. Accessed 15 Jan 2017.
- Elith J, Graham C, Anderson R, Dudík M, Ferrier S, Guisan A, et al. Novel methods improve prediction of species’ distributions from occurrence data. Ecography. 2006;29:129–51.View ArticleGoogle Scholar
- Merow C, Smith MJ, Silander J. A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. Ecography. 2013;36:1058–69.View ArticleGoogle Scholar
- Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Model. 2006;190:231–59.View ArticleGoogle Scholar
- Franklin J. Mapping species distribution. Cambridge: Cambridge University Press; 2010.View ArticleGoogle Scholar
- Phillips S, Dudik M. Modeling of species distributions with Maxent: new extensions and a comprehensive evaluation. Ecography. 2008;31:161–75.View ArticleGoogle Scholar
- Bean WT, Stafford R, Brashares JS. The effects of small sample size and sample bias on threshold selection and accuracy assessment of species distribution models. Ecography. 2011;35:250–8.View ArticleGoogle Scholar
- Escalante T, Rodríguez-Tapia G, Linaje M, Illoldi-Rangel P, González-López R. Identification of areas of endemism from species distribution models: threshold selection and Nearctic mammals. Revista Especializada en Ciencias Químico-Biológicas (TIP). 2013;16:5–17.Google Scholar
- Pearson RG, Raxworthy CJ, Nakamura M. Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar. J Biogeogr. 2007;34:102–17.View ArticleGoogle Scholar
- DeLong ER, DeLong DM, Clarke-Pearson DL. Comparing the areas under two or more correlated receiver operating characteristic curves: a nonparametric approach. Biometrics. 1988;44:837–45.View ArticlePubMedGoogle Scholar
- Hanley JA, McNeil BJ. A method of comparing the areas under receiver operating characteristic curves derived from the same cases. Radiology. 1983;148:839–43.View ArticlePubMedGoogle Scholar