We constructed our network model using a broad-scale data set from the U.S. National Recreation Reservation Service (NRRS), which operates an online reservation system for campgrounds and related facilities operated by the U.S. Forest Service, National Park Service, Bureau of Land Management, and other federal agencies. The available data spanned a period of greater than five years (January 2004-September 2009) and documented ~7.2 million visitor reservations, including reservations from Canada, at campgrounds and recreational facilities throughout the continental U.S. Although this is a large data set, there are many public and private campgrounds outside the NRRS system. Hence, we assumed the data to be a representative sample of all camper travel.

Each reservation record documented the visitor’s origin location (i.e., by ZIP code or Canadian postal code), the destination campground, and the date of the visit. The initial processing of the NRRS data set is described in greater detail in Koch et al. (2012). Notably, the data did not indicate camper firewood usage, so our model instead depicts the general travel behaviour of all campers as documented in the data, a consistent proportion of whom we assumed were carrying potentially infested firewood (for additional discussion regarding this proportion, see Haack et al. 2010; Jacobi et al. 2011; Koch et al. 2012). We should also note that we ignored the visit date for this study, although we acknowledge that a visit’s timing (e.g., in terms of a particular season) can affect the likelihood of pest emergence from firewood and subsequent establishment.

For each individual reservation, we calculated the geographic (Euclidean) distance between the visitor’s origin location (i.e., the centroid of the polygon depicting the visitor’s ZIP or postal code) and the destination campground. We then parsed the data into a set of unique pathway segments. Some data aggregation was necessary during this step to ensure model tractability. Conceptually, our network model is formulated as a first-order transition matrix (Karlin and Taylor 1975) that estimates the probability of travel between every possible origin-destination combination. Each reservation record in the NRRS data represented a single trip of some specified distance between a pair of origin and destination locations. Altogether, the data featured >50, 000 visitor origin locations and >2500 destination locations, which translated to >973, 000 pairwise combinations involving at least one visitor reservation (i.e., at least one trip) during the study period. To build a transition matrix using this many pairwise combinations was computationally impractical, so we aggregated the data at a coarse spatial resolution (into 16 × 16 km cells), ending up with a networked set of ~14, 000 unique map units. These cells provided complete spatial coverage of our study area, which included both the continental U.S. and most of Canada.

Any two map cells, designated i and j, in this networked set were connected by a unique pathway segment, ij. Each ij had a value, m_ij, representing the number of trips from i to j documented in the aggregated NRRS data during the study time period t (i.e., January 2004-September 2009), and another value, m_ji, representing the number of return trips from j back to i. We constructed a matrix, M, from these m_ij and m_ji values. The size of the matrix was n × n, corresponding to the number of unique map units in the aggregated data (i.e., n = ~14, 000):

      (1)

The non-diagonal elements of M were the m_ij and m_ji values for the pathway segments connecting each (i, j) pair of map cells. The diagonal elements of M, where i = j, were set to 0.

Based on M, we constructed a probability matrix, P_t, of camper travel (and by extension, transport of potentially infested firewood) along the pathway segments during the study period t. We assumed that the probability, p_ij, of travel along any given segment ij was linearly related to the number of trips from i to j during period t:

p_ij = m_ij λ_t,      (2)

where λ_t is a scaling parameter. Ideally, the parameter λ_t would define the total likelihood (i.e., over the study period t) of camper transport of pest-infested firewood from i to j in the pathway matrix. In fact, a precise value of λ_t would be necessary to estimate an exact probability of a pest being moved from location i to j. Our objective, however, was not to derive exactly estimated output probabilities, but the more conservative goal of testing the relative sensitivity of model outputs to changes in the inputs, such that the output probabilities simply serve as a relative measure of risk. For this purpose, we only needed λ_t to be sufficiently small to ensure that the sum of the p_ij probabilities in each P_t matrix row was below 1.

Each row of P_t included another variable, p_{i term}, representing the probability that no camper travel (i.e., to any j) would originate at i:

      (3)

If p_{i term} was equal to 1 for any row, then the map cell i associated with that row did not serve as a point of origin in the model. This did not preclude the cell from serving as a potential destination j. The probability matrix P_t of camper travel along each pathway segment (i.e., during period t) was thus specified as follows:

      (4)

We used P_t as the basis for 2 × 10⁶ stochastic pathway simulations of camper travel from each cell i to other map cells. Setting i as the origin location, the model extracted the vector of travel probabilities associated with i (i.e., the probabilities in row i) from P_t in order to simulate subsequent camper travel from i through the pathway network to other map cells. Briefly, in each simulation run for a given cell i, the model performed a uniform random draw against the associated row of probabilities in P_t in order to select the next map cell defining the path of an individual trip. The simulation of this path continued until reaching a final destination map cell (i.e., a cell with no outgoing travel), or instead, when a terminal state (i.e., no further travel) was selected based on the p_{i term} value. Regardless of the number of individual segments comprising a simulated path, we assumed that the path was completely traversed within the study period t. Consequently, for a given pathway segment ij, the probability of camper travel from location i to location j during period t was estimated as follows:

φ_ij = M_ij/M,      (5)

where M_ij is the number of times pathway travel from location i to j was simulated to occur, and M is the total number of simulations (M = 2 × 10⁶ in this study).

The network model permitted us to generate maps where each individual cell in the study area could be set as the origin or destination location of interest. While cell-specific maps might have value when addressing specific invasion scenarios (e.g., identifying the probable source locations for a single cell found to be invaded), they do not provide a comprehensive picture of invasion risk for use in setting management priorities. So, we chose to summarize the model results in a broader context. In a previous analysis of the NRRS data (Koch et al. 2012), we found that the majority (~53%) of camping trips involved travel distances of less than 100 km. However, we also found that 10% of trips involved distances of greater than 500 km. This is consistent with other research (Haack et al. 2010; Jacobi et al. 2011; USDA Animal and Plant Health Inspection Service 2011) suggesting that campers frequently transport firewood across state or provincial borders. Furthermore, much of the regulatory decision-making with respect to invasive species, including the implementation of firewood restrictions to limit forest pest spread, happens at the state or provincial level (Filbey et al. 2002). For these reasons, we opted to develop model outputs that characterized those locations (i.e., map cells) outside a state or province of interest that were most strongly linked to locations inside the state or province.

For each U.S. state and Canadian province, we generated a map where, for each cell i outside the target state (province), we summed the model-derived travel probabilities (i.e., the φ_ij values) for all pathways between i and any destination cell j within the target state (province). Essentially, each map depicts the most likely external source locations if a firewood-transported forest pest were to be discovered within the state (province) of interest. (Note that for the Canadian provinces, we assume that campers returning from U.S. locations may be transporting forest pests, despite border biosecurity measures.) A principal feature of these maps is that they highlight key “bottleneck” locations where surveillance, awareness campaigns, or quarantine procedures have the best potential to be cost-effective in terms of protecting the state (province) of interest from invasion (Hauser and McCarthy 2009). These maps served as a set of baseline maps that we could then compare to the maps generated for a set of sensitivity scenarios, described below.

We performed sensitivity analyses to evaluate how uncertainty in key aspects of the network model’s structure would influence the φ_ij probabilities estimated in the baseline maps. The primary goal of these analyses was to assess whether uncertainties in any of the tested model aspects yielded a set of output maps with spatial patterns that departed in unexpected ways from the patterns depicted in the corresponding baseline maps. We adopted a Monte Carlo simulation approach (Morgan and Henrion 1990; Crosetto and Tarantola 2001; Li and Wu 2006) for the analyses, involving three basic steps: random sampling from an input distribution defined for each aspect of interest; repeated model runs using the randomly sampled values; and summarization of these results for comparison to the baseline maps.

We focused on three model aspects for this study. Our first sensitivity scenario tested the impact of uncertainty in the network configuration by randomly removing up to 30% of nodes (i.e., map cells and their associated vectors of p_ij probabilities) from the network prior to each new model run. Our second sensitivity scenario tested the impact of uncertainty in the scaling parameter λ_t (see Eq. 2). In this case, we added uniform random error within the range ±0.15 to λ_t prior to completing each new model run. Our final sensitivity scenario evaluated the impact of uncertainty in the p_ij values comprising P_t. To test this aspect, we added uniform random error of up to 0.01 to the p_ij values before each new model run. Ideally, we would have repeated the three scenarios with series of bounding values. However, because of the computational complexity of this exercise, we chose a single bounding value for each sensitivity scenario that, based on preliminary work with the model, we were confident would provide meaningful contrast to the baseline scenario.

We completed 2 × 10⁶ pathway simulations for each sensitivity scenario. This process yielded scenario-specific maps of the φ_ij probabilities (see Eq. 5) for each state and province that we could compare to the corresponding baseline maps. For this comparison, we calculated and mapped the differences, in percentage terms, between the output probabilities under the baseline scenario and under each sensitivity scenario. We also calculated a Moran’s I (Moran 1950) value for each difference map. Moran’s I is a measure of the global (i.e., map-wide) spatial autocorrelation of values in a spatially referenced data set (Getis and Ord 1992; Perry et al. 2002). It ranges in value from -1 to 1, with values close to -1 indicating a high degree of dispersion in the data (e.g., a checkerboard spatial pattern), values close to 1 indicating a high degree of clustering, and values near 0 indicating a random spatial pattern.

We employed a simple rank-difference approach for further comparison of the baseline results to those from our third sensitivity scenario dealing with uncertainty in the p_ij values. After calculating the per-cell differences between the baseline and sensitivity maps for each state and province, we converted both the baseline probabilities and the calculated differences into ranks. Each cell in a given map was assigned a rank from 1 to 6 based on a global percentile distribution of the values that occurred in the maps for all states and provinces; cell values that fell within the bottom 25% of this global distribution received a rank of 1, cells in the 25–75% range received a rank of 2, cells in the 75–90% range received a rank of 3, cells in the 90–95% range received a rank of 4, cells in the 95–99% range received a rank of 5, and cells in the top 1% of this distribution received a rank of 6. We ranked each cell twice in this fashion, first according to its percentile value under the baseline scenario, then according to its percentile value in the global distribution of differences under this sensitivity scenario. Next, we calculated the change in rank between the two, which yielded an index value ranging from -5 to 5. Finally, we mapped the index values for each state and province in order to identify any spatial trends.

Across all individual state and province maps created for the baseline scenario, we saw three general spatial patterns of out-of-state (or out-of-province) origin risk. First, as illustrated by the state of Alabama (Fig. 1a), there were cases where the highest origin probabilities (i.e., the per-cell φ_ij values) were primarily limited to a narrow fringe zone surrounding the target state. This pattern was common among states in the southeastern and northeastern U.S. The second general pattern, exemplified by Missouri (Fig. 1b), similarly featured a localized zone of high probabilities around the target state, but this was augmented by several high-probability hotspots associated with major urban areas located outside this zone. For instance, the map for Missouri included hotspots associated with the cities of Chicago (IL), Dallas (TX), Denver (CO), all of which are at least 400 km away. This mixed origin pattern was typical of many states in the U.S. Midwest. In the third pattern, exemplified by Utah (Fig. 1c), there was also a subtle local fringe zone, but the majority of map cells with high probabilities were associated with major urban areas that are fairly distant from the target state. With respect to Utah, the most prominent urban areas (i.e., those displaying uniformly high probabilities) were generally in the western U.S., but a number of eastern U.S. cities also had probabilities in the top 5% of the global distribution of probabilities across all states and provinces. This sort of dispersed pattern of origin risk was exhibited by other states such as Arizona and California that, like Utah, feature several popular recreational destinations (e.g., national parks) that draw campers from throughout the continent.

The maps for the Canadian provinces, especially the more populous ones (i.e., Quebec, Ontario, Alberta, and British Columbia), usually exhibited a dispersed spatial pattern similar to that observed for Utah, although the origin probabilities in each map were low compared to their U.S. counterparts. Indeed, the highest probabilities in these maps very rarely reached the top 10% of the global distribution. The only exceptions were a few high-probability (top 5%) cells in each of the maps for the most populous Canadian provinces, which appeared to be associated with specific recreational destinations such as Grand Canyon National Park (AZ) and Zion National Park (UT).

Two factors seem to govern the observed patterns of origin risk. Foremost is the population of the state or province of interest. Because the model is bi-directional (i.e., it also simulates return travel by campers from their destinations to their origin locations), the population of the target state or province influences the range of output probabilities that will be portrayed in its map, particularly the frequency at which high probabilities occur. In short, the maps of heavily populated states and provinces tended to have higher probabilities overall than sparsely populated states and provinces. The second factor is the nature of the recreational opportunities available in the target state or province. States or provinces containing several high-profile recreational destinations displayed a dispersed pattern of origin risk, while states or provinces with numerous but comparatively low-profile recreational destinations usually displayed a more localized origin risk pattern.

Figure 2 shows three example maps of the mean percent difference (i.e., over all model runs) in the φ_ij output probabilities under this sensitivity scenario. The example states correspond to those highlighted under the baseline scenario: Alabama (Fig. 2a), Missouri (Fig. 2b), and Utah (Fig. 2c). Note that only map cells in the top 10% of probabilities under the baseline scenario are shown. For all states and provinces, the percent change in the values of cells in this top 10% typically ranged from -24% to -36%, with a map-wide mean (i.e., across all cells in a given map) near -30%. This level of difference is consistent with the proportion of nodes removed for this scenario. More importantly, as illustrated by the difference maps for all three example states, there was no obvious spatial trend in the pattern of change. Essentially, each map has a salt-and-pepper appearance, such that cells with higher probabilities under the baseline scenario (see Fig. 1) do not appear to have greater (or smaller) percent changes than cells with lower probabilities under the baseline. This observation is supported by the Moran’s I values calculated for all states and provinces, which ranged between -0.016 and 0.021 (median = 0.005) and thus indicated a minimal level of spatial autocorrelation.

Figure 3 shows three example maps of the mean percent difference in the φ_ij output probabilities under this sensitivity scenario. As in the previous scenario, the example states correspond to those highlighted under the baseline scenario: Alabama (Fig. 3a), Missouri (Fig. 3b), and Utah (Fig. 3c). Only map cells in the top 10% of probabilities under the baseline scenario are shown. For all states and provinces, the percent change in the values of cells in this top 10% typically ranged from -15% to 15%, with a map-wide mean slightly greater than 0. This degree of change is consistent with the alterations to λ_t under this scenario.

None of the difference maps for the example states appears to show a strong spatial trend in the pattern of change. There may be a slight tendency for map cells with higher probabilities under the baseline scenario to exhibit small but positive percent changes, while cells with lower probabilities under the baseline may tend toward small negative changes. However, the existence of a subtle spatial trend is not really supported by the Moran’s I values for each state and province, which ranged from -0.005 to 0.133 (median=0.008). Indeed, with the exception of a single Canadian province, Nova Scotia, which has relatively few connections to outside locations in our network model, none of the maps had an I value greater than 0.0324. Similar to the previous sensitivity scenario, this suggests a low level of spatial autocorrelation.

Figure 4 shows three example maps of the rank differences (see Methods) between the baseline scenario and this sensitivity scenario, which tested the effects of adding uniform random error to the p_ij probabilities. The maps display an obvious spatial trend in the pattern of change. (Again, only cells in the top 10% of probabilities under the baseline scenario are shown.) This trend can be summarized in general terms: cells with comparatively higher probabilities under the baseline scenario were more likely to exhibit a decrease or no change in rank under this sensitivity scenario, while cells with comparatively lower probabilities under the baseline were more likely to exhibit an increase in rank. Thus, for states and provinces like Alabama (Fig. 4a), where the highest origin probabilities under the baseline scenario were primarily limited to a localized zone adjoining the state or province (also see Fig. 1a), any decreases in rank were also usually limited to this zone. Alternatively, in states and provinces like Missouri (Fig. 4b) and Utah (Fig. 4c), where high origin probabilities under the baseline were regularly associated with major urban areas outside the fringe zone, it is primarily map cells in these areas that exhibited a decrease or, just as commonly, no change in rank. For all states and provinces, the map cells that exhibited increases in rank under this sensitivity scenario were usually found in rural areas or the peri-urban zones (Allen 2003) that surround large cities.

The existence of a spatial trend in the pattern of change is supported by the Moran’s I values calculated for the difference maps of each state and province. With the exception of Nova Scotia, the values ranged from 0.057 to 0.208 (median = 0.135), indicating a degree of positive spatial autocorrelation (i.e., clustering) that varied according to whether the state or province of interest had a well-defined fringe zone of high probabilities under the baseline scenario. Typically, states and provinces with such a zone exhibited larger I values. This variability, however, does not contradict the notion of a generally consistent spatial trend across all states and provinces.