A safety rule approach to surveillance and eradication of biological invasions

Uncertainty about future spread of invasive organisms hinders planning of effective response measures. We present a two-stage scenario optimization model that accounts for uncertainty about the spread of an invader, and determines survey and eradication strategies that minimize the expected program cost subject to a safety rule for eradication success. The safety rule includes a risk standard for the desired probability of eradication in each invasion scenario. Because the risk standard may not be attainable in every scenario, the safety rule defines a minimum proportion of scenarios with successful eradication. We apply the model to the problem of allocating resources to survey and eradicate the Asian longhorned beetle (ALB, Anoplophora glabripennis) after its discovery in the Greater Toronto Area, Ontario, Canada. We use historical data on ALB spread to generate a set of plausible invasion scenarios that characterizes the uncertainty of the beetle’s extent. We use these scenarios in the model to find survey and tree removal strategies that minimize the expected program cost while satisfying the safety rule. We also identify strategies that reduce the risk of very high program costs. Our results reveal two alternative strategies: (i) delimiting surveys and subsequent tree removal based on the surveys' outcomes, or (ii) preventive host tree removal without referring to delimiting surveys. The second strategy is more likely to meet the stated objectives when the capacity to detect an invader is low or the aspirations to eradicate it are high. Our results provide practical guidelines to identify the best management strategy given aspirational targets for eradication and spending.


Introduction
Programs that are designed to stop the spread of invasive species often involve the allocation of resources to survey and eradicate the invading individuals [1][2][3][4][5][6]. Although a considerable share of the available resources is often devoted to surveys [7], such efforts rarely reveal complete information about the presence of the species of interest. The resulting uncertainty about the spread and extent of the invader means that the costs-and likelihood of success-of the a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 management actions to control the invasion are also uncertain. This creates problems with allocating resources appropriately to control the invasions, because critical management decisions (e.g., when and where to prioritize surveys or eradication efforts) must be made under uncertain expectations of the likelihood and outcomes of invasion. A common approach in this situation is to estimate the expected cost of eradicating the species based on probabilistic expectations of the species' spread in the area of concern. These expectations are uncertain but can be represented as a large set of plausible stochastic scenarios that help estimate the bounds of uncertainty on those expectations [8][9][10][11][12][13][14]. However, the true costs of the actions that ultimately must be taken to manage the invasion remain unknown, because the decisions based on the expected costs could be wrong and lead to extremely high management costs.
Detection and control of biological invasions can be greatly improved through application of spatial-dynamic optimization models that predict economically optimal strategies for surveillance and eradication of invasive species. The models depict the allocation of resources to control an invasion as an optimization problem, with some important parameters and decision variables depicted in temporal and spatial domains. Epanchin-Niell et al. [15] developed a dynamic model of pest colony establishment and growth and designed optimal long-term equilibrium surveillance and eradication programs to minimize program costs. They used the model to optimize long-term surveillance effort across heterogeneous landscapes subject to region-wide surveillance budgets. Hauser and McCarthy [16] proposed a static model to optimize one-time surveillance effort across multiple sites when a species' presence is uncertain prior to detection and probability of occurrence varies across sites. In contrast to the equilibrium analysis of Epanchin-Niell et al. [15], the static model of Hauser and McCarthy [16] is appropriate for optimizing surveillance when many local populations are thought to have established prior to the initiation of a surveillance program. Horie et al. [17] and Yemshanov et al. [14] developed models to optimize one-time surveillance effort across multiple sites given uncertainty about the extent (rather than simply the presence) of the infestation in each site. They handled this uncertainty by splitting the management decision into two stages. In the first stage, sites are selected for surveillance given their expected levels of infestation, and in the second stage, eradication treatments are prescribed within the surveyed sites contingent on the levels of infestation found. The objective is to minimize the expected growth of the infestation subject to the total budget for surveillance and treatment. Epanchin-Niell et al. [18] developed a mechanistic bioeconomic model that relates surveillance intensity and invasion size to probabilities of detection and control. Their model determined, in a geographic domain, the optimal investment in surveillance, in terms of the numbers and distributions of traps, to minimize the total invasion impact. Moore and McCarthy [19] proposed a model that optimizes the allocation of surveillance efforts in both spatial and temporal domains and accounts for stochastically varying detection rates in geographical space.
In this study, we address a problem in which a decision maker must select a program for delimiting surveys and eradication that minimizes overall program costs and attains a desired level of eradication success despite uncertainty about the current and future extent of an invasion. Similar to Horie et al. [17] and Yemshanov et al. [14], we split the management decision into two stages representing the placement of surveillance in the first stage and the intensity of treatment given the outcome of surveillance in the second stage. Rather than attempting to minimize expected invasion expansion, we include probabilistic constraints for attaining eradication success as a way to limit the damage from the invader populations that have been established in the area of interest. These constraints are consistent with a safety-rule approach to addressing uncertainty in environmental regulation [20][21][22]. A safety rule includes a risk standard representing a minimum probability of attaining a desirable environmental outcome (e.g., eradication of an invasive species population). The safety rule includes a probability that management action(s) will not achieve a desired outcome. This probability value (also known as margin of safety) represents the decision-maker's aversion to uncertainty about the attaining the risk standard.
We apply our safety-rule approach to the case of managing the invasion of a forest pest, the Asian longhorned beetle (ALB), Anoplophora glabripennis (Motschulsky), in the Greater Toronto Area of Ontario, Canada. This woodborer is a native of China and the Korean Peninsula and is listed among the world's worst invasive species [23]. In invaded landscapes, its preferred host is maple (Acer spp.); other suitable hosts include, but are not limited to, birches (Betula spp.), poplars (Populus spp.), willows (Salix spp.), and elms (Ulmus spp.) [24][25][26][27][28][29]. Eastern North America seems especially vulnerable to infestations by this beetle because of its abundance of maple and of suitable habitat conditions [30]. Breeding populations of ALB have been found at several locations in eastern North America [28,31,32] and in Europe [33][34][35]. Most of these introductions involved beetles that had originated from the species' native range [36][37][38]. In most countries outside of its native range, discovery of ALB leads to the implementation of an eradication program [39,40]. Such programs typically feature a delimitation phase of the invaded area, which leads to the establishment of a quarantine zone to limit the spread of the beetle as well as a treatment phase [41]. Removal and destruction of all infested trees and removal or treatment of high-risk suitable host trees, a strategy that takes advantage of ALB's slow rates of spread [42,43], represents an effective method to eradicate ALB. Indeed, this strategy has led to successful eradications in North America [44,45].

Scenario-based invasion management model
We develop a two-stage, spatial optimization model in which uncertainty about the presence and extent of the invader (i.e., ALB or another forest pest) is represented by a set of probabilistic scenarios. We assume at the beginning of the first stage that the invader is known to be present in a single area and a quarantine zone is established around that invaded area. The area under quarantine contains J sites where the invader may also be present, but which of these sites are actually infested is unknown. Defining a quarantine zone ("regulated area" hereafter) surrounding known infestations is a common phytosanitary practice aimed at containing expanding invader populations [41]. The size of the area J is defined by decision makers based on expectations of future spread and new introductions. Each site j, j 2 J, in the regulated area J contains N j suitable host trees that can be used by the invader to complete its development and reproduction (see Table 1 for a summary of model parameters). The proportion of host trees at site j that are infested is θ j , θ j 2 [0,1]. While the true proportion of infested trees, θ j , is unknown, the θ j value at each survey site is estimated based on the invader's historical patterns of spread and the numbers of infested trees found in surveyed sites during previous surveillance campaigns. These estimates are used to develop a large set of probabilistic scenarios, S, where each site j in scenario s, s = 1,. . ., S, is characterized by the proportion, θ js , of trees that are infested at the site under that scenario. We assume each invasion scenario s has an equal probability of occurrence, 1/S. The model has one set of decision variables in each stage. In the first stage, the binary decision variables x j , x j 2{0,1} and j 2 J, represent whether or not to select site j for survey. While sites are selected for survey with knowledge that they could potentially be invaded, the selection is done without knowing which sites are actually invaded, and if so, to what extent. In the second stage, the decision variables R js , j 2 J and s 2 S, represent the number of trees to remove in each site j under scenario s, which is contingent on the proportion of infested trees that are found to be present under the scenario.
The manager's objective is to minimize the expected cost of survey and tree removal in the regulated area J over S invasion scenarios, i.e.: where c j and t j are the costs of inspecting and removing host trees in site j and β, β 2 [0;1], is the proportion of a site's area (in our formulation, equivalent to the proportion of host trees at a site) that has been surveyed. Note that our model considered a single level of survey intensity at all sites. Ideally, the proportion of a site's area that is surveyed (β) should be a decision variable, but the model formulation would become non-linear in that case. While it may be possible to reformulate the model to address this non-linearity it would significantly increase the computational complexity of the problem. We define two constraints on the number of trees that can be removed from each site. The number of removed trees cannot exceed the total number of trees at a surveyed site j: Trees are only removed (R js ! 0) at sites which are selected for surveys (x j = 1). If a tree is infested, it can be found with the detection rate γ, γ 2 [0,1], and each tree that is found to be infested must be removed. Because the survey covers only a proportion β of a site, the minimum number of infested trees to remove from site j in scenario s is: Eq 3 implies that all detected infested trees must be removed. A special case with β = 0 describes a situation when a decision-maker chooses to skip the delimiting survey stage and proceed with preventive tree removal at a selected site instead. We next define a probabilistic constraint for successful eradication of the pest from all of the sites in area J under a given scenario s. The probability that one or more of the remaining host trees in area J is infested after R js trees have been removed at sites j, j 2 J, is: The term ð1 À y js ½ð1 À bgÞ=ð1 À bgy js ÞÞ ðN j À R js Þ defines the probability that the remaining trees in site j are not infested in a scenario s. Then, the product sign in Eq 4 represents the probability that the remaining trees in all of the sites are not infested, and one minus that probability equals the probability that removal fails to eradicate the infestation from area J. The derivation of Eq 4 can be found in S1 File. A requirement that the probability of eradication failure be below a chosen threshold, D, is: ð1 À y js ½ð1 À bgÞ=ð1 À bgy js ÞÞ Rearranging terms, the probabilistic constraint for successful eradication of the pest is: ð1 À y js ½ð1 À bgÞ=ð1 À bgy js ÞÞ where the threshold d, d = 1 -D, defines the minimum probability for successful eradication in area J. Eq 6 is linearized with respect to the decision variables R js : ðN j À R js Þlnð1 À y js ½ð1 À bgÞ=ð1 À bgy js ÞÞ Eq 6 represents a decision-maker's aspiration to satisfy the eradication success threshold for each and every invasion scenario. When uncertainty about the extent of the invasion is significant, achieving successful eradication in all possible scenarios may be cost-prohibitive because most of the sites would require tree removal. We address this issue by defining a margin of safety p, representing the minimum proportion of the scenarios that must satisfy the eradication success threshold d. To create a constraint for the margin of safety, we first define a binary indicator variable, g s , g s 2{0,1}, that specifies, for each scenario s, whether eradication succeeds with probability d (i.e., g s = 1) or fails (g s = 0; the failure condition is defined with a probability of eradication success, d 0 , equal or close to zero. Since the d 0 value is under the logarithm sign it was set to a very low positive value.). We have used the following constraint to set the g s value: Then, a constraint that requires attainment of the margin of safety is: Note that Eqs 8 and 9 are consistent with a safety-rule approach to dealing with uncertainty in environmental regulation (e.g., [20][21][22]). Eq 8 defines whether or not a risk standard is satisfied for each scenario, where the risk standard is the required probability of eradication, d. Eq 9 defines the margin of safety, which is the required proportion of the scenarios that satisfy the risk standard. The margin of safety p can be interpreted as a confidence level for a hypothesis test with the values typically set to 0.9 or 0.95. In summary, the problem is to determine which sites to survey in the first stage and how many trees to remove in the second stage contingent on the outcome of the first stage, in order to minimize the expected cost of surveys and tree removals over S scenarios (Eq 1), subject to upper and lower bounds on tree removals (Eqs 2 and 3) and constraints for satisfying the risk standard for eradication (Eq 8) and the margin of safety (Eq 9).

Controlling the risk of extreme project costs
Due to the uncertainty about the invader's spread, the actual costs of survey and eradication in individual scenarios may vary. Some scenarios could be severe in terms of invasion extent and impact, and therefore may require costly eradication. These scenarios are located in the right tail of the cost distribution (Fig 1), which the decision-maker tries to avoid. The problem of controlling the worst (i.e., most cost-intensive) outcomes essentially amounts to controlling the costs in the right tail of the program cost distribution. This can be achieved with a metric that characterizes the upper percentile of the distribution.
Percentile-based metrics such as maximum loss [46,47], Value-at-Risk [48] and Conditional Value-at-Risk [49,50] have been widely used to quantify risks of extreme losses in many disciplines. In particular, Value-at-Risk (VaR) and Conditional Value-at-Risk (CVaR), also known as expected shortfall or conditional tail expectation (CTE), are commonly used by financial institutions to evaluate the potential for extreme losses in investment portfolios [51][52][53][54][55]. With respect to our invasive species example, VaR α is defined, with a confidence level α, α 2 [0;1], as the value in the distribution of the surveillance and eradication program costs that is exceeded only in (1 -α)×100% of the worst scenarios (Fig 1). In turn, CVaR α , for a confidence level α, is defined as the expected value of the cost distribution over (1 -α)×100% of the worst scenarios, or alternatively, as the expected value above VaR α for confidence level α. Incorporating VaR in an optimization framework is difficult [56] because VaR is a nonconvex function for discrete distributions and may not account for properties of the distribution beyond the confidence level α. The CVaR is more useful for optimization-based models [49,51] because it is a coherent risk metric (cf. Artzner et al. [57]) and continuous with respect to the confidence level α. The most appealing property of CVaR is its convexity with respect to the decision variables [49,56]. Optimization of CVaR with respect to linear decision variables for discrete scenario-based distributions can be expressed by a set of linear equations [56]. In comparison, optimizing VaR in the same problem setting could be numerically intractable.
In our model, we used the CVaR to control the right tail of the program cost distribution and reduce the cost uncertainty. Ideally, both the expected cost value (i.e., the mean value across the entire cost distribution) and the CVaR should be minimized, however this requires a two-objective formulation. We reformulated the objective function equation as a weighted average between the expected cost value and the CVaR α of the program cost, i.e.: where F is the weighting parameter that defines a decision-maker's preferences towards minimizing the expected cost versus minimizing the expected value in the right tail of the cost distribution (i.e., the CVaR α ). For F = 0, the objective function minimizes the expected cost, and for F = 1, the CVaR α is minimized. Eq 10 adopts an approach of combining multiple objectives in a single objective function equation via weighted averaging [58,59], with F and 1-F values representing the objective weights. By altering the F values, a trade-off between the objectives can be explored. The objective function in Eq 1 is linear with respect to the decision variables x j and R js , so we applied a linearized formulation of the CVaR α from [49,60]. For a discrete distribution of S scenarios with equal probability of occurrence, 1/S, the CVaR of the program cost, at a confidence level α, can be approximated with the following equivalent set of S + 1 auxiliary decision variables and S + 1 inequality constraints: where P J j¼1 ðbN j c j x j þ t j R js Þ in Eq 12 is the total program cost in a scenario s, z and w s are auxiliary decision variables and z is a member of a set of real numbers. We have modified the objective function to follow the idea of Eq 10 as: ðbN j c j x j þ t j R js Þ is the expected program cost, subject to model constraints in Eqs 2, 3, 8, 9 and auxiliary constraints in Eqs 12 and 13.
By adding the CVaR α to the objective function equation, the uncertainty in the right tail of the program cost distribution can be controlled and the risk of incurring high program costs can be reduced. Reducing the risk of high program costs increases the expected program costs because more resources will be required to reduce the cost uncertainty. We have explored this aspect by evaluating the optimal solutions for a range of F values between 0 and 1, and report the solutions with F values between 0.5 and 1.0 that yielded the greatest reduction of the cost at CVaR α .
Case study: Assessing the program costs for controlling the outbreak of Asian longhorned beetle (ALB) in the Greater Toronto Area (Ontario, Canada) Assessing the human-mediated spread of ALB in an urban setting. Two infestations of ALB have been found in the Greater Toronto Area (GTA). The first population was discovered in 2003 in Toronto [61], and the second one in 2013 in Mississauga [62]. Discovery of each population led to the implementation of an eradication program [28,62,63]. The regulated (i.e., quarantine) area to manage the first 2003 infestation was declared pest-free in 2013 [44]. The Mississauga quarantine area was about 46 km 2 and outside that of the Toronto quarantine area (Fig 2). Treatment for both eradication programs consisted of the removal of all infested trees as well as host trees categorised as suitable and at high risk of being infested. We applied our model to manage the most recent ALB incursion in Mississauga. Several techniques have been tested in an attempt to improve early detection of ALB [28,64,65], but visual inspection of trees for signs of attack remains the most practical detection method [35,63]. The possibility of new introductions and high costs of eradication necessitate that surveillance activities linked to the eradication of ALB extend beyond the initial regulated area, as well as the potential expansion of the regulated area and an assessment of the total eradication costs were this to happen in the future.
ALB is known to have slow spread rates [42,43]: 80% of the population at a given site is expected to spread less than 300 m per year [66]. Most recent ALB introductions have been attributed to human activities [38]. Growing anecdotal evidence suggests that the pest may hitchhike on slow-moving vehicles that have been previously parked near suitable host trees [32,Turgeon,pers. obs.], similar to the documented spread of another invasive forest insect, the emerald ash borer (Agrilus planipennis Fairmaire) [67]. Local street traffic accounts for a large portion of the movement of people and goods in urban settings, and has been recognized as a proxy for a variety of local economic activities [68]. We used volumes of local road traffic as a measure of activities that could facilitate ALB spread. We utilized a dataset on local street traffic volumes [69] that had been linked to the GTA portion of the ESRI Street Map geospatial database [70,71] to estimate probabilities of ALB movement from previously invaded locations (see description in S2 File). We divided the GTA street network into 400×400-meter blocks, each representing a potential survey site, and then used the local traffic volume data to estimate a matrix of probabilities of ALB movement from block to block via the network. We also adjusted the probability values by the annual projected rates of traffic volume increase, which we estimated from the available traffic volume data over the last 10 years. This matrix was used to simulate 5×10 6 randomized pathways of ALB spread between sites in the study area, and to estimate arrival rates for each of the destination sites (see next section). The ALB spread model then was calibrated to match the historical spread rates of ALB in the GTA as determined from previous survey campaigns prior to the current eradication effort (S2 File).
Parameterizing the optimization model. The optimization model required a large set of plausible invasion scenarios. As noted before, each scenario had an associated set of proportions of trees that are infested, θ js . We generated the invasion scenarios in two steps. First, we used our calibrated pest spread model to estimate the probability of ALB arrival, p j est , for each site j in the study area (Fig 2). For each spread scenario s, we generated a stochastic pattern of invaded sites via uniform random draws against the p j est values. Next, each invaded site j in a scenario s was assigned a number of infested trees that was randomly sampled from an empirical distribution of the number of infested trees. This distribution was based on previous records of historical ALB detections in the GTA prior to the current eradication campaign. The proportion of infested trees, θ js , at a site j in a scenario s was then found by dividing the number of infested trees sampled from the distribution by N j , the total number of host trees at site j. To estimate the total number of host trees (N j ), we first estimated the area of tree cover at each survey site j from the SOLRIS land cover dataset for the GTA [72]. Subsequently, we converted the area of tree cover into a corresponding number of host trees by multiplying by tree density and the host species proportion; the tree density values came from the SOLRIS data, while the estimates of host tree species proportions came from the City of Toronto's Every Tree Counts survey [73], which provided a detailed summary of Toronto's urban forests.
The model also required estimates of the costs of survey (c j ) and tree removal (t j ) as well as the pest detection rate value (γ). Because surveys are conducted within an accessible street network, we assumed equal survey costs on a per-tree basis. We estimated the survey cost from contractor rates paid to do visual tree inspections in previous ALB survey campaigns. This yielded an average survey cost of $6.83 per tree. The cost of tree removal was based on current tree disposal costs in the ongoing ALB eradication program and was set at $1000 per tree. The high cost of tree removal was due to regulatory requirements when disposing of a tree, which require costly chipping operations at designated disposal sites. The baseline probability value of detecting ALB by inspecting a host tree, γ, was set to 0.7. This estimate is based on experience gained during the previous ALB surveillance programs and assumes that inspections are performed by trained personnel [63]. We have also evaluated the model behaviour for a range of detection rates between 0.3 (inspections by untrained personnel [37]) and 0.95 (inspections in idealized conditions). We evaluated optimal solutions with different aspirational targets of eradication success, d = 0.5 and 0.95, and safety margins, p = 1 and 0.95.
Computing bounds on the objective function value. We also assessed the suitable range for the number of stochastic spread scenarios, S, that should be used in the model. Ideally, the optimization model would be supplied with a very large set of invasion scenarios, but the number of scenarios is limited by computational capacity. Optimal solutions based on a finite number of random scenarios provide an approximation of the true optimal solution. To better understand how close our model solutions were to the solution with a complete set of scenarios, we estimated the upper and lower bounds on the optimal objective function value using concepts from Mak et al. [74] and Lee et al. [75]. We computed the bounds on the objective function for an area J of sufficient size to cover the majority of plausible ALB spread patterns in the GTA over a short-term planning horizon (i.e., 3208 sites). This was done to ensure that the number of scenarios was sufficient to depict the variation of possible spread outcomes. The lower bound ( " L) was estimated as the mean of the objective function values for the solutions to 20 independent replicate problems with S scenarios, S = 400, 800, 1200, 1800, 2400 and 3000 scenarios. For each of the 20 solutions based on non-overlapping sets of spread scenarios, obtained to compute the lower bound, we re-computed the objective function value using a set of 6000 scenarios, and then estimated the upper bound ( " U ) as the mean of the objective function values in those 20 sets. The optimality gap was estimated as the relative difference between the upper and lower bounds, i.e.:ð " U À " LÞ= " U [15]. A summary of the model parameters and variables is shown in Table 1. We prototyped the model in SolverStudio [76] and GAMS environments [77] and solved the MIP problem using the GUROBI linear programming solver [78].

Number of spread scenarios and the optimality gap
We estimated the upper and lower bounds on the objective function value for problem solutions with 400, 800, 1200, 1800, 2400 and 3000 spread scenarios ( Table 2). The optimality gap was around 8% for problem solutions with 400 scenarios and around 1% for problem solutions with 2400 scenarios. We also examined the program costs (i.e., averaged over 20 independent replicates) for the optimal solutions based on different numbers of scenarios (Table 3). Specifically, the worst-case and upper percentile (α = 0.95) cost values indicate how well the scenarios depict the most damaging outcomes with the highest cost. As the number of scenarios S exceeds 2400, the worst-case cost value and the CVaR both stabilize, which indicates that further increase of the number of scenarios does not add much information about extreme invasion events. It does not appear that solving problems with more than 2400 scenarios would Table 2

. Upper and lower bounds on the objective function value for differing numbers of spread scenarios. The mean values, "
Uand " L, of the objective function over a specified number of scenarios S, as well as their 95% confidence intervals, were calculated for a base case with the desired probability of successful eradication d = 0.5 in all scenarios (p = 1), the detection rate γ = 0.7 and the proportion of the site covered by surveys β = 1. The objective function minimized the expected program cost (i.e., F = 1) and was computed with sets of 20 independent replicates with increasing numbers of scenarios, S = 400, 800, 1200, 1800, 2400 and 3000.  Table 3. Basic description of the program costs and survey and tree removal components for differing numbers of spread scenarios. The program cost values are mean estimates based on 20 independent replicates with increasing numbers of scenarios, S = 400, 800, 1200, 1800, 2400 and 3000. See Table 1 for the model parameter settings. yield much benefit in terms of precision or accuracy, but would greatly increase computing time.

Scenario-based program costs, $M Survey cost, $M Number of removed trees
The number of scenarios affected the spatial patterns of the surveyed sites. Fig 3a shows the survey selection pattern (averaged for 20 optimization runs) based on 400 scenarios, while Fig 3b shows the difference between the average selection patterns based on 400 and 2400 scenarios. More scenarios increased the number of surveys established at distant locations with low risk of infestation. The survey pattern at distant locations was stabilized for 2400 or more scenarios, which indicates that 2400 scenarios is sufficient to represent the majority of plausible outcomes from our model formulation, predicted by the stochastic spread model.
The number of spread scenarios also influenced key survey characteristics (Table 3). Increasing the number of scenarios from 400 to 3000 increased the proportion of sites surveyed from 39.4% to 71.3%, and also increased the proportion of the total budget devoted to survey from 45% to 59%. By including more long-distance dispersal events, a larger number of scenarios provides a more complete depiction of where ALB is likely to spread through time, and in turn, requires more sites to be surveyed.

General model behaviour
The proportional allocation of the program budget to surveys and tree removal is influenced by the eradication success constraint in Eq 8. The summation sign in Eq 8 makes it dependent on the size of the regulated area (J), the detection rate (γ) and the total number of host trees that are left after removal (N j −R js ). This indicates that the model's optimal solutions are likely to be influenced by the parameter values, so we explored the model's general behaviour in the parameter space {J, γ, β, θ js , N j }.
We first evaluated the optimal solutions for different values of β, which defines the proportion of the area of each selected site that is surveyed. We tested a range of possible values from 0 to 1; β = 1 implies that all host trees (i.e., the entire area) at the selected sites are surveyed, while β = 0 indicates that a decision-maker chooses preventive tree removal at the selected sites without performing surveys of the regulated area. Although preventive removal of trees without prior surveys may seem counterintuitive, this strategy was adopted for previous ALB eradication efforts in the GTA to reduce costs. Fig 4 depicts  When the regulated area J is small and the detection rate γ is high (as in Fig 4a), surveys of each site's entire area (i.e., β = 1) with subsequent tree removal yield the lowest expected cost. In this case, the surveys provide information about the presence of infested trees, which, in turn, helps reduce the cost of tree removal. When a larger proportion of a site is surveyed, there is less chance of infested trees going unsurveyed. Consequently, fewer susceptible host trees have to be removed to protect against the possibility that an undiscovered infestation will facilitate future spread.
While information gained from the surveys can thus help reduce the number of trees that must be removed, the capacity of the surveys to find infested trees depends on the detection rateγ. When γ is low, more infested trees are overlooked, and so more susceptible (i.e., possibly infested) trees must be removed to satisfy the eradication success constraint in Eq 8. With respect to our ALB example, when the detection rate is set to a comparatively low value (γ = 0.7), the lowest expected cost occurs at β = 1, meaning that the optimal solution in this case is to allocate all of the available budget to tree removal only (Fig 4b). An increase in the size of the regulated area J has similar impact: When J becomes very large, it is no longer optimal to survey the sites before tree removal, so instead the entire program budget should be allocated to preventive tree removal. In those conditions, information gained from the surveys provides only a marginal reduction of the total tree removal cost and therefore does not justify the cost of the surveys.
Switching between the "survey-and-remove" and "remove" policies The model behavior illustrated in Fig 4 suggests two alternative optimal management policies depending on the combination of the model parameter values. The "survey and remove" policy prescribes delimiting surveys in the regulated area J, with 100% coverage of the selected survey sites (i.e., β = 1) and subsequent removal of host trees based on the outcomes of the surveys. The "remove only" policy allocates the entire program budget to preventive tree removalaccording to the model optimal solution-without undertaking prior surveys. We further explored the model parameter combinations that cause the policy to switch. We depicted the switch between the "survey and remove" and "remove only" as a boundary curve in the dimensions of the size of the regulated area J and the pest detection rate γ (Fig 5). The area above the boundary curve corresponds to the parameter combinations for which the "survey and remove" policy is optimal, and the area below the boundary corresponds to the parameter combinations for which the "remove only" policy is preferred.
A lower pest detection rate restricts the optimality of the "survey and remove" strategy to small areas. However, the size of the regulated area for which "survey and remove" is optimal increases when the safety margin p or the eradication success threshold d is relaxed (Fig 5), i.e., when a decision-maker has lower aspirations about eradicating the pest from regulated area J. For instance, for the eradication success threshold d = 0.9, the boundary curve stabilizes around γ = 0.94, which implies that "survey and remove" policy is optimal for detection rates above that rate. For the empirical ALB detection rate in the GTA, γ = 0.7, the "remove only" policy is the only optimal choice. Similar to what occurs with lower values of the proportion of each site surveyed (β), lower detection rates translate to a higher number of infested trees that  are left undetected. When the number of infested trees that are left undetected is too high, surveys are not sufficiently able to reduce the number of trees that must be removed, since a large number of susceptible host trees have to be removed to prevent spread from undetected infestations.
Other model parameters, such as the expected proportion of trees that are infested (θ js ) and local host density (N j ), also influenced the optimal policy choice (S1 Fig). In general, lower infestation rates make the "survey and remove" policy attractive for larger regulated areas (J) or, alternatively, for lower detection rates (γ) (S1 Fig). The impact of decreasing the host density is similar: Lower host densities make the "survey and remove" policy attractive for larger areas because fewer trees will need to be removed in all cases. Lower host densities (or infestation rates) shift the policy switch curve down (S1 Fig). Note that this impact is only noticeable when the detection rate is high (i.e., γ > 0.9), in which case surveys can detect most of the infested trees. At lower detection rates (such as γ = 0.7 and below), the "remove only" policy is the only preferred choice.

Impacts of controlling the extreme program cost on the optimal policy
Recall that the objective function formulation in Eq 14 minimizes a combination of two objectives, the expected program cost and the CVaR α , which depicts the extreme cost in the right tail of the program cost distribution. Different decision-making preferences between minimizing the expected cost versus the extreme cost can be explored by changing the weighting factor F. When F = 0, the objective function minimizes the expected cost only. Increasing the F value above 0 places more emphasis on minimizing the cost in the right tail of the cost distribution. A decrease of the extreme cost causes the expected cost to increase. This penalty is expected: The only way to decrease the cost with respect to worst-case scenarios is to survey more sites. The total reduction of the extreme cost depends on the chosen values of the detection rate γ, eradication success threshold d and the safety margin p in the constraint Eq 8: a higher rate γ (or lower d and p values) lead to a greater reduction of the extreme cost. In a geographical context, the additional surveys in the solutions with CVaR α selected distant locations where the probability of spread is very low (Fig 6). Although these additional sites had considerably lower probabilities of ALB spread, they had higher host densities than the rest of the surveyed sites (Fig 7). These sites represent low-risk locations in terms of spread, but because of their abundant host, infestations could be severe and require costly eradication if they were to be invaded.
Controlling the extreme cost also influences the optimal management policy. In the current model setup, the only way to reduce the extreme cost is to survey more sites with respect to the worst-case scenarios. Thus, controlling the extreme program cost is only possible by increasing the number of surveys, i.e., by implementing the "survey and remove" policy. Nevertheless, the choice of the optimal policy also depends on how a decision-maker perceives the risk of incurring high program costs. Controlling the extreme cost adds a penalty to the expected cost, which makes the "survey and remove" policy less attractive when a decision-maker is risk-neutral and therefore aspires to minimize the expected cost (Fig 8, dotted lines). In this case, the boundary between the "survey and remove" and the "remove only" policies shifts towards higher detection rates, meaning that a risk-neutral decision-maker probably will not adopt the "survey and remove" policy unless the likelihood of detection is very high. In contrast, when a decision-maker perceives the worst-case costs as a reasonable proxy of the true program costs, the boundary shifts downward (Fig 8a and 8b, dashed lines), making the "survey and remove" policy attractive for lower detection rates and larger regulated areas. For the specific example of ALB in the GTA, with a detection rate γ = 0.7 and high expectations of eradication success (i.e., with an eradication success threshold d = 0.9 and safety margin p = 0.95; see Fig 8b), the "survey and remove" policy remains a preferred choice for regulated areas approximately 160 ha or less in size. However, for regulated areas larger than 500 ha, the "survey and remove" policy is optimal only if the detection rate γ is 0.8 or greater (Fig 8b, dashed line).

Discussion
Eradication of invasive species often requires costly investments. Agencies tasked with managing biological invasions face various budget and regulatory constraints that prevent them from undertaking full-scale eradication. When the need to eradicate an invasive species conflicts with a poor capacity to detect it, or when decision-makers face limited budgets, management policies can be developed that attempt to attain a desired level of eradication while meeting the budget and detection capacity constraints. The primary methodological contribution of our work is a resource allocation model that can help decision-makers working in terrestrial settings to develop cost-effective site survey and host removal strategies in areas under threat of invasion. Our model setup generally conforms to the current decision-making environment for managing ALB in the GTA, and depicts the management program as a two-step decisionmaking process. First, delimiting surveys are allocated over sites in a defined management area at the beginning of the survey season. The distribution of the pest within this management area, and the extent of damage, is uncertain at that point in time, and is modeled with a large set of stochastic scenarios, where each scenario depicts one possible outcome of the invasion. In the second stage at the end of the survey period, subsequent actions are applied to eradicate  infestations that are discovered at the surveyed sites. The model involves a set of constraints that focus regulatory decisions on two key parameters in particular: the desired level of eradication success and the margin of safety. These parameters are value judgments and define a decision-maker's aspirations regarding eradication and tolerance for program failure. Our model offers a practical approach to estimating the costs of these value judgments. The model formulation is generalizable and can be applied to other geographical regions and species of concern.
Our analyses revealed two alternative optimal management policies for ALB. The first choice prescribes delimiting surveys within a regulated area, with subsequent removal of host trees at infested sites based on the outcomes of the surveys, while the second choice prescribes preventive removal of infested and susceptible host trees at selected sites in the regulated area without prior delimiting surveys. Which policy is optimal depends not only on the value judgments mentioned above, but also on the combination of key assumptions, such as the capacity to detect the pest, the size of the regulated area, the expected infestation rates and local host densities. It turns out that the policy of preventive tree removal is optimal for the typical sizes of regulated areas that have previously been established to manage the ALB outbreak in the GTA (i.e., on the order of 46 km 2 ). Moreover, our results agree with past and current management practices for ALB in the GTA, where removal of infested and nearby susceptible host trees was initiated without delimiting surveys.
The optimal policy also depends on how decision-makers perceive the risk of eradication failure. In our ALB example, if a decision-maker strives to minimize the cost with respect to the worst-case invasion scenarios, delimiting surveys appear to be optimal for a larger size of regulated area than if the decision-maker perceives the expected program cost to be a reasonable representation of the true cost. The only way to reduce the cost with respect to the worst-case invasion scenarios is to survey more sites and more trees, which provides more opportunities to find infested trees and reduces the number of susceptible trees that must be removed to protect against overlooked infestations. However, the policy of minimizing the chance of a worst-case scenario imposes extra costs, thus making the expected program cost, on average, appreciably higher than the policy that strives only to minimize this expected cost.
The need for adequate control of the risk of high program costs Generally, decision-makers tasked with containing the spread of an invasive species desire to achieve successful eradication of the invader with minimum costs. Yet, there is always some risk that the costs of eradication could be very high. This aspect can be problematic for decision--makers, who may have low tolerance for incurring high costs and may agree that control of the risk of high program costs is important, even if this comes with a necessary increase in the overall costs. The level of risk that can be tolerated often depends on the type of organism, its capacity to spread and cause damage to economically viable hosts, as well as the objectives of the management program. We believe that adequate control of the risk of high program costs is an important part of a successful pest management program. A typical approach to manage risk in financial and national security applications is to estimate and control value at risk of high losses, costs or damages with a specified confidence level, such as 95% [53,79].
Adequate control of the risk of high program costs also depends on an analyst's ability to predict the outcomes of future invasions and management actions and to quantify other key factors that influence the species' ability to invade novel habitats. While the outcomes of ALB management actions in our GTA study (i.e., complete removal of host trees in an urban environment) were straightforward and well-understood, assessment of the outcomes of ALB invasion in natural ecosystems would require better understanding of factors that control the species' capacity to invade and establish a viable population, such as an ability to displace other species from their habitats or interspecific interactions, e.g., competition with native woodand bark-boring insects for host resources [6,80].
Our study was primarily focused on quantifying human-mediated dispersal of ALB. Human activities are known to be major contributors to rates of spread for many invasive insects [81]. In our case, the natural spread capacity of ALB by biological means is known to be poor [42,43] and the pattern of recently discovered ALB infestations in North America has been attributed largely to human-assisted movement [64,66], so we felt justified in our focus on tracking the human-mediated spread. Accounting for biological spread would be critical for an invasive pest species with strong flight capability that could cover long distances by its own means.
Spread rates may also be affected by other abiotic factors, such as climatic variation [82,83]. In our case, the study area was very small and the forecast horizon was relatively short (i.e., less than five years within GTA city limits), so therefore we did not account for site-to-site or temporal climatic variation and its potential impact on the rate of spread. Accounting for the impact of spatial and temporal climatic variation on the spread rates would be important for cases with long forecast time horizons when species are expected to spread over long distances.
Our study illustrates how the safety-rule approach with controlling the conditional value-at risk for program costs can be formulated as a mathematical programming problem. In our case, a linearized formulation of CVaR enabled solving the model with a linear programming solver for a large number of spatial cases and spread scenarios. Notably, minimizing the CVaR requires substantially greater computational effort than minimizing the expected cost. The CVaR metric is sensitive to the shape of the distribution tail above the chosen confidence level α, and the tail configuration affects the solution times. Nevertheless, our approach appears to be suitable for large-scale geographical applications and enables representation of the uncertainty of a pest's future spread in large regions via a large set of discrete scenarios.