A comparative modeling study on non-climatic and climatic risk assessment on Asian Tiger Mosquito (Aedes albopictus)

Study area selection

The study site is located between 75°30′00″W and 92°00′00″W, and 25°00′00″N and 36°30′00″N in USA (Fig. 1). In selecting the study site, we looked for an area exhibiting variations of each conditioning factor, as well as Ae. albopictus presence. For example, in terms of altitude, the study area should display a range of altitudes. Our selected study area had an altitude range from 0 m to 2,031 m above sea level. For geology, the study area had 80 different geological fractures.

Spatial datasets

Conditioning factors

The six geographical conditioning factors (i) altitude, (ii) slope, (iii) aspect, (iv) distance of locality from road, (v) distance of locality from river and (vi) geology, with a grid cell size 90 × 90 m were used to run the EBF model. The quantile classification scheme was used for all conditioning factors, as recommended by Tehrany, Pradhan & Jebur (2013). Altitude, slope and aspect layers were generated from DEM data obtained from EarthExplorer (2014) as shown in Figs. 2A–2C respectively. Distances from road and river layers were generated by Euclidean Distance tool and divided into ten classes using the quantile method, as shown in Figs. 2D and 2E respectively. The geology layer, obtained from the United States Department of Agriculture (USDA) (United States Department of Agriculture, 2017) contained 80 different types of lithology as shown in Fig. 3. The elevation layer was included as it depicts climate variations and the physical barriers limiting dispersion. The road, river and geological layers were included as the greatest densities of Ae. albopictus occur in urban environments (Rochlin et al., 2013a).

Non-climatic modeling

Evidential Belief Function (EBF), which is also called Dempster-Shafer theory of evidence, was developed by Dempster (Dempster, 1967), based on the Bayesian theory of subjective probability. Its advantages are the relative flexibility with which it accepts uncertainty and its ability to aggregate beliefs from many sources of evidence (Thiam, 2005). Rather than estimating the validity of probabilities, the Dempster-Shafer technique calculates the nearness of the evidence in proving the validity of a hypothesis (Pearl, 1990). Applications of EBF have been effective in many fields of research that utilize GIS data (Malpica, Alonso & Sanz, 2007).

To produce a hazard index of presence of Ae. albopictus, the conditioning factors were expressed individually as acquired weights and then aggregated (Eq. (1)).

Assuming a set of Ae. albopictus presence conditioning factors C = (C_i, i = 1, 2, 3, …, n), consisting of mutually exclusive and exhaustive factors C_i. C represents the frame of discernment, and a fundamental probability assignment is represented by the function m:P(C) → [0, 1].

The set P(C) includes all subsets of C, as well as C itself and the empty set. m:P(C) → [0, 1] is described as a mass function and satisfies m(Φ) = 0 and ∑_ACm(A) = 1, in which Φ represents the empty set and A represents any subset of C. The m(A) estimates to what degree the evidence supports A, which is denoted by Bel (A), a belief function.

There are four basic evidential belief functions attributable to a proposition, based on evidence. These four functions establish the degree of: (i) Belief (Bel), (ii) Disbelief (Dis), (iii) Uncertainty (Unc) and (iv) Plausibility (Pls). Bel represents the lower bound and Pls represents the upper bound of probability (Althuwaynee, Pradhan & Lee, 2012; Awasthi & Chauhan, 2011). Unc is established by the difference between Bel and Pls, and represents the ignorance. Dis represents the degree of probability that the proposition is false.

Dis = 1 − Pls or 1 − Unc − Bel, such that Bel + Unc + Dis = 1. For a case of C_ij zero presence of Ae. albopictus, implying that Bel = 0, Dis is reset to zero, whether that is the case or not (Carranza, Hale & Faassen, 2008). EBF can be estimated on the basis of a subjective judgment or calculated on the input of data (Srivastava, Mock & Gao, 2011). By superimposing the inventory map (L) of Ae. albopictus onto the six individual conditioning factor maps, we ascertained the number of pixels with Ae. albopictus presence and absence, for each separate conditioning factor. Assuming that N(L) represents the total of presence pixels and N(C) the total pixels comprising the study site, C_ij represents the jth class attribute of Ae. albopictus presence conditioning factors C_i ( i = 1, 2, …, n), N(C_ij) is the total of pixels for class C_ij, and N(L∩C_ij) is the Ae. albopictus presence pixels in C_ij. According to (Carranza & Hale, 2003), estimation of EBFs based on data is represented by: (1)${\displaystyle Bel\left(C_{ij}\right)=\frac{W_{C_{ij}\left(Ae.Albopictus\text{presence}\right)}}{\sum _{j=1}^n{W}_{C_{ij}\left(Ae.Albopictus\text{absence}\right)}}}$ (2)${\displaystyle W_{C_{ij}\left(Ae.Albopictus\text{presence}\right)}=\frac{\frac{N\left(L\cap C_{ij}\right)}{N\left(L\right)}}{\frac{\left[N\left(C_{{}_{ij}}\right)-N\left(L\cap C_{ij}\right)\right]}{\left[N\left(C\right)-N\left(L\right)\right]}}}$ (3)${\displaystyle Dis\left(C_{ij}\right)=\frac{W_{C_{ij}\left(Ae.Albopictus\text{absence}\right)}}{\sum _{j=1}^n{W}_{C_{ij}\left(Ae.Albopictus\text{absence}\right)}}}$ where, (4)${\displaystyle W_{C_{ij}\left(Ae.Albopictus\text{presence}\right)}=\frac{\frac{\left[N\left(C_{ij}\right)-N\left(L\bigcap C_{ij}\right)\right]}{N\left(L\right)}}{\frac{\left[N\left(C\right)-N\left(L\right)-N\left(C_{ij}\right)+N\left(L\bigcap C_{ij}\right)\right]}{\left[N\left(C\right)-N\left(L\right)\right]}}.}$ In Eq. (2) the numerator represents the proportion of Ae. albopictus presence pixels occurring in factor class C_ij, while the denominator represents the proportion of Ae. albopictus absence pixels in factor class C_ij. In Eq. (4) the numerator represents the proportion of Ae. albopictus absence pixels in factor class C_ij, while the denominator represents the proportion of Ae. albopictus absence pixels in attributes excluding factor class C_ij. The parameter W_Cij (Ae. albopictus presence) represents the weight of C_ij supporting the belief that Ae. albopictus presence exceeds Ae. albopictus absence. Parameter W_Cij (Ae. albopictus absence) is the weight of C_ij that supports the belief that Ae. albopictus absence exceeds presence.

After calculation of the EBF function for all Ae. albopictus presence conditioning factors, Dempster’s combination rule was introduced to produce the four integrated EBFs (Dempster, 1967). The formulae for combination of two Ae. albopictus presence conditioning factors _C1 and _C2 are as follows (Carranza, Woldai & Chikambwe, 2005): (5)${\displaystyle {\mathrm{Bel}}_{C_1{C}_2}=\frac{{\mathrm{Bel}}_{C_1}{\mathrm{Bel}}_{C_2}+{\mathrm{Bel}}_{C_1}{\mathrm{Unc}}_{C_2}+{\mathrm{Bel}}_{C_2}{\mathrm{Unc}}_{C_1}}{1-{\mathrm{Bel}}_{C_1}{\mathrm{Dis}}_{C_2}-{\mathrm{Dis}}_{C_1}{\mathrm{Bel}}_{C_2}}}$ (6)${\displaystyle {\mathrm{Dis}}_{C_1{C}_2}=\frac{{\mathrm{Dis}}_{C_1}{\mathrm{Dis}}_{C_2}+{\mathrm{Dis}}_{C_1}{\mathrm{Unc}}_{C_2}+{\mathrm{Dis}}_{C_2}{\mathrm{Unc}}_{C_1}}{1-{\mathrm{Bel}}_{C_1}{\mathrm{Dis}}_{C_2}-{\mathrm{Dis}}_{C_1}{\mathrm{Bel}}_{C_2}}}$ (7)${\displaystyle {\mathrm{Dis}}_{C_1{C}_2}=\frac{{\mathrm{Dis}}_{C_1}{\mathrm{Dis}}_{C_2}+{\mathrm{Dis}}_{C_1}{\mathrm{Unc}}_{C_2}+{\mathrm{Dis}}_{C_2}{\mathrm{Unc}}_{C_1}}{1-{\mathrm{Bel}}_{C_1}{\mathrm{Dis}}_{C_2}-{\mathrm{Dis}}_{C_1}{\mathrm{Bel}}_{C_2}}.}$

Integrated EBFs of the Ae. albopictus presence conditioning factors are applied in sequence by means of Eqs. (5)–(7). Table 1 shows the estimated EBFs for the six Ae. albopictus presence conditioning factors.

Climatic data, future scenarios and climate models

Baseline climate was represented by the WorldClim current climate dataset of BIOCLIM variables (http://www.worldclim.org). WorldClim is a high-resolution climate average for the period 1961 to 1990, with global coverage and spanning the time period over which the majority of occurrence records were collected. Possible future climates at global scale incorporated four IPCC5 greenhouse gas concentration (GHC) trajectories, which differ in terms of GHC emission peaks. The lower the number of the trajectory, the earlier in the century it peaks.

We purposefully chose the worst (extreme) RCP8.5 (peak 2080) (Stocker et al., 2013) for incorporation into the future climate scenario in the model projections as it is not yet possible to determine which estimates of the climate change and RCPs of 2.6, 4.5, 6.0 and 8.5 are the most reliable (Randall et al., 2007). RCP8.5 is a representative concentration pathway that includes relatively high emissions of greenhouse gases. Other factors assumed in RCP8.5 are high demographic development, relatively slow economic growth, with modest progress in technology and the introduction of novel sources of energy. These factors culminate in an increased demand for energy and higher GHG emissions over the long term, without a more radical approach to the projected impact of climate change (Riahi et al., 2011).

There are 19 General Circulation Models (GCMs) in WorldClim database and we have selected GCMs of Miroc3.2 and CSIRO-MK30, which have higher reputations and have been used for projections of many invasive species, agricultural crops and pests (Da Silva et al., 2017a; Da Silva et al., 2017b; Lamsal et al., 2017; Paterson et al., 2017; Ramirez-Cabral, Kumar & Shabani, 2017b; Ramirez-Cabral, Kumar & Shabani, 2017a; Shabani, Kumar & Ahmadi, 2016; Shabani, Kumar & Ahmadi, 2017; Shabani, Kumar & Taylor, 2012).

Climatic modeling

MaxEnt desktop version 3.3.3k, with modified parameters, was used to construct the climatic model (Phillips, Anderson & Schapire, 2006; Phillips & Dudík, 2008). MaxEnt requires a user-defined background of geographical data (Guillera-Arroita, Lahoz-Monfort & Elith, 2014) in order to compare the climate of a set of grid cells representing the presence of a species, with the reference set representing the climate of the sampled cells. The selected geopgraphical data is a significant determinant of the results of the model (Elith et al., 2011) and it is important that it reflect all environmental variations covering the areas representing the presence of the species (Elith, Kearney & Phillips, 2010). The algorithm in MaxEnt estimates the maximum entropy probability distribution that approximates uniformity, based on a comparison of presence and background location interactions with a set of variables, limited by parameters imposed by the observed spatial distributions and environmental factors. Optimisation of the maximum entropy probability distribution is achieved by minimisation of the relative entropy between presence and background point data (Phillips, Anderson & Schapire, 2006). MaxEnt, with inbuilt MESS analysis tool, has the capacity to predict future distributions, generated from two datasets of environmental variables (Elith et al., 2011). In our study, the current conditions are used to generate the model, with a set of variables utilized for projection of the future scenario (in this case 2055).

Using jackknife analysis and Pearson correlation technique to correlate coefficients, we selected the most influential variables showing low correlation (R² < 0.5) for this modeling study. Here, BIO11 (Mean Temperature of Coldest Quarter), BIO16 (Precipitation of Wettest Quarter) and BIO17 (Precipitation of Driest Quarter) were selected for the modeling. To achieve greater consistency of background data and overcome the potential for finding fewer records representing areas more recently experiencing invasions, as well as those incompletely sampled, we assigned greater prominence to the records representing less geographical proximity. However, it should be noted that without information on actual survey returns, there is no method of separating unsuitable and under-sampled areas, and that the weighting of prominence cannot overcome the fusion of these two categories of data. After using the Gaussian kernel method to establish deviations from the ArcGIS default values, the formula applied for weighting is to divide total weighted records by the weighted number of land cells of the specific area, to exclude coastal region edge effects. By adjusting the resulting grid to a range of 1–20, extreme values were excluded. This weighting method, as advocated by Elith, Kearney & Phillips (2010), reduces bias that gives prominence to the records of more densely sampled areas. Background training points were generated from the kernel density layer for the species using Hawths Tools extension (Beyer, 2004).

Model validation

The non-climatic modeling analysis was executed and validated using known Ae. albopictus presence (Fig. 1). Using training and testing Ae. albopictus presence data, validation was carried out using the area under curve (AUC) method. While training presence data was used to generate the model, the results using this data does not fully represent the model’s total efficiency. The prediction rate was measured to establish how efficiently the model and selected conditioning factors predicted Ae. albopictus presence. AUC can assess prediction accuracy qualitatively by arrangement of the calculated values of all cells of the study locations into descending order, providing an individual hierarchical ranking of the accuracy of each prediction. Thereafter, the values of cells were divided into 100 classes with accumulation intervals of 1%.

Non-climatic modeling

We examined, and assessed individually, six geological variable factors that directly impact the presence of Ae. albopictus in a specific locality. The altitude EBFs indicated that localities of 0 to 97 m above sea level had a high probability presence of Ae. albopictus. The belief value (Bel) peaked at 24 with altitudes from 15 to 27 m, while it was 5 and 0 at altitudes from 127 to 167 m and 167 to 228 m respectively (Table 1). Slope EBF indicated that classes of 2.72 to 7.84°, 10.56 to 13.86°and 33.76 to 76.85°produced Bel values of 12, 14 and 4, respectively. The three highest Aspect Bel values of 15, 13 and 12 related to the classes of South, Northwest and Southwest, respectively (Table 1). Distance from locality to road was included as a conditioning factor, as motor vehicles have been shown as a means of Ae. albopictus transmission. The highest Bel values for this factor were 22 and 19, representing the classes of 0 to 252 m and 252 to 757 m, respectively. For Distance of locality from river, EBF estimated the probability of Ae. albopictus presence, for all ten classes (Table 1). For Geology, the classes of unconsolidated deposit, calcarenite and delta scored Bel values of 16 and 14 and 41, indicating the probability of Ae. albopictus presence in these geological formations (Table 1).

Probability index and suitability map

The verification results for EBF model are shown in Fig. 4. Probability index maps of Ae. albopictus presence produced by EBF method are shown in Figs. 5A and 5B respectively. The range is from 0 to 1, where 0 represents zero probability and 1 represents 100% probability. For producing suitability maps, as well as improving the visual interpretation of locational suitability, probability maps require some form of classification (Umar et al., 2014). There are a variety of classification techniques such as equal interval, natural break, standard deviation and quantile, the selection of which should be based on the research data characteristics and study objectives. Equal interval is suitable when the data displays a normal distribution, while standard deviation arranges the data into a fixed number of classes. Natural break suits a dataset exhibiting a sudden or big jump. Here, in order to have a reliable judgment regarding the impact of every class of each factor on species occurrence, we attempted to reduce the influence of classification algorithm on the conditioning factors classes as much as possible. In some population analysis projects where the goal is to find a big jump in the data, natural break technique is highly recommended (Hui et al., 2010; Umar et al., 2014), while in this research, this method would not be efficient. Hence, quantile-based classification technique was found to be more appropriate to classify the factors in this study. This method groups equal number of pixels (area) into each group without any interference in the separation of the pixels. We thus selected the quantile method to produce the suitability classes. The verification results for EBF model are shown in Fig. 4. The AUC results showed 0.73 success rate and 0.70 prediction rate (Fig. 4) and these values are high enough and satisfactory for model prediction as documented in Umar et al. (2014). Our probability indexes were into five zones of suitability: very low, low, moderate, high, and very high, for EBF output (Fig. 5B).

Climatic modeling

The climatic model produced by MaxEnt, using two GCMs, Miroc3.2 (Fig. 6A) and CSIRO-MK30 (Fig. 6B), under the RCP 8.5 scenario, shows virtually the whole study site is highly suitable for Ae. albopictus and that this condition will persist until at least 2055. Comparing the GCM projections for 2055, CSIRO-MK30 produced a more moderate pattern of climatic suitability than Miroc3.2. Both GCM response curves show the highest probabilities of Ae. albopictus presence in areas with Coldest Quarter Mean Temp (bio11) from 16 to 23 °C, Wettest Quarter Precipitation (bio16) of 430 mm and Driest Quarter Precipitation (bio17) of 350 to 450 mm (Figs. 6A and 6B). The Miroc3.2 mean AUC was 0.868, while CSIRO-MK30 indicated 0.864.