Using a

Landscape Ecol (2010) 25:353–369 DOI 10.1007/s10980-009-9434-9 RESEARCH ARTICLE Modeling the invasive emerald ash borer risk of spread using a spatially explicit cellular model Anantha M. Prasad • Louis R. Iverson • Matthew P. Peters Jonathan M. Bossenbroek • Stephen N. Matthews • T. Davis Sydnor • Mark W. Schwartz • Received: 30 June 2009 / Accepted: 23 November 2009 / Published online: 11 December 2009 Ó US Government 2009 Abstract The emerald ash borer (EAB, Agrilus planipennis) is decimating native ashes (Fraxinus sp.) throughout midwestern North America, killing millions of trees over the years. With plenty of ash available throughout the continent, the spread of this destructive insect is likely to continue. We estimate that the insect has been moving along a ‘‘front’’ at about 20 km/year since about 1998, but more alarming is its long-range dispersal into new locations facilitated by human activities. We describe a spatially explicit cell-based model used to calculate Electronic supplementary material The online version of this article (doi:10.1007/s10980-009-9434-9) contains supplementary material, which is available to authorized users. A. M. Prasad (&)
L. R. Iverson
M. P. Peters
S. N. Matthews USDA Forest Service, Northern Research Station, 359 Main Road, Delaware, OH 43015, USA e-mail: [email protected] L. R. Iverson e-mail: [email protected] M. P. Peters e-mail: [email protected] S. N. Matthews e-mail: [email protected] risk of spread in Ohio, by combining the insect’s flight and short-range dispersal (‘‘insect flight’’) with human-facilitated, long-range dispersal (‘‘insect ride’’). This hybrid model requires estimates of EAB abundance, ash abundance, major roads and traffic density, campground size and usage, distance from the core infested zone, wood products industry size and type of wood usage, and human population density. With the ‘‘insect flight’’ model, probability of movement is dependent on EAB abundance in the source cells, the quantity of ash in the target cells, and the distances between them. With the ‘‘insect-ride’’ model, we modify the value related to ash abundance based on factors related to potential human-assisted movements of EAB-infested ash wood or just T. Davis Sydnor School of Environment and Natural Resources, Ohio State University, 367 Kottman Hall, 2021 Coffey Road, Columbus, OH 43210, USA e-mail: [email protected] M. W. Schwartz Department of Environmental Science and Policy, University of California-Davis, One Shields Avenue, Davis, CA 95616, USA e-mail: [email protected] J. M. Bossenbroek Department of Environmental Sciences & Lake Erie Center, University of Toledo, 6200 Bayshore Rd., Oregon, OH 43618, USA e-mail: [email protected] 123 354 hitchhiking insects. We attempt to show the advantage of our model compared to statistical approaches and to justify its practical value to field managers working with imperfect knowledge. We stress the importance of the road network in distributing insects to new geographically dispersed sites in Ohio, where 84% were within 1 km of a major highway. Keywords Emerald ash borer
EAB
Agrilus planipennis
Spread model
Stratified dispersal
Spatially explicit cellular model Ohio
Gravity model
Fraxinus
Ash
Roads networks
Invasive
Highway traffic
Insect flight model
Insect ride model Introduction The emerald ash borer (EAB), Agrilus planipennis Fairmaire (Coleoptera: Buprestidae), poses a serious threat to all native ash trees in North America, especially in the eastern United States. The larvae feed on phloem, producing galleries that kill large trees in 3–4 years and small trees in as little as 1 year (Poland and McCullough 2006; Wei et al. 2007). A native of northeastern China, Korea, Japan, Mongolia, Taiwan, and eastern Russia, the species was first identified in the United States near Detroit, MI, in July 2002 (Haack 2006; Poland and McCullough 2006). The borer was thought to be imported into Michigan in the early 1990s via infested ash crating or pallets (Herms et al. 2004). Since its initial establishment in the early to mid 1990s (Siegert et al. 2007), it has spread at an accelerating pace from that position. So far, all attempts to stop the spread of the organism have failed. In the initial years after the EAB was first identified, many eradication efforts in Michigan, Ohio, Maryland, and Ontario were attempted. Typically, all ash trees within an 800-m radius of the initial detection tree(s) were being cut, chipped to very small pieces, and incinerated. This expensive program was halted in Ohio by early 2006 because of funding shortages and because of numerous newly discovered infestations. The primary hope to slow the spread now lies with the introduction of 123 Landscape Ecol (2010) 25:353–369 specialized natural enemies (Gould 2007; Bauer et al. 2007), monitoring and education programs, highly targeted ash removals, and regulation. The impact of EAB may be enormous. An estimated 8 billion ash trees exist in the United States, comprising roughly 7.5% of the volume of hardwood sawtimber, 14% of the urban leaf area (as estimated across eight US cities), and valued at more than $300 billion (Poland and McCullough 2006; Sydnor et al. 2007). Spread models The population dynamics during the establishment of any pest is influenced both by the Allee effect (population size is highly correlated with survival of individuals) and demographic stochasticity (Liebhold and Tobin 2008). If the population overcomes the Allee effect, spread occurs with the growth of population and subsequent dispersal. The dispersal mechanism has been studied extensively by various researchers and there is considerable literature in this field. We will not attempt a detailed review of the spread models here, but will provide a brief overview so that we can characterize our hybrid modeling approach. The models can be broadly divided into mathematical and cellular approaches. Mathematical models One of the first attempts to tackle spread characterized the initial dispersal by continuous random movement defined by reaction-diffusion partial differential equations (Skellam 1951). To keep the interface simple, we prefer not to reproduce the equations involved here—the interested readers can refer to the sources for mathematical exposition. The magnitude of the dispersal depends on the diffusion coefficient (the standard deviation of dispersal distance), often derived from empirical data. The simple diffusion models are deterministic—hence, the front forms a concentric ring of expanding range at constant speed. Although simple diffusion models can accommodate interspecific interactions like competition, predation etc., they cannot account for the long-range dispersal of the insect which are mostly facilitated by weather and human agents. These longrange movements are often characterized by fat-tailed Landscape Ecol (2010) 25:353–369 distributions (e.g., a higher probability than a normally distributed variable of extreme values). The combination of simple diffusion and long-range dispersal mechanisms is termed ‘stratified dispersal’ (Hengeveld 1989; Shigesada et al. 1995). The stratified dispersal mechanism can be defined by integro-difference equations which are timediscrete and space-continuous models which consist of two parts: population growth (defined through difference equations) and dispersal in space (defined with an integral operator; Kot et al. 1996). Like reaction-diffusion equations, they can generate constant-speed traveling waves. In addition, they can account for invasions whose spread rates increase with time (longer range dispersal). The dispersal kernel, which is a probability density function, determines the speed and movement of individuals between two points in space (Lewis 1997). While the shape of the dispersal kernel is an important factor in determining spread, by itself it tends to overestimate the dispersion speed if demography is not considered (Neubert and Caswell 2000). Efforts have been underway to incorporate stochastic variation in dispersal and reproduction into density-dependent models (Clark et al. 2001; Snyder 2003) and to differentiate effects of stochasticity from the effects of nonlinearity (Kot et al. 2004). However, it is likely that inherent uncertainty rather than parameter sensitivity makes long-distance dispersal difficult to predict (Clark et al. 2003). Cellular models A major limitation of using reaction-diffusion and integro-difference models is that they do not consider how dispersal and demography are influenced by spatial pattern within the landscape, as they treat the landscape as homogeneous for mathematical tractability (With 2002). Several theoretical explorations indicate that spread rates are affected by habitat fragmentation and other aspects of the spatial arrangement of favorable habitat (Lonsdale 1999; With 2002). Therefore, it is desirable to incorporate spatial structure into the dispersal models. This is where the spatially explicit cellular approach based on cellular automata is appropriate. It involves fitting transition state models to real systems where the transition between one time step and the next depends on the empirical relationship between the current 355 state of the target cell and that of its neighbors (Molofsky and Bever 2004). There have been previous attempts to model the spread of EAB with spatial dynamic models using STELLA and Spatial Modeling Environment (SME) in conjunction with high resolution landscape data (BenDor and Metcalf 2006; BenDor et al. 2006). However, these studies were limited to one county in Illinois and based on limited data available at the time. Our approach uses a simple fat-tailed power function to model the spread (Schwartz 1992) by taking into account stochasticity, demographics, and anthropogenic factors that influence the spread of insects (including the results of a ‘gravity model’ where movements are not random but biased by the attractiveness of destinations) in a spatially explicit cellular framework. We also take fragmentation of habitat into account because the spread depends on the ash basal area in each cell which reflects the fragmented nature of ash and EAB distribution across current landscapes. We therefore address the limitation of purely mathematical approaches by taking landscape heterogeneity into account. Specifically, we incorporate both high-resolution ash quantity maps and human-assisted components to build a model of EAB risk for Ohio. We modified a previously developed spatially explicit cellular spread model called SHIFT (Schwartz 1992; Iverson et al. 1999; Schwartz et al. 2001) with an integrated approach that combines the insect’s short-distance movement patterns, which we call the insect flight model (IFM), with a number of avenues for human-assisted transmittal including road networks, campgrounds, wood product industries, and human population density, which we call the insect ride model (IRM). It should be noted that because we combine human facilitated risk factors in our hybrid approach, we are modeling the ‘risk of spread’ (Drake and Lodge 2006). Overall approach and SHIFT modeling scheme The modified SHIFT model calculates the probability of EAB infestation of currently unoccupied cells (270 9 270 m) based on the abundance of EAB in the occupied cells (as assessed by the years since infestation), the habitat availability of ash (as assessed by the ash basal area), and the distance 123 356 Landscape Ecol (2010) 25:353–369 between all occupied and unoccupied cells within a search window. The model uses the inverse power law (Gregory 1968) to calculate the probability of an unoccupied cell becoming infested during each generation. While this relation describes the asymptotic distribution of seeds, spores, or pollen (Okubo and Levin 1989), we have modified it to suit our objective of tracking the risk of EAB spread. The relation is: Pi;t ¼ HQi n  X j¼1   HQj  Fj;t  C=DXi;j ! ð1Þ where Pi,t is the probability of unoccupied cell i being infested to a detectable level at time t; HQi and HQj are habitat quality (i.e., ash abundance) scalars for unoccupied cell i and occupied cell j, respectively, that are based on the total basal area of ash (m2/ha) in each 270-m cell; Fj,t, an abundance scalar (0–1), is the current estimated abundance of EAB in the occupied cell j based on the years since infestation at time t; Di,j is the distance between unoccupied cell i and an occupied cell j; and n is the number of cells in the search window. The value of C, a rate constant, is derived independently through trial runs to achieve a dispersal rate of approximately 20 km per year, a value estimated from empirical data of EAB finds (Iverson et al. 2008a). The value of X, or dispersal exponent, determines the rate at which dispersal declines with distance. Being in the denominator, this decreases infestation risk with distance as an inverse power function. For this model, we make the simplifying assumption that a generation for the insect can be accomplished in 1 year, although we realize that the cycle can last 2 years, especially in healthy, newly invaded stands or those that are located in far north zones with less than about 150 days of growing season length (Poland and McCullough 2006; Siegert et al. 2007; Wei et al. 2007). The infestation probability for each unoccupied cell, a value between 0 and 1, is summed across all occupied cells at each generation. Thus, an unoccupied cell very close to numerous occupied cells may end up with an infestation probability [1.0. These cells are modeled with a 100% probability of being infested. For cells with summed infestation probabilities \1, a random number \1.0 is chosen and all cells with a probability of infestation that exceeds the 123 random number are infested in that model step. This adds an element of stochasticity to the otherwise deterministic model. Those ‘‘newly infested’’ cells then contribute to the infestation probability of unoccupied cells in later model time steps. Further discussion of the dispersal function can be found in Schwartz et al. (2001). The SHIFT model advances the ‘‘front’’ based on the current front location, the abundance of EAB behind the front, and the quantity of ash ahead of the front. Based on reports from the early infested zones in the Detroit area (e.g., Siegert et al. 2007) and early data, we assumed a 10-year span from when EAB is initially detected to the death of all ash trees within the 270 m 9 270 m cell. Recent field evidence indicates that the beetle can kill individual trees within 3–6 years of detection (Kathleen Knight, US Forest Service, personal communication; Peters et al. in press), although it takes longer to kill all the trees in a cell. Because we do not start the year counter until the cell has detectable EAB (visual symptoms have occurred, or EAB densities are high enough to be trapped), the model is insensitive to the early years of population increase following initial colonization. EAB abundance in the cell was assumed to form a modified (skewed left) bell-shaped curve, with the multipliers picked from the curve. The maximum abundance (multiplier = 1) occurred in years 3, 4, and 5, a 0.61 multiplier in years 2 and 6, a 0.14 multiplier in years 1 and 7, and a 0.01 multiplier in years 8, 9, and 10. The assumptions for this curve include a slow EAB population increase for the first few years after initial infestation (but not yet detected or counted in this 10-year cycle), followed by rapid population growth once detected, followed by peak infestation for 3 years starting with year 3, followed by a rapid decline as all the ash trees in the cell die off in years 7–10. We realize that small ash saplings remain viable and new ash germinate and establish within the cell subsequent to the major die off, but our model ignores these in subsequent generations because of the difficulty in keeping track of seedling and EAB demographics. Research underway now will determine if regenerating ash can sustain a viable EAB population in the wake of the killing first wave. For each cell, the insect flight model (IFM) calculates the probability of a new infestation, based on our empirically derived spread rates (*20 km/ year, Fig. 1) which include a small probability that Landscape Ecol (2010) 25:353–369 357 Fig. 1 Estimate of EAB spread 1998–2006, defined by positive detection trees emanating out from the original infestation the insect will fly or be transported from an occupied cell to an unoccupied cell, for all surrounding cells within a specified search window of 40 km. This distance includes the EAB’s historical and potential flight patterns (Taylor et al. 2005, 2007), but also some movement other than flight (e.g., wind and other short-distance transport mechanisms). For the insect ride model (IRM), we use the same model modified for long-distance dispersal (see later section) up to 400 km, based on observed patterns of long-distance dispersal and size of the state of Ohio. Once selected for infestation, the cell starts the 10-year cycle of EAB increasing and then decreasing as ash dies out. In our modeling strategy, we separately model the broad insect flight movement with the IFM and the human-facilitated long-distance spread of risk through various agents like roads, campgrounds, wood products industries, and human population density with the IRM. The spatial addition of these two models captures the relative risk of infestation of any 270-m cell within Ohio by EAB. Since we are modeling risk of spread by combining two models, it is difficult to assign a time element to the prediction; however, based on our assessment of 2005–2007 EAB detections, a reasonable time frame to evaluate the utility of our map would be 2–4 years. To develop the IRM, we used GIS data to weight factors related to potential human-assisted movements 123 358 of EAB using various agents/factors that assist in its rapid movement: traffic on roads, various wood products industries, population density, and campgrounds. Each of these four factors was converted into weighting layers that became multipliers for the ash basal area component of the IRM: in effect increasing the probability of EAB infestation according to the insect ride factors by boosting the amount of ash available in those cells. Thus, if no ash exists in a particular cell, it does not matter whether there is an escaped EAB from one of the human-assisted factors, but if there is an ash component, an escaped EAB would find the cell ‘more’ attractive to colonize in proportion to the basal area of ash. We describe the methods and the weighting schemes used in the following sections. It should be kept in mind that with our ‘spatially explicit’ cellular model, we are inevitably balancing the tradeoffs between generality, precision, and realism that limit all model development in ecology (Levins 1966). In our approach, generality is somewhat sacrificed because environmentally dependent rates of spread make it difficult for models parameterized in one location to make accurate predictions in another (Hastings et al. 2005). Therefore, the parameters and weights presented here are specific to Ohio, and, though the same approach should be appropriate for other regions, a re-evaluation of the model parameters and weights would likely be necessary. Methods Ash quantity estimation We have compiled estimates of ash supply at two scales. To visualize the amount of ash available over all the eastern United States, we used forest inventory and analysis (FIA) data (Miles et al. 2001) compiled at a 20 9 20 km scale (see Prasad and Iverson 2003; Prasad et al. 2007; Iverson et al. 2008b; Supplementary Figs. 1, 2, 3, 4). For model input, we created a detailed estimate of ash basal area for Ohio (Iverson et al. 2008a). This estimate, at a scale of 270 9 270 m, was created by combining estimates of ash basal area per FIA plot with a Landsat TM-based classification of forest types, resulting in the total ash basal area in m2/ha (Supplementary Fig. 5). 123 Landscape Ecol (2010) 25:353–369 Collection of EAB positives Between 2005 and 2007, the Ohio Department of Agriculture set up nearly 20,000 detection trees throughout the state to monitor for EAB. Detection trees were then felled and peeled following adult emergence (after August) and assayed for EAB. Both the detection tree locations and the positive finds were provided to us for analysis. Core area of EAB infestation To model risk of potential future spread as well as establish a core area from which the front is spreading, we estimated historical spread from 1998 to 2006 (Fig. 1). For this estimate of spread that included years prior to first identification of the insect (1998–2002), multiple data sets, GIS processes, and assumptions were used, and for the most part, represent the spread of visible damage to ash trees rather than the initial infestation of EAB. First, we set the 2002 and 2003 EAB boundaries by using webpublished pest maps (accessed from www.michigan. gov/dnr in 2005, no longer online) of the Michigan Department of Natural Resources. The primary source for pre-2002 estimates were web-published maps by Smitley and others (internet maps no longer on the web but published in Smitley et al. 2008). These maps were based on interstate highway exit surveys for 2003 and 2004, with ten ash trees per exit tallied for death/life (irrespective of what killed the ash). The authors then mapped percent ash dead in classes of 0–10, 11–40, 41–80, and 81–100%. We assumed it took 5 years for a newly infested tree to die, so that the location of the class 80–100% dead in 2003 was assumed to be infested in 1998. Similarly, the additional area of the class 80–100% dead in 2004 was assumed infested in 1999, and so on. These data allowed us to manually delineate estimated infestation zones for 1998–2002. For boundary estimates after 2002, several other sources were used. A Michigan ash damage survey from September 2004, along with the actual locations and density of known EAB locations as of December 2005, were obtained from the Cooperative Emerald Ash Borer Project. In addition, multiple dates of national EAB positive maps were acquired to detect additional finds temporally (http://emeraldashborer.info/surveyinfo.cfm). Finally, our own field work on ash tree assessment in Landscape Ecol (2010) 25:353–369 northern Ohio and southern Michigan during the summers of 2004 and 2005 yielded additional spatial information, particularly on ash not yet visually affected by EAB. Further details with maps are found in Iverson et al. (2008a). We define the ‘front’ by the density of newly identified positive detection trees each year. Recent studies of tree rings in the initial zone of infestation have indicated that initial death of ash trees occurred in 1997 (Siegert et al. 2007). Although EAB was certainly present before 1998, we assume that EAB abundance was low and not easily detected before that time. Model inputs Weighting the parameters for the insect ride model was difficult because we did not have reliable, independent means of arriving at weights to the inputs of the model. Iteratively arriving at a fit was prohibitive in terms of time because the SHIFT model is computationally intensive. We selected an alternative, decision tree-based approach called ‘‘RandomForest’’ (see Prasad et al. 2006) to test our hypothesis that the ‘non-core’ locations of geographically dispersed EAB infested locations (which we henceforth term ‘outlier’) were correlated with road networks, campgrounds, wood industries, and distance from the core-infested zone, and to assign weights to these factors in the IRM based on the importance of the predictors. RandomForest, which is based on averaging numerous decision trees by perturbing both data and predictors have been shown to be a superior modeling technique for many applications where imputed predictions are desired (Hudak et al. 2008), but may not be an appropriate modeling technique for applications with a strong spatial dependency (e.g., EAB dispersing from a small section of the territory). However, RandomForest can provide a list of variable importance scores for the predictors that can be used as a guide to assign weights to the IRM. In order to determine the importance of the predictors in the RandomForest (RF) model, we generated the average distance to roads, wood product industries, campgrounds, and occupied core zone from both the EAB positive trees and the non-positive detection trees for the outlier zone, and then randomly sampled by county from the outlier positives (from a total of 255) and the 359 non-positive detection trees (from a total of 9,964) in a 1 positive to 3 non-positive ratio. We stratified by county to reduce the over-weighting of some counties that had a disproportionate number of positives. The most important predictor, expectedly, was distance to occupied core zone followed by average distance to major roads, wood products industries, and campgrounds. We also derived a imputed map (with an error estimate of 8%) for the RF model to compare with our modified SHIFT model (Supplementary Fig. 6). We discuss this in the summarizing model outputs section. Our weighting scheme consisted of using the relative values of RF-derived importance scores, along with expert opinion, to assign weights which essentially ‘‘boost’’ the ash basal area estimates and hence ‘‘attractiveness’’ of the cell to EAB. While we realize that the weightings described next may not be appropriate across regions, extensive details for Ohio are provided to facilitate application for other regions (Fig. 2). Roads Intuitively from studying maps of the reported positives, and also through our predictor importance scores derived from the RandomForest model, it was clear that the network of major highways was highly coincident to the locations of new outlier infestations (see Fig. 3). Evidence is mounting that many of the new introductions may be due to insect hitchhiking, in addition to the transport of materials like firewood (Buck and Marshall in press). For example, a recent detection (June, 2009) in Cattaraugus County, New York, is very close to a busy exit on Route I-86. To register the increased probability of insects riding on windshields, radiators, or otherwise attached to vehicles moving down the road, we assigned weights to two widths (1 and 2 km) of major road corridors. We used the National Highway Planning Network data for the first half of 2007 (Federal Highway Administration 2007), which reported average daily traffic (ADT) values calculated as the number of vehicle miles traveled divided by the number of center line miles of highway. Classes of ADT values, reported at the county level for all but the smallest highways, were assigned weights to the inside buffer (1 km), ranging from 1 to 60 as a linear proportion to their traffic density, which ranged from 2,000 to 123 360 Landscape Ecol (2010) 25:353–369 Fig. 2 Inputs to the model of EAB risk (see text for weighting schemes): a major roads rated by traffic density in Ohio; b campground weights; c wood product weights; and d human population density by zip code 164,000 cars per day. Naturally, traffic was generally much lower in rural zones, with ADT on rural interstate highways of \50,000 (Fig. 2a). Weights were reduced by half for the outside buffer of 1–2-km distance from each road. Campgrounds Campgrounds have often been considered likely destinations of human-assisted EAB transport, primarily through the movement of firewood. There are many documented or highly suspected cases for this mode of movement in Michigan and more recently into West Virginia (West Virginia Division of Forestry 2007). Because the general public is involved, it is much more difficult (compared to wood industries and nurseries) to achieve education, regulation, compliance, and enforcement goals 123 related to stopping EAB spread. Campground locations were acquired from Dunn & Bradstreet and the AAA Travel and Insurance Company. We weighted the campgrounds with a ‘gravity model’ (described next) combined with a campground coverage buffered and scored according to number of campground sites. Gravity model Unlike typical dispersal and spread models, gravity models explicitly assume that movements are not random but biased by the attractiveness of destinations. They allow for the prediction of long-distance dispersal events by considering not only the nature of source populations, but also the spatial configuration and nature of potential colonization sites. Because of this, gravity models have the potential to more Landscape Ecol (2010) 25:353–369 361 Fig. 3 Locations of a detection trees in a 2006 and b 2007, and detected EAB positives in c 2006 and d 2007 accurately forecast species movement through heterogeneous landscapes than do diffusion models, which do not explicitly consider the spatial pattern of distant sites (Bossenbroek et al. 2001). Our gravity model considers traffic volumes and routes between EAB source areas and various distances to campgrounds. Muirhead et al. (2006) earlier used a gravity model approach to movement of EAB; our current model expands on that work by incorporating more detailed data on both public and private campgrounds. Because empirical data on the use of campgrounds, i.e., reservation data, are available only for public campgrounds, a modeling framework is necessary to incorporate private campgrounds; in this case to predict the relative number of campers traveling from EAB-infested areas to Ohio campgrounds. Gravity models calculate the number of individuals (e.g., campers) who travel from location i (i.e., zip code) to destination j (i.e., a campground), Tij, as estimated as Tij ¼ Ai Oi Wj ca ij ð2Þ where, Ai is a scalar for location i (see below), Oi is the number people at location i, Wj is the ‘‘attractiveness’’ (=number of camp sites) of location j, cij is the distance from location i to location j, and a is a distance coefficient, or distance-decay parameter, which defines how much of a deterrent distance is to interaction. Ai is estimated via XN W ca ð3Þ Ai ¼ 1= j¼1 j ij where N represents the total number of destinations and j represents each destination in the study region. We recorded a total of 241 public and private campgrounds in Ohio; in the absence of actual usage data, we assumed a direct relationship between the 123 362 number of camp sites and the usage of each campground (W). The distance between a zip code and a campground (c) was calculated as the road network distance between these locations; for simplification, we used only roads with either a state or federal designation, excluding local roads. The number of campers in each zip code (Ai) was assumed to be proportional to the total population, thus we used the population of each zip code based on 2000 census data. The point of origin within each zip code was the road location nearest the centroids of the zip code region, while the campground locations were determined as the nearest point to a state or federal road. The area of EAB infestation consisted of all the zip codes that had their centroids within the 2006 front. To estimate the distance coefficient (a), we compared our gravity model with reservation data obtained from the Ohio Division of Parks and Recreation for 58 state parks. These records contained the number of reservations for each campground summed by zip code of the camper’s residence. We used sum of squares to measure goodness-of-fit between model predictions and the observed data (Hilborn and Mangel 1997). By fitting the model to the reservation data for Ohio state parks, we assume that campers using private and public campgrounds behave in the same manner, i.e., the distance to and size of a particular campground affect their travel decisions in the same manner. Once the gravity model was parameterized, we used the estimated distance coefficient value to determine the expected number of campers that would travel to all 241 campgrounds within Ohio. To give a relative estimate of risk we reported the percentage of campers coming from EAB-infested zip codes traveling to each campground in Ohio. Each campground in Ohio had a relative score derived from the gravity model calculations ranging from 51 to 174,800. These were divided into quartiles with breakpoints of 1,640, 3,902, and 8,552, and used with increasing (but arbitrary) multipliers of 1.0, 1.33, 1.67, and 2.0, respectively. In this way, those campgrounds with many visitors arriving from the core zone were given double weight. The spatial influence assumed for each campground in this analysis was a radius of *4 km around the campground headquarters. 123 Landscape Ecol (2010) 25:353–369 Weighting by size of campgrounds The weighting (10 points for inner buffer, 5 points for outer buffer) and buffer size (ranging from 0.5 km to 4 km) for campgrounds were based on the number of camp sites listed in Dunn & Bradstreet data (www.dnb.com/us/). The smallest campgrounds (e.g., \50 camp sites) had buffers of only 0.5 and 1 km, while the largest campgrounds ([600 camp sites) had buffers of 2.5 and 5 km. These scores of 5 and 10 were then multiplied by the multipliers from the gravity model (1.0–2.0) to yield a maximum possible total of 20 points for campgrounds (Fig. 2b). Wood products Wood products industries, including nurseries, also have been responsible for some EAB movement (including an outbreak near an ash handle factory in northwestern Ohio and another via landscaping a new restaurant in central Ohio), so a scheme was developed to weight buffers around individual businesses dealing in wood products. Because the wood products industries are regulated, we expect that the risk is relatively low now and will become lower in the future, but pre-regulation infestations or accidental (or not) introductions still might occur. First, we assessed industries for their potential to use ash logs, based on the listing of standard industrial classification (SIC) codes from Dunn and Bradstreet. We scored each industry for likelihood of EAB getting to the site and emerging based on our estimate of the amount and status of ash used in the industry: 0 = none; 2 = small likelihood; 4 = somewhat likely; 6 = higher likelihood. For example, wood pallet industries scored a 6, while manufacturers of pressed logs of sawdust or woodchips scored a 2. Movement of material from nurseries historically has been a source for several infestations, but is a minor component in this model of future risk because of the close scrutiny on nursery stock movements. Next, buffer distances around the businesses, and corresponding weights, were created based on the number of employees (surrogate for volume of wood) working at the facility. For 1–10 employees, the buffer of 0–1 km scored 8 and the 1–2 km buffer scored 3; for 11–50 employees, the buffer of 0–1.5 km scored 9 and the 1.5–3 km buffer scored 4; and if the facility had more than 50 employees, the Landscape Ecol (2010) 25:353–369 0–2 km buffer scored 10 and the 2–4 km buffer scored 5. Any overlapping buffers were summed with a maximum weight of 10. Because of the lower probability of movement via this mechanism (resulting from these industries being highly regulated), the total value of the wood products component was only 10% in the composite model (Fig. 2c). Population density Although road densities and ADT values are good surrogates for population densities, we included population by zip code to capture areas where humans could contribute to the spread through alternative activities or even through vehicular movement on minor roads, which were not included in our road network. This factor creates an imputed score over the entire state and distinguishes rural from more urbanized areas. Data were acquired for 2000 from the US Census Bureau and were divided into six classes with scoring as follows: 1 = 1–100 people/km2; 2 = 101–200; 4 = 201–800; 6 = 801–2,000; 8 = 2,001 –4,000; 10 = 4,001–16,582 (Fig. 2d). 363 raw ash in their industrial process (e.g., forest nurseries and wood pallet industries); and (4) the maximum score for human population density (by zip code) was 10, which was achieved when the population density exceeded 4,000 people/km2 in that particular zip code. With 60% of the potential risk attributed to roads, we underscore that road networks are the quickest and most likely mode of dispersal at present and increasingly in the future as campgrounds and wood product industries become increasingly well-regulated, and humans are better indoctrinated on the negatives of moving wood material. The actual scores in our data ranged from 0 to 77 (100 being the maximum possible). These final weights were then used as multipliers on the existing basal area in each 270 9 270-m cell throughout Ohio. Because the range of ash basal area per cell was 0.016–5.11 (plus we added 1 to each value to keep values [1 for GIS and quantitative purposes), the maximum possible score for the final insect ride classification was 611 (actual data ranged from 1 to 440). These values were scaled back to 1–100 for relative scoring and mapping of ‘risk’. Development of risk map Quarantine zones The Ohio Department of Agriculture quarantines a county when EAB is confirmed from anywhere in it. Once the county is quarantined, ash is legally allowed to be moved within the county or to any other adjacent quarantined counties. Thus, the possibilities for EAB spread are enhanced within those counties. We effectively modeled the spread risk within quarantined counties to be twice that of non-quarantined counties by finally reducing the scores achieved from each of the other four factors by half in nonquarantined counties. Overall weighting scheme The final weighting scheme was as follows: (1) the maximum score possible for roads was 60, meaning the traffic density was as much as 164,000 cars per day; (2) the maximum score for campgrounds was 20, meaning the campgrounds with a gravity model score [8,553 and the campground had [600 campsites scored the highest possible; (3) the maximum score for wood products industries was 10 when the number of employees exceeded 50 and they used considerable The final risk map was derived by summing the scores for the insect ride and the insect flight components of the model. As such, the flight model affects only zones within about 40 km from the core area, beyond which the ride model is responsible for the additional risk. The combination logically creates maximum risk near the present core zone of infestation. Because of the hybrid approach and the combining of two models, it was difficult to assign a prediction time for the risk. We estimate that the risk map should be useful for about 2–4 years, although it will be useful to rerun, and potentially improve, the model periodically as new data on outbreaks and other information becomes available. Verifying and summarizing model outputs Next we verified and summarized the model outputs according to the confirmed EAB positive finds as of December 2007. The number of positives that were recorded and that lay outside the 2006 core zone from December 2003, to December 2007, number 255: 1 in 2003, 7 in 2004, 78 in 2005, 110 in 2006, and 59 in 2007. Although visual inspection of the positives 123 364 influenced the model building, the actual locations were not associated with the model. Assessment of positives versus detection trees During late 2005 through 2006, the Ohio Department of Agriculture designated a total of 9,670 ash trees as EAB detection trees, and another 9,964 trees were designated in 2007 (Fig. 3a, b). Detection trees were girdled to make them more attractive to EAB and were set with sticky traps during the growing season. In the following months, the trees were removed, peeled, and inspected for the presence of EAB. If EAB was confirmed, the tree was considered positive. These positive tree locations were used to analyze the accuracy of our model (Fig. 3c, d). Some people have proposed that the higher proportion of new infestations near major roads was merely due to greater sampling of ash trees near the highways, rather than EAB establishing more frequently near Fig. 4 Risk map for EAB in Ohio 123 Landscape Ecol (2010) 25:353–369 highways from transport of materials containing EAB and hitchhiking insects. To test this possibility, we analyzed the relationship between all positive trees and the 2006 and 2007 detection trees. By comparing the ratio of 2006/2007 positive trees to 2006/2007 detection trees for a series of variables (e.g., average basal areas of ash and average distances from roads, campgrounds, wood products, and the occupied zone), we could assess the relative importance of individual variables in detecting potential positives. Results and discussion Risk map of EAB spread The risk map visually shows extremely high risk near the current core zone, due to the additive risk of the flight and ride models, and the availability of nearby EAB (Fig. 4). Next in risk are the metropolitan Landscape Ecol (2010) 25:353–369 regions of Columbus, Lima, and Cincinnati/Dayton. These elevated risks are primarily due to all the humans moving around, some of which inevitably would carry EAB accidentally or via wood material. There is a wide variability of risk elsewhere, with areas near major roads showing more risk than minor roads or the rural areas. These relative risk levels will change as the EAB spreads, new counties are quarantined, and new invasions provide new sources of the insect to spread from. We provide three output maps of risk: (1) risk of the insect to fly into new zones (Supplementary Fig. 7); (2) risk of the insect to ride with humans (Supplementary Fig. 8); and (3) overall risk as a combination of the flight and ride models (Fig. 4). The risk maps show the large influence of several factors. First, it is readily apparent that distance from current EAB centers is very important—that regardless of how the insects move, there is much higher risk near where EAB is currently present. Second, the road network is apparent in the risk map, especially where traffic is higher, such as along major interstates and US highways. Third, population density matters with risk: as more people move about, there is an increasing chance of EAB moving with the humans. And finally, quarantined counties show up as important in the risk map (Fig. 4). Once EAB is initially detected in the county, the entire county is quarantined. When that happens, it is legal to freely move wood materials around within the county and between adjacent quarantined counties. Thus, the risk within quarantined counties increases. 365 Summarizing model outputs Of the 255 outlier locations, 82 (32%) fell in our highest risk class (extreme) that primarily captures those zones very near the core with high risk from both the ride and flight models, 76 (30%) fell in the high classes, 89 (35%) in the medium classes, 8 (3%) in the low classes, and 0 in the least class (Figs. 4, 5). In comparing these percentages of positives to area within those classes, 97% of the positives were captured in land representing the medium, high, and extreme classes which occupy only 36.5% of the modeled part of the state, and 62% of the positives occurred in the high and extreme classes, which occupy only 14% of the state (Fig. 5). Although these are promising aspects of model outputs, we realize that these statistics do not really validate the model because our model results are artifacts of a weighting scheme influenced by EAB-positive locations, even though we used RandomForest’s predictor importance as a guide to weight the ride model. The models of this type are inherently difficult to validate because we do not have future data available and even the present data set is far from ideal. When long-distance chance dispersal events dominate as in the case of EAB spread, accurate predictions of spread rate are very difficult or impossible to calculate (Lewis 1997). We therefore ascribe to the philosophy that establishing confidence in the usefulness of the model, especially its relevance to managers, is more important than traditional model validation, which in this case is practically unachievable. Therefore our model aims to provide a better vision of risk for the Fig. 5 Percentages of EAB positive finds and land occupied by risk class from the model shown in Fig. 4 123 366 managers on the ground that are presently faced with the unenviable facts of dwindling budgets, uncertain information, and continued decimation of ash resources. However, we do demonstrate the utility for our risk model by comparing it to known EAB positives and to a statistical imputed map of EAB risk using RandomForest. The RandomForest model produced a map (Supplementary Fig. 6) that identified only 24% of the known outlier positives, while our SHIFT model map of medium to extreme risk identified 97% of the known outlier positives. The RandomForest map also modeled 81% of Ohio as having the least risk for colonization, an underestimate compared to the SHIFT map, which modeled 37% of the state as having the least risk. Thus, the SHIFT model delineates the locations of known EAB positives quite well, and we also have the added information on gradations of risk from low to high giving better information to the managers on the ground. We also tested how well the model, when using the 2005 EAB-occupied map, predicted subsequent positive EAB locations in 2006 and 2007. While the 2005 model captured the trends fairly well for the high to extreme classes (41% of the total EAB positives in 2006 and 2007 captured), a higher proportion of positives fell in the medium risk class (54% of the total; Fig. 6). However, this analysis suffered from the following drawbacks that were mostly beyond our control: (1) Unknown sampling distribution for positives before 2006; (2) Quarantined counties were the same in both 2005 and 2007 models; (3) Decrease in sampling intensity over the years due to lack of resources; and (4) Uncertain lag time between infestation of ash and detection of EAB. Fig. 6 Outlier EAB positive finds (2006–2007) classified by the 2005 combination of the IFM and IRM 123 Landscape Ecol (2010) 25:353–369 However, this analysis points to the need for managers to also consider the medium risk class important especially in proximity to other risk factors, including their individual expertise and experience. Because the extreme class results from the addition of relatively high model scores from both IFM and the IRM models, it will be infested by EAB very quickly due to its proximity to the core infested zone boundary. Therefore, the extreme class is not as valuable for managers compared to the high and medium risk classes. Our narrowing down the possibility of infestation to about 25% of the total area (i.e., the medium and high risk classes) re-emphasizes the utility of our modeling effort. Our EAB SHIFT model outputs also tend to be consistent with the ongoing spread of the insect as reported by the Cooperative Emerald Ash Borer Project (http://emeraldashborer.info/surveyinfo.cfm). It places a level of high risk around the cities of Columbus and Dayton/Cincinnati where population, traffic density, and wood product industries are key influences. In 2006, 24 (22%) and in 2007, 26 (46%) of the outlier positive trees were found either in Columbus and Dayton/Cincinnati or in the neighboring counties where these cities are located. To explore the importance of the primary traffic patterns in relation to the spread of outlier colonies, eleven of the major routes from Detroit, MI, to major cities in Ohio were subset from the total road network; these routes also included one turn only onto an adjoining highway. The analysis of positives within various buffers around this subset of roads revealed that 52% fell within 1 km, 64% fell within 2 km and 81% fell within 4 km of these few roads that represent only 34.7% of the total road length Landscape Ecol (2010) 25:353–369 367 Fig. 7 Ratio of proportion of detection trees to proportion of EAB positives within various distances from highways used in the modeling (Fig. 2a). This result highlights the importance of the role of major highways that are in a connected road network in spreading EAB. Assessment of positives versus detection trees An analysis of the ratios of EAB positives to detection trees validates the extra importance of roads to EAB risk of spread. Even though a higher proportion of positive trees fell within 2 km of the major highways than beyond 2 km, an even greater proportion of the detection trees were located beyond 2 km of a major highway. For example in 2006, 46% of the detection trees fell within 2 km of the highways, but 83% of the positives were within that zone; in 2007, 52% of detection trees and 74% of positives were within the 2-km zone. Thus, overall, there is a 4.5-fold probability for a given ash tree to be attacked by EAB if the tree is within 2 km of a highway as compared to a distance of 4-6 km (Fig. 7), highlighting the importance of road networks in the spread of EAB. Conclusions There is a great deal of ash resource in the eastern United States, especially in the northern half. The EAB is just now entering zones of extremely high amounts of available ash—northeast Ohio, northwest Pennsylvania, and southwest New York. Even in the zones of lower ash availability, like northwest Ohio and northeast Indiana, plenty of ash is available because of high ash per unit of forest in small woodlots, riparian woods, small wetlands, and field borders. We estimate the expansion of the front from 1998 to 2006 to be roughly 20 km per year. This rate of expansion would necessarily have to include both the biological dispersal capacity of the insect and some short-distance movement assisted by humans (e.g., on or in vehicles, plant material, wood material). The stratified nature of the dispersal makes it more difficult to identify a specific front, increasingly so with time. We believe that our hybrid modeling approach adequately captures the dynamics of the EAB spread. The spatially explicit cell-based approach takes into account landscape heterogeneity that mathematical models of spread ignore, and by a novel combination of the insect’s flight characteristics and humanfacilitated movement, addresses both short and long range dispersal. It results in a map of spread that estimates risk areas over approximately the next 2–4 years with much better accuracy than simple imputed statistical maps. We are able to outline degrees of risk in our maps that agree reasonably well with the positive EAB locations so far. This mapping effort should help managers better anticipate future risk from EAB based on uncertain information by locating areas of higher risk and thus allow them to focus where infestations are most likely to occur. It may also help state and county agencies in the placement of a limited number of traps or detection trees, or in sample plot design for researchers. In addition, our approach can be applied to other regions, although reassessment of the core area and re-weighting of the insect ride components may be needed. To sum, we hope our modeling effort results in a better understanding of the risk associated with the 123 368 spread of this destructive insect, and, better informed decisions can be made to detect, monitor and slow the spread. Acknowledgments Thanks to Lindsey Vest of the Ohio Department of Agriculture for providing data on EAB infestation locations. We are grateful to Elizabeth LaPoint from the FIA GIS Support Center for overlaying the Ohio FIA plots with the Ohio GAP data to estimate ash BA. We thank the Ohio Center for Mapping, especially Lawrence Spencer, for creating and providing the Ohio GAP data. We thank Mary Brown, FHWA Office of Highway Policy Information, for providing traffic data. Thanks to Dan Kashian, Denys Yemshanov, John Pedlar and Rueben Keller for the friendly reviews and the two anonymous reviewers, and John Stanovich for the statistical review. We thank Lucy Burde for editing the manuscript. This work was funded in part by the PREISM Program of the USDA (awarded to JMB & LRI). This is publication No. 200X-XXX from the University of Toledo Lake Erie Center. References Bauer LS, Liu HP, Gould J, Reardon R (2007) Progress on biological control of the emerald ash borer in North America. Biocontrol News Inf 28:51N–54N BenDor TK, Metcalf SS (2006) The spatial dynamics of invasive species spread. Syst Dyn Rev 22:27–50 BenDor TK, Metcalf SS, Fontenot LE, Sangunett B, Hannon B (2006) Modeling the spread of the emerald ash borer. Ecol Model 197:221–236 Bossenbroek JM, Kraft CE, Nekola JC (2001) Prediction of long-distance dispersal using gravity models: zebra mussel invasion of inland lakes. Ecol Appl 11:1778–1788 Buck J, Marshall J (in press) In consideration of a secondary dispersal pathway for emerald ash borer adults. Great Lakes Entomol Clark JS, Lewis M, Horvath L (2001) Invasion by extremes: population spread with variation in dispersal and reproduction. Am Nat 157:537–554 Clark JS, Lewis M, Mclachlan JS, Hillerislambers J (2003) Estimating population spread: what can we forecast and how well? Ecology 84(8):1979–1988 Drake JM, Lodge DM (2006) Allee effects, propagule pressure and the probability of establishment: risk analysis for biological invasions. Biol Invasions 8:365–375 Federal Highway Administration USDoT (2007) The national highway planning network (NHPN). http://www.fhwa. dot.gov/planning/nhpn/ (Last Accessed 9–20/2007) Gould J (2007) Proposed release of three parasitoids for the biological control of the emerald ash borer (Agrilus planipennis) in the continental United States. Environmental Assessment, April 2007. United States Department of Agriculture Plant Protection and Quarantine, Otis ANGB, MA Gregory PH (1968) Interpreting plant disease dispersal gradients. Annu Rev Phytopathol 6:189–212 123 Landscape Ecol (2010) 25:353–369 Haack RA (2006) Exotic bark- and wood-boring coleoptera in the United States: recent establishments and interceptions. Can J For Res 36:269–288 Hastings A, Cuddington K, Davies K, Dugaw C, Elmendorf S, Freestone A et al (2005) The spatial spread of invasions: new developments in theory and evidence. Ecol Lett 8:91–101 Hengeveld R (1989) Dynamics of biological invasions. Chapman and Hall, London Herms DA, McCullough DG, Smitley DR (2004) Under attack. American Nurseryman (October), pp 20–26 Hilborn R, Mangel M (1997) The ecological detective: confronting models with data. Princeton University Press, Princeton Hudak AT, Crookston NL, Evans JS, Hall DE, Falkowski MJ (2008) Nearest neighbor imputation of species-level, plotscale forest structure attributes from LiDAR data. Remote Sens Environ 112:2232–2245 Iverson LR, Prasad AM, Schwartz MW (1999) Modeling potential future individual tree-species distributions in the Eastern United States under a climate change scenario: a case study with Pinus virginiana. Ecol Modell 115:77–93 Iverson LR, Schwartz MW, Prasad A (2004) How fast and far might tree species migrate under climate change in the eastern United States? Glob Ecol Biogeogr 13:209–219 Iverson LR, Prasad A, Bossenbroek J, Sydnor D, Schwartz MW (2008a) Modeling potential movements of an ash threat: the emerald ash borer. In: Pye J, Raucher M (eds) Advances in threat assessment and their application to forest and rangeland management. www.threats. forestencyclopedia.net Iverson LR, Prasad AM, Matthews SN, Peters M (2008b) Estimating potential habitat for 134 eastern US tree species under six climate scenarios. For Ecol Manage 254:390–406 Kot M, Lewis MA, van den Driesshe P (1996) Dispersal data and the spread of invading organisms. Ecology 77:2027– 2042 Kot M, Medlock J, Reluga T, Walton BD (2004) Stochasticity, invasions, and branching random walks. Theor Popul Biol 66:175–184 Levins R (1966) The strategy of model building in population biology. Am Sci 54:421–431 Lewis MA (1997) Variability, patchiness, and jump dispersal in the spread of an invading population. In: Tilman D, Kareiva P (eds) Spatial ecology: the role of space in population dynamics and interspecific interactions. Princeton University Press, Princeton, pp 46–69 Liebhold AM, Tobin PC (2008) Population ecology of insect invasions and their management. Annu Rev Entomol 53:387–408 Lonsdale WM (1999) Global patterns of plant invasions and the concept of invasibility. Ecology 80:1522–1536 Miles PD, Brand GJ, Alerich CL, Bednar LR, Woudenberg SW, Glover JF, Ezzell EN (2001) The forest inventory and analysis database: database description and users manual version 1.0. GTR NC-218. North Central Research Station, USDA Forest Service, St. Paul, MN, p 130 Molofsky J, Bever JD (2004) A new kind of ecology? Bioscience 54(5):440–446 Landscape Ecol (2010) 25:353–369 Muirhead JR, Leung B, van Overdijk C, Kelly DW, Nandakumar K, Marchant KR, MacIsaac HJ (2006) Modelling local and long-distance dispersal of invasive emerald ash borer Agrilus planipennis (Coleoptera) in North America. Divers Distrib 12:71–79 Neubert MG, Caswell H (2000) Demography and dispersal: calculation and sensitivity analysis of invasion speed for structured populations. Ecology 81:1613–1628 Neubert MG, Kot M, Lewis MA (2000) Invasion speeds in fluctuating environments. Proc R Soc Lond Ser B 267:1603–1610 Okubo A, Levin SA (1989) A theoretical framework for data analysis of wind dispersal of seeds and pollen. Ecology 70:329–338 Peters MP, Iverson LR, Sydnor TD (in press) Emerald ash borer (Agrilus planipennis): towards a classification of tree health and early detection. Ohio J Sci Poland TM, McCullough DG (2006) Emerald ash borer: invasion of the urban forest and the threat to North America’s ash resource. J For 104(April/May):118–124 Prasad AM, Iverson LR (2003) Little’s range and FIA importance value database for 135 eastern US tree species. http://www.fs.fed.us/ne/delaware/4153/global/littlefia/ index.html Prasad AM, Iverson LR, Liaw A (2006) Newer classification and regression tree techniques: bagging and random forests for ecological prediction. Ecosystems (9):181–199 Prasad AM, Iverson LR, Matthews SN, Peters MP (2007) A climate change atlas for 134 forest tree species of the Eastern United States [database]. Northern Research Station, USDA Forest Service, Delaware, Ohio. http://www.nrs.fs.fed.us/atlas/tree Schwartz M (1992) Modelling effects of habitat fragmentation on the ability of trees to respond to climatic warming. Biodivers Conserv 2:51–61 Schwartz MW, Iverson LR, Prasad AM (2001) Predicting the potential future distribution of four tree species in Ohio, USA, using current habitat availability and climatic forcing. Ecosystems 4:568–581 Shigesada N, Kawasaki K (2002) Invasion and range expansion of species: effects of long-distance dispersal. In: Bullock JM, Kenward RE, Hails RS (eds) Dispersal ecology. Blackwell, Malden, pp 350–373 369 Shigesada N, Kawasaki K, Takeda Y (1995) Modeling stratified diffusion in biological invasions. Am Nat 146: 229–251 Siegert NW, McCullough DG, Liebhold AM, Telewski FW (2007) Resurrected from the ashes: a historical reconstruction of emerald ash borer dynamics through dendrochronological analyses. In: Mastro V, Reardon R, Parra G (eds) Proceedings of the emerald ash borer/Asian longhorned beetle research and technology development meeting. FHTET-2007-04. USDA Forest Service, Morgantown, pp 18–19 Skellam JG (1951) Random dispersal in theoretical populations. Biometrika 38:196–218 Smitley D, Tavis T, Rebek E (2008) Progression of ash canopy thinning and dieback outward from the initial infestation of emerald ash borer (Coleoptera: Buprestidae) in southeastern Michigan. J Econ Entomol 101(5):1643–1650 Snyder RE (2003) How demographic stochasticity can slow biological invasions. Ecology 84(5):1333–1339 Sydnor TD, Bumgardner M, Todd A (2007) The potential economic impacts of emerald ash borer (Agrilus planipennis) on Ohio, US communities. Arboric Urban For 33(1):48–54 Taylor RAJ, Bauer LS, Miller DL, Haack RA (2005) Emerald ash borer flight potential. In: Mastro V, Reardon R (eds) Emerald ash borer research and technology development meeting. FHTET-2005-15. USDA Forest Service, Morgantown, pp 15–16 Taylor RAJ, Poland TM, Bauer LS, Windell KN, Kautz JL (2007) Emerald ash borer flight estimates revised. In: Mastro V, Reardon R, Parra G (eds) Proceedings of the emerald ash borer/Asian longhorned beetle research and Te. FHTET-2007-04. US Department of Agriculture Forest Service, Forest Health Technology Enterprise Team, Morgantown, pp 10 Wei X, Wu YUN, Reardon R, Sun T-H, Lu M, Sun J-H (2007) Biology and damage traits of emerald ash borer (Agrilus planipennis Fairmaire) in China. Insect Sci 14:367–373 West Virginia Division of Forestry (2007) Pest threatens west Virginia’s ash trees. http://www.wvforestry.com/eab.cfm (Last Accessed 2/4/2008) With KA (2002) The landscape ecology of invasive spread. Conserv Biol 16:1192–1203 123