Establishment of an expansion-predicting model for invasive alien cerambycid beetle Aromia bungii based on a virtual ecology approach

The pragmatic management of invasive alien species should integrate two essential items: 1) management interventions and 2) a spatially explicit management plan. Predicting the future expansion of target species in a region at the early invasion stage is an important step toward the establishment of a spatially explicit management plan. However, information regarding the distributions of target species is limited, making it challenging to predict range expansions. In the present study, we established a simulation model that could predict the future expansion of the invasive insect Aromia bungii, which is harmful to Prunus trees (including cherry trees [Cerasus × yedoensis]), in Japan. We employed a virtual ecology approach that simulated species dynamics based on a simple model in Saitama Prefecture, which is in the Kanto region of Japan. Since the first record of the species in this region of Japan in 2013, its range has expanded dramatically. Three candidate pathways and combinations of these for the range expansion of A. bungii were tested to identify the major proxies of expansion for this species, followed by the validation of these results using occurrence records for the species through 2019. Both the river density model and combined river and road density models showed good predictive performance. Using these models, we established a predictive map of the future expansion of this species in the wider range of the simulation area. Based on the results, we recommend concentration of management efforts in the mid-northeast region of the Saitama Prefecture.


Introduction
Invasive alien species can cause environmental and economic losses, which have been well-documented (Pimentel et al. 2001;Pyšek and Richardson 2010;Vilà et al. 2010). For example, Diagne et al. (2021) estimated that the annual mean cost of biological invasions worldwide could reach US$162.7 billion as of 2017. To minimize such losses, substantial time and effort must be dedicated to managing populations of invasive alien species and preventing their re-introduction and further spread (Pichancourt et al. 2012;Osawa et al. 2016Osawa et al. , 2019. Although studies have been conducted with respect to management strategies for invasive alien species, pragmatic management should integrate two essential items: 1) management interventions, i.e., the establishment of specific actions for the target species, and 2) spatially explicit management plans, i.e., the prioritization of target areas within and/or potentially within the invasion range of the target species (Foxcroft et al. 2009;Osawa 2020). An adequate management plan would lead to the optimal allocation of limited resources across target areas (Januchowski-Hartley et al. 2011;McDonald-Madden and Chades 2011;Kumschick et al. 2012;Osawa et al. 2019). Hence, land managers and/or practitioners need to establish strategic management plans against invasive alien species based on a combination of specific, effective actions and spatial prioritization plans (Foxcroft et al. 2009;Osawa 2020). After establishing a strategic management plan, efforts are most effective when invasion is detected at an early stage and comprehensive control measures are rapidly implemented (Simberloff 2003;Kaplan et al. 2014).
The longhorn beetle Aromia bungii (Coleoptera: Cerambycidae) is a pest insect that is native to eastern Asia from Vietnam to Mongolia and the Russian Far East (Kong 2015;Zou et al. 2019); however, it is not native to Japan because of Japan's geographic isolation. Since its first adult record in 2011 (Adachi 2017) and first adult emergence record in 2012 (Kiriyama et al. 2015) in Japan, this species has caused substantial damage to peach orchards and a species of cherry tree (Cerasus × yedoensis "Somei-yoshino") that is one of the most popular ornamental trees in the country. In response to these serious problems, several studies have addressed various aspects of this species, including their basic ecology (Iwata 2018) and the establishment of methods to eradicate them (physical [Nakano 2018], biological [Sunamura et al. 2020], and chemical [Xu et al. 2017;Yasuoka 2017] approaches) in recent years. To apply these established management actions effectively, a spatially explicit management plan is required for the allocation of limited resources. However, based on our knowledge, there are few or no studies on the establishment of a spatially explicit management plan for this species.
One of the most effective approaches for establishing a spatially explicit management plan is the prediction of future expansion on the basis of ecological traits of the target species (Ficetola et al. 2007;Fukasawa et al. 2009;Koike and Iwasaki 2011;Osawa and Ito 2015). Although species distribution models (SDMs) are frequently applied for predicting potential future expansion of species, these types of techniques assume that species are at equilibrium with their environments, thus there are serious challenges to apply such techniques for range-shifting species (Menke et al. 2009;Elith et al. 2010). Gallien et al. (2010) indicated that mechanistic simulation models are the tools of choice to predict range expansion because they can explicitly incorporate ecological processes and dynamics. In fact, to estimate geographical range expansions, several simulation models, from complex to simple, have been developed (Gallien et al. 2010;Koike and Morimoto 2018). However, in the early stages of an invasion, such as that of A. bungii in Japan, ecological information regarding the target species in target areas is severely limited (Koike 2006), making it necessary to predict its range expansions based on limited data (Koike 2006;Osawa and Ito 2015;Koike and Morimoto 2018). Under such conditions, the virtual ecology approach, which simulates species dynamics using simple models together with "virtual observation" of the species dynamics (Zurell et al. 2010;Pagel and Schurr 2012), is an effective way to predict species expansion (Skarpaas et al. 2005;Osawa and Ito 2015;Osawa et al. 2016Osawa et al. , 2018. One of the merits of this approach is that it can apply any model type from complex to simple (Zurell et al. 2010); thus, it can be used even with simple models that do not require detailed ecological information. Using this approach, a hypothetical simulated dispersal can be evaluated using actual distribution records, even with limited data (Osawa and Ito 2015;Osawa et al. 2016Osawa et al. , 2018. In the present study, a virtual ecology approach was used to predict the future expansion of the pest A. bungii, a winged insect with flying ability (Iwata 2018). Using a simple concept, we predicted that this species expands based on habitat continuity, particularly using its host trees (i.e., Prunus and/or Cerasus) as a food resource at the larval stage. The dispersal process essentially depends on local factors, such as habitat quality (Chen et al. 2008). We used a cellular automata (CA)-based virtual ecology model established previously by Osawa et al. (2018) for the mirid bug (Stenotus rubrovittatus), a winged insect with flying ability similar to A. bungii (Osawa et al. 2018). Although S. rubrovittatus and A. bungii are not related to one another, the model can be applied to both because of their common trait of flying ability. In this model, a population can expand into a neighboring landscape matrix per the resistance of the landscape elements in the simulation. In this case, the resistance is defined as the cost of dispersing across a particular landscape feature relative to the cost of dispersal across a reference condition (McRae 2006;Graves et al. 2014). Thus, we defined landscape features that contained abundant food resources to have low resistance, i.e., high habitat quality, in this model (Osawa et al. 2018). In our model, invasion patterns of the target species depend only on the spatial locations of matrix types; thus, the model does not require fitting species-specific dispersal parameters. Using this model, we simulated the dispersal process of A. bungii with the help of several candidate proxies: river density, road density, and the number of planted cherry trees, which could reflect the amount of food resources and habitat because of the distribution of cherry tree-lined areas (described in Methods). We subsequently evaluated their simulation output results using real distribution records to identify the effective proxies of expansion for this species among candidate proxies, and combined these. On the basis of this approach, we predicted an area with a high invasion risk over a wide area using the model with the most important proxies identified.

Study species
The invasive wood-boring beetle A. bungii (Faldermann) (Coleoptera: Cerambycidae, subfamily Cerambycinae, tribe Callichromatini) is native to eastern Asia from Vietnam to Mongolia and the Russian Far East (excluding the Japan Archipelago); however, it has recently become established in Japan, Germany, and Italy, and it has also been intercepted in cargo entering the USA, Great Britain, and Australia (Ma et al. 2007;Garonna et al. 2013;Zou et al. 2019). In Germany and Italy, eradication protocols are in effect, including the removal of infested trees, public education campaigns, and intensive surveillance around any known sites of infestation (Kong 2015). While there are no known established populations in the latter three countries as of yet, detection and quarantine measures have been initiated to prevent their establishment (Zou et al. 2019).
Aromia bungii infests trees of the genus Prunus, which includes a number of economically important fruits, including peaches, plums, cherries, and apricots (Xu et al. 2017), as well as ornamental cherry trees. The feeding-associated damages of the vascular tissues and weakening of the trunk and branches by larval tunneling lead to the death of the host tree over several years. The Center for Agriculture and Bioscience International (CABI, Invasive Species Compendium, https://www.cabi.org/isc/datasheet/ 118984#EF51C788-F9D6-4F79-B667-F7DD03B02F65, accessed 1 July 2021) lists A. bungii as a significant risk to all stone fruit-growing countries in Europe and neighboring countries (CABI, Invasive Species Compendium, https://www.cabi.org/isc/datasheet/118984#EF51C788-F9D6-4F79-B667-F7DD03B02F65, accessed 1 July 2021).

Study area
The study area was comprised of the northern section of Saitama Prefecture, including parts of Gunma and Tochigi Prefectures, Japan ( Figure 1). The area has a mean annual precipitation of 1,346 mm and mean annual temperature of 14.8 °C (Japan Meteorological Agency; Saitama City. https://www.data.jma.go.jp/obd/stats/etrn/view/nml_amd_ym.php?prec_n o=43&block_no=0363&year=&month=&day=&view=, accessed 1 July 2021). The first detection of damage on cherry trees by A. bungii in Saitama Prefecture was recorded in 2013 in its southern area (Shoda-Kagaya 2015), followed by identification of damage in 2017 in the northern areas (EKTDI 2018). The distribution range of this species has since expanded drastically, with serious damage to Prunus species in this area.

Figure 1. Study area and simulation unit
We used a grid size of approximately 1 km (termed 1 km cell hereafter) as a study unit (Figure 1). The 1 km cell grid system is a standard Japanese unit utilized for several types of statistics (Statistics Bureau of Japan). This unit is also consistent with the potential annual expanded range of A. bungii reported previously in Japan; specifically, new larval frass in a cherry tree (Cerasus × yedoensis) was found about 1 km from a cherry tree that was infected the prior year in Saitama Prefecture (Shoda-Kagaya 2015). The total number of 1 km cells in the target area was 1,140 ( Figure 1).

Insect data collection
The distribution records of A. bungii were obtained by the authors from local land managers, such as park managers and citizens, within the study area. From the latter, we requested A. bungii occurrence records of both captured adults and larval frass from trees. In the larval stage, infestations can be detected in later stages by the accumulation of larval frass (Shoda-Kagaya 2015; Xu et al. 2017). To request information from the general public, we published a press release on the activity and method of reporting. More specifically, we requested observation records including date and location information (such as address and latitude/longitude) for both adults and larval frass along with photographs containing contact information for the provider. To support these data collection activities, we published a brief manual on data collection and identification of the target species on the Internet (currently version 5 is available, https://www.pref.saitama.lg.jp/ documents/117809/kubiakamanual5.pdf, accessed 1 July 2021). After receiving this information, we confirmed all collected records of A. bungii based on both field surveys and photographs captured by the providers. If there were one or more occurrence records within the 1 km cell, we used that the cell as an occurrence cell in this study. The occurrence records through 2020 were released on the website of the Center for Environmental Science in Saitama (https://cessgis.maps.arcgis.com/apps/webappviewer/index.html? id=67495ef407974aa7bdcbfcb42749471b, accessed 1 July 2020).

Virtual ecology model
We utilized a simple CA model that can predict the expansion of a target species based on candidate habitats/pathways (Osawa et al. 2018). The digital space in the CA model comprises a rectangular grid of square cells that represents the target area ( Figure 1). The cell is set to the same size as the unit of the predicted range expansion. The model yields a theoretical number of invasions (described below) in each cell, which serves as an invasion probability. Each cell has a unique ID and parameter, which is the expansion path vector (Figure 2a). The cell ID indicates the spatial location of the cell within the area of analysis, and the expansion path involves four variables that represent the four directional vectors in adjacent cells: top, left, bottom, and right ( Figure 2b). Each variable is a probability value, i.e., 0-1, which is independent from the three other variables and indicates the probability of successful invasion for the adjacent cells (described below in detail). That value is determined by the matrix resistance value in this study.
The analysis unit of A. bungii was not an individual insect; rather, it was the local population in a grid cell. Our model simulates the capability for range expansion within the target area; thus, if a population invaded the neighboring cells, the occurrence cell number would increase and not just simulate a move, i.e., not disappear in the starting cell. A population can expand its range according to the expansion probability values and not according to population growth (described below in detail).

Expansion probability between cells
In this CA model, the expansion probability was defined as the probability of expansion into the adjacent cells: top, left, bottom, and right ( Figure 2b). The expansion probability was defined as a vector with four probability values; i.e., numerical sequence: where e1, e2, e3, and e4 indicate the probability of expansion success to the top, left, bottom, and right cells, respectively ( Figure 2b). These values are independent of one another. If all expansion path values are 0, the population in this cell cannot expand to any other cells. By contrast, if all expansion probability values are 1, the population in this cell can expand to all adjacent cells. We used three candidate variables to represent potential habitats/pathways (river density, road density, and the number of planted cherry trees) for setting the expansion probability values (Supplementary material Figure S1). We selected these three candidates with respect to the feeding habits of A. bungii at the larval stage; it feeds on various species of Prunus. Among Prunus species, the cherry tree (Cerasus × yedoensis) is the most widely distributed in Japan, including the study area. In general, there are lines of cherry trees growing along rivers (Kanebishi 2010), motorways, and footpaths (Kyoya et al. 2003) in Japan. Human movement along these features, such as walking or driving, could also support the expansion of A. bungii (Akasaka et al. 2015). The river line data were derived from the National Land Numerical Information Download Service of Japan (Ministry of Land, Infrastructure, Transport and Tourism, Japan). The river lines were derived from GIS line data on national primary rivers (Ministry of Land, Infrastructure, Transport and Tourism, Japan, accessed 1 July 2021). Although we did not consider smaller rivers (secondary rivers, small streams, or ditches), this is not a severe problem because large rivers in Japan are more likely to be cherry tree-lined than small rivers. Moreover, the road data were derived from the National Land Numerical Information Download Service of Japan. The data were derived from total road length per 1 km cell, which does not include other road characteristics (such as width). The number of planted cherry trees was derived from local governments and calculated per 1 km cell. This information included only trees in public facilities and not those in private sectors. These records have also been provided on the CESS website (https://cessgis.maps.arcgis.com/apps/View/index.html?appid= 78aef6d0168a4d178f6b248750bf8dd9, accessed 25 August 2021). Before the analysis, we tested the statistical correlations among the three variables using Pearson's method. The results showed that the river and road densities were significantly correlated (Table 1). However, all absolute correlation values were < 0.1 (Table 1), indicating that the simulation results would not be similar. The expansion probability was defined based on the amount of a candidate habitat/pathway per unit area. The probability values (e1, e2, e3, and e4) were defined as the amount of a candidate habitat/pathway in a destination cell divided by the maximum amount per cell in the target area. That value could change depending on the target area. Hence, these were relative values throughout the study area and reflect the spatial position of the landscape matrix in the study area. If the landscape matrix shows a sequence of low resistance; i.e., if there is an existing continuous habitat/pathway, then the insect population can expand via the pathway.

Expansion simulation
Simulation, validation and application procedures are briefly shown in Figure 3. First, we conducted simulations using three variables (river density, road density, and the number of planted cherry trees) that were run independently (Figure 3a). In addition, we conducted a simulation using a null variable; namely, all probability values (e1, e2, e3, and e4) for all cells were the same. The probability values of the null model were determined from the average values of probability of the river density model (0.27). We predicted and compiled the expansion ranges from all cells in the CA field; i.e., we started the simulation from one cell to its steady state, repeated that simulation for all cells, and summed them all, which we termed a "run". Thus, each "run" consisted of n simulations, where n is the number of cells in the CA field. In each iteration, we started with 0's in all cells except for one initial invasion site, which was set to 1. The initial invasion site differed in each simulation within the run. Each simulation ran for 999 steps, which is almost the same as the number of cells with the target area (1,140), to reach the steady state. At each step, each invaded cell attempted to invade the cells to the top, left, bottom, and right of itself based on these cells' invasion probabilities from equation (1). A successful invasion converted a cell from a value of 0 to 1. The next step was then initiated from any newly invaded cells and the previously invaded cells. Cells with high values indicated the possibility of invasion by virtual populations from several other cells. We used that value as our metric of expansion success. Second, we conducted simulations using a combination of expansion variables that showed high performance during independent simulation (Figure 3b). We conducted two types of combination models. In combined model 1, we defined the probability values (e1, e2, e3, and e4) by their product among variables that were weighted based on their performances during single variable simulation and normalized these. For example, the combination of river density (e_river1, e_river2, e_river3, and e_river4) and road density (e_road1, e_ road2, e_road3, and e_road4) was determined as follows: Expansion probability with river and road = w1(x_river) * w2(x_road) * N (2), where x_river indicates the probability values of (e1, e2, e3, and e4) on the river, and x_road indicate the probability values of (e1, e2, e3, and e4) on the road, respectively. Both w1 and w2 were weighted functions that determined the area under the curve (AUC) value (described later), which was used as a model performance index. Both w1 and w2 were determined as follows: w1 = AUC on river model/(AUC on river model + AUC on road model) (3) w2 = AUC on road model/(AUC on river model + AUC on road model) (4) If the AUC of the river model was 0.9 and that of the road model was 0.6, then w1 is 0.6 and w2 is 0.4.
Each probability was produced in each direction: N in formula (2) indicates the normalized value to avoid the reduction of probability caused by products of probabilities. We determined N using the sum of probability values in the single model that showed the best performance among the combined variables. In the case of the combination of river density and road density, if the river model showed the better performance, N = sum (e_river1, e_river2, e_river3, and e_river4) (6), whereas if the road model showed better performance, N = sum (e_road1, e_road 2, e_road 3, and e_road 4) (7).
Therefore, the sum of expansion probability from one cell in the combined model was conserved.
In combined model 2, we conducted a simple combination method using addition among results of independent simulations with a weighted function that decided the AUC value, which was determined by equations (3) and (4). We used this combined method because there was a possibility that A. bungii could expand using proxies such as rivers or roads. For example, the combination of river density and road density could be described as:

Result of combined model 2 = w1(result of river) + w2(result of road) (8),
where the river results indicate the virtual invasion number in the river model, and road results indicate the virtual invasion number in the road model. The functions w1 and w2 were the same for equations (3) and (4).
The simulation models were implemented using C language and compiled by gcc (GCC).

Validation using the current distribution range
We validated the results of the CA model using the merged occurrence records of A. bungii that were available up to 2019. These records contain data collected from 2017 to 2019, reflecting the current distribution record. First, we compared the average virtual invasion values between grid cells that did and did not contain occurrence records. We also compared the average virtual invasion values between grid cells containing occurrence records and all cells in the study area. If the model prediction accurately represents A. bungii dispersion in the real world, the invasion values should be higher in cells containing occurrence records than in those without and should be higher than the average for all analyzed grid cells. We used the Mann-Whitney U test to compare model predictions with field data. For the combination model, we combined CA fields that accurately showed a relatively high prediction in this procedure.
We subsequently evaluated the performance of the model using a receiver operating characteristic curve (ROC) analysis with area under the curve (AUC), sensitivity, and selectivity values (Lobo and Jiménez-Valverde 2008). AUC values use presence/absence data and range from 0 (for an inverse model) to 0.5 (for a random model) and finally 1 (for a perfect model). To calculate the AUC value, we designated grid cells containing occurrence records as "presence" and cells without occurrence records as "absence." Sensitivity was the ratio of correctly predicted "presence" events, whereas selectivity was the ratio of correctly predicted "absence." We decided that the cutoff point was optimal when the sensitivity + selectivity value was maximal. As a predictor value, we used the number of invasions in the cell divided by the maximum number of invasions in the target area.
All statistical analyses were performed using R ver. 3.6.1 (R Development Core Team).

Application of the model to a wider range
We applied a virtual ecology model that showed high performance based on the evaluation results for the prediction of a wide range. We used the two models that showed the best performance among the combined models

Distribution records of A. bungii
We acquired and merged the distribution records of A. bungii in the target area from 2017 to 2019. A total of 128 1 km cells within the 1,140 cells were designated as occurrences in the target area (Figure 4a).

Accuracy of the CA model simulation
The results of the model accuracy analysis using the Mann-Whitney U test are shown in Table 2. All three independent models showed a higher number of virtual invasions in cells with occurrence records than in those without occurrence and cells both with and without occurrences (i.e., all cells), whereas the null model showed reverse results; namely, a lower number of virtual invasions in cells with occurrence records than cells both with and without occurrences (  Table 2). Thus, our virtual ecology models using river density showed the best performance among the three independent models. The second-best performance model was the road model, in which cells with occurrence records (403.23 ± 503.64) exhibited a higher number of invasions than those without occurrence (323.38 ± 473.59) and all cells (332.34 ± 477.499), even when these were not statistically significant (Mann-Whitney U test, P = 0.05 and 0.08, respectively). The models based on the number of cherry trees had relatively low performance because cells with occurrence records (0.047 ± 0.21) did not exhibit a significantly higher number of invasions than those without occurrence (0.038 ± 0.21) and all-analyzed cells (0.039 ± 0.21) (Mann-Whitney U test, P = 0.45 and 0.50, respectively; Table 2). Accordingly, we conducted combined simulations using both the river and road density. Both river and road combined models 1 and 2 showed a higher number of virtual invasions in cells with occurrence records than in those without occurrence and all cells (Table 2). In combined model 1, the cells with occurrence records (587.94 ± 367.58; mean ± S.D.) exhibited a significantly higher number of invasions than those without occurrence (310.16 ± 393.469) and all cells (341.35 ± 400.23), respectively (Mann-Whitney U test, P > 0.001;  Table 2). Additionally, in combined model 2, the cells with occurrence records (522.78 ± 328.09; mean ± S.D.) exhibited a significantly higher number of invasions than those without occurrence (292.76 ± 324.49) and all cells (318.59 ± 332.78), respectively (Mann-Whitney U test, P > 0.001; Table 2). The results of the ROC analysis using all six models are shown in Table 3. Among these, the river and road combined model 1 had the best performance, with an AUC value of 0.68, and combined model 2 showed the secondhighest value of AUC (0.67). The river model showed the third-highest value of AUC (0.65). The AUC values of the other three models were all < 0.6 ( Table 3).

Prediction of invasion probability in Saitama Prefecture
Using the two developed models (river and road combined model 1, which showed the best performance among all models, and the river model, which showed the best performance among single variable models), we performed an expansion simulation using 100 steps for the entire area of Saitama Prefecture (which has 7,700 1 km 2 cells in total; Figure 5). The simulations indicated relatively low invasion probabilities in both models for the western areas of the Prefecture, whereas the mid to eastern areas showed relatively high invasion probabilities ( Figure 5). In the combined model, there were linear low-probability zones in the mid-eastern areas (Figure 5a). On the other hand, the river model showed an especially high probability in the western-north area (Figure 5b). We conducted 20, 50, and 100 steps on that procedure, with almost identical result trends, i.e., relatively low in the western region to relatively high in the eastern region. However, we could not perform 999 steps, which is the number of steps used in the model evaluation process, because of memory limitations of the PCs used. Thus, we conducted simulations using these model for the tested region in Saitama Prefecture with both 100 and 999 steps; we obtained basically the same trends in 100 steps (cells with occurrence records [509.75 ± 278.27] exhibited a significantly higher number of invasions than those without occurrence [220.94 ± 294.39] and all cells [253.369 ± 306.40]; Mann-Whitney U test, P > 0.001 for both) in the combined model 1, and cells with occurrence records (362.75 ± 288.89) exhibited a significantly higher number of invasions than those without occurrence (155.38 ± 246.77) Figure 5. Map of predicted occurrence units for the whole of Saitama Prefecture using both the river density model and river single model. The degree of shading reflects the theoretical invasion number predicted by each simulation model. and cells both with and without occurrences (178.66 ± 260.09; Mann-Whitney U test, P > 0.001 for both) for the river model, respectively.

Discussion
In the present study, we developed a dispersal simulation model of A. bungii based on CA with matrix resistance theory to establish a spatially explicit management plan. The simulation prediction based on combined river and road density per 1 km 2 cell as an expansion pathway exhibited prediction accuracy; virtual invasion values were higher in cells that contained occurrence records than those that did not and showed a higher AUC value among all simulation models for predicting the current distribution range of the species. Additionally, the simulation prediction based on river density only exhibited a prediction accuracy near the level of the combined model even when a single variable was used. Therefore, with these virtual ecology models, we were able to generate risk maps that can form the basis of a spatially explicit management plan for A. bungii.
Among the three candidate habitat/pathway proxies for our independent virtual ecology model, the river model had the highest accuracy, i.e., it reflected the real distribution range. In Japan, rows of cherry trees are often present along rivers (Kanebishi 2010) and roads (Kyoya et al. 2003). Thus, our findings suggest that A. bungii expands along those cherry tree rows as a dispersal pathway, i.e., continuous habitat. Although not significant, the road model could also predict the invasion of this species to some degree. The relatively low accuracy of the road model compared with the river model might be explained by variations in the size and type of roadway. There are several types of roads (such as highways and small paved roads) that have different widths and the presence/absence of a footpath. The decision on whether to plant rows of cherry trees is affected by such differences. In addition, among roads, there are differences regarding the density of cherry trees in the tree line (Kyoya et al. 2003). In this study, we used only road lengths without considering their type and size. The inclusion of road type and size in the model of information, which would reveal the tendency to have cherry tree rows, would probably improve its prediction accuracy.
By contrast, the model using the number of cherry trees did not predict the actual distribution range of A. bungii. One possible explanation for this result is lack of data. Our distribution data on cherry trees were mainly reported from public facilities, such as parks and schools, managed by municipal governments. Thus, our data did not include continuous rows of trees, including riversides and roadsides. Moreover, the data did not include private lands, such as gardens. Therefore, the simulation results might be inconsistent with actual distributions. In this context, the use of the proxy of cherry tree rows, rather than the number of trees, is reasonable for a virtual ecology approach.
Results of simulations using combined models of both river and road density, which exhibited relatively high performance in the independent model, showed better performance compared with single models. In our model, both river and road density were a proxy of cherry tree rows; thus, there was a possibility that A. bungii could expand using either option. Therefore, simulation using combined proxy models showed that better performance is basically reasonable. However, interestingly, model 1, which used multiplication of proxies, showed better performance than model 2, which used addition of results, and even the difference in AUC value was small. In model 1, if either value was low, invasive success could not occur even if the other value was high. This result suggests that river density is more essential than road density for expansion of A. bungii. The performance of the river single model, which showed the best performance of the three single variable models, also supports that. One possible explanation for this is tree density other than cherry trees. In general, there are often tree lines along rivers and roads, but the density of trees is high along rivers compared to roads in Japan. Trees other than cherry could not be used as feeding habit but could serve as dispersal corridors. A. bungii might use tree lines along rivers as expansion corridors even if the tree species are not food resources. From this perspective, we consider river density the most useful variable for prediction of A. bungii in the current situation, i.e., at a very early stage of invasion and expansion in its range.
We used the distribution records of A. bungii in the target area from 2017 to 2019 to evaluate model performance. However, the distribution range of this species is expanding (CESS, https://www.pref.saitama.lg.jp/cess/ center/kubiaka.html, accessed 25 August 2021); as such, it did not reach an equilibrium stage. Moreover, our method of data collection was not exhaustive; thus, our distribution records could include underestimation. This could influence the performance of the model and is one of the essential challenges in the prediction of the potential ranges of invasive species (Menke et al. 2009;Elith et al. 2010;Václavík and Meentemeyer 2012). Nevertheless, the combined model and river model, which used these imperfect data, yielded an AUC value greater than 0.65. This suggests that an area with a high number of theoretical invasions in the model could reflect invasion success at a relatively early stage. Therefore, for management practice, areas without detected empirical invasions that possess high numbers of simulated invasions in our model could be prioritized in surveillance efforts.
The results of the expansion simulation for the whole area of Saitama Prefecture indicate that the general trend of invasive risk is gradually higher from West to East, even though there were some anomalies ( Figure 5). These results are plausible because the western area comprises mainly mountainous and hilly areas that are covered with various deciduous and coniferous trees and, thus, a small number of planted cherry trees exist along rivers and roadsides. In contrast, the eastern area is urbanized and has a relatively large number of planted cherry trees along the rivers and roadsides, as well as in public spaces (i.e., parks and schools). Therefore, the probability of A. bungii invasions might be higher in the eastern than in the western areas of this Prefecture. The linear regions in the eastern portion of the Prefecture that showed a consistently low probability from the result of the combined model 1 are basically comprised of agricultural zones, which have low road density; thus, invasion probability could decrease. In combined model 1, if neither the river nor road value was low, invasive success could not occur even if the other value was high. Although its AUC value was the highest among all models, this could cause underestimation when we applied that model for unknown areas. This is one application consideration of our model. To establish a management plan, decision-makers should adopt conservative, precautionary methods; thus, we should employ the most severe result when simulation results between models do not agree. In addition, the river single model is simple compared with the combined model. Therefore, we would recommend the river single model for practical use.
Some locations in the Saitama Prefecture rely on mature cherry trees for tourism income, particularly in its eastern area, such as Gongendo-tsutsumi in Satte city (Gongendo-tsutsumi), where a dead male adult was collected in our survey in 2020 (CESS, https://cessgis.maps.arcgis.com/apps/ webappviewer/index.html?id=67495ef407974aa7bdcbfcb42749471b, accessed 1 July 2020). The economic damage to such tourism spots might be severe if they are invaded by A. bungii. Based on our prediction map, land managers and/or practitioners should focus on prevention and eradication efforts in such high-risk areas and/or the rivers and roads that connect with such areas. More specifically, we recommend allocating eradication efforts in the mid-northeast region of the Saitama Prefecture, where both simulation results showed high invasion probability. This area is a border zone in Gunma Prefecture close to Tatebayashi city, Itakura town, Chiyoda town, Meiwa town, Oizumni town, and Oura town, which has already been severely damaged by A. bungii and a council for countermeasures established (Tatebayashi city). For effective management, efforts should be concentrated there, and collaboration with these local governments is an important step.

Conclusion
In the present study, we used a virtual ecology approach to establish an expansion prediction model of A. bungii that could predict areas with a high invasion risk. One of the most cost effective approaches to the management of invasive species is the implementation of control measures at the early stage of invasion (Simberloff 2003;Kaplan et al. 2014). Predicting future species expansion is an important step toward achieving this goal (Ficetola et al. 2007;Fukasawa et al. 2009;Koike and Iwasaki 2011;Osawa and Ito 2015). Thus, our prediction model supports the establishment of cost-effective management plans. Moreover, our approach actualizes the prediction of high-risk areas without the use of detailed information on the target species in the target area. Based on the results of this approach, decision-makers can pre-emptively establish effective management plans by identifying areas in which to focus their efforts at the early stages of invasion.

Funding declaration
This study was partially supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Number 16H05061 (TO).