Target the Source: Optimal Spatiotemporal Resource Allocation for Invasive Species Control

Determining how to cost‐effectively allocate resources for managing invasive species is a notoriously difficult problem. Invasive species control problems are always spatiotemporal, but much of the current theory about control strategies is either spatial or temporal. This article uses a deterministic spatiotemporal model of invasive species dynamics and identifies the optimal management strategy across a range of situations. The optimal solution points a principle of targeting the source of the population, which in many cases is the region of the landscape where the invader is most abundant. The analysis presented here is the first capable of solving for optimal strategies for invasive species over large and irregular environments. Thus, it is an important step forward for both the understanding of control strategies and the application to management scenarios.


Introduction
Invasive species cause massive environmental and economic damage (Duraiappah et al. 2005;Pimentel et al. 2005), and their management is a critical issue in conservation. Where possible, eradicating an invasive species in the early stages of an invasion is preferred (Courchamp et al. 2003), as allowing invasions to proceed unchecked can lead to disastrous consequences (Larson et al. 2011). However, determining cost-effective management strategies about how to allocate resources through time and space is particularly a difficult problem.
The spatiotemporal dynamics of a species is an integral part of devising a management strategy. Given the difficulty of solving spatiotemporal problems, many authors make simplifying assumptions that are useful, but leave some aspects of the problem unanswered. For example, focusing solely on the temporal dynamics (Hastings et al. 2006;Baxter et al. 2008) means that even though one can draw conclusions about how to allocate resources though time, there is no indication of how to allocate across space within the patch.
Despite the difficulties, numerous analyses have sought to find optimal management strategies to control spreading invaders. A prominent question in this literature is about allocating resources between main and outlier populations. Early work suggested that resources should be directed toward outlier populations (Moody & Mack 1988), and there has been much work that has focused on controlling these outlier populations to slow the spread of an invasion (Sharov & Liebhold 1998;Sharov et al. 2002). However, more recent work has shown that optimal policies should split resources between main and outlier populations (Whittle et al. 2007;Blackwood et al. 2010). Results that advocate allocating all resources to only main or only outlier populations likely stem from model assumptions, or optimization constraints, and are probably not a true reflection on reality, as in almost all situations it would be best to split resources between invaded regions (Hastings et al. 2006). Nevertheless, immediate removal of satellite populations is often regarded as the optimal strategy in both scientific publications and management guidelines (Murphy 2008;Ciosi et al. 2011;Gorchov et al. 2014). Even a good understanding of how to split resources between these different regions is only part of the picture, as control effort still needs to be allocated spatially within regions. There has been some work that allocates control resources across a region, but these often simplify grid cells to be invaded or empty (Epanchin-Niell & Wilen 2012), thus losing any density dependence of growth rate and control. This is important because newly invaded sites initially spread slower than established sites (Shigesada et al. 1995) and models that allow density to take any value predict this (see Appendix S1).
This article seeks to address these shortcomings by using a model that allows for density-dependent population dynamics and control, and multiple dispersal mechanisms, to solve for optimal management in a variety of scenarios. To provide useful and interesting results, this article takes two tracks. Firstly, providing optimal strategies in simple situations, seeking to understand the qualitative qualities of optimal management; and secondly, solving a complex, real-world example to show how the method is applied over broad regions and that the optimal solution agrees with simple cases.

Methods
Partial differential models are a flexible method for modeling the spread of a species. They can cope with coefficient values changing through space and time as well as various types of dispersal dynamics. Here, the canonical reaction-diffusion equation for species dispersal (Okubo & Levin 2001) is generalized to include advection and the effect of control efforts (Baker & Bode 2013). The invasive species abundance, n(x, t), at position x and time t is given by: The equation is defined over a region , whose boundary is labeled ∂ . The first term on the right-hand side is a diffusive term, which describes collective motion of randomly moving individuals throughout a landscape (Murray 2002). The coefficient D is the diffusivity, which governs how quickly the species disperses. The second term is advection, which is controlled by the coefficient v. This allows for dispersal to be biased in a certain direction, which is important for modeling dispersal that is influenced by external forces, such as currents or winds (Cousens et al. 2008). The third term is a logistic population growth term, where the coefficient r is the intrinsic growth rate of the population and K is the carrying capacity. The final term represents the species mortality due to a control effort E(x, t), which is restricted to nonnegative values. The effort is multiplied by the effectiveness, μ, and raised to the power q ∈ (0, 1). The parameter q represents diminishing marginal returns on increasing control effort; if q is small, then the marginal effectiveness of control reduces very quickly as control efforts are increased. This term is multiplied by the species abundance, n, meaning that a given level of control effort causes a proportional rate of decline of the population. The coefficients D, v, r, d, and μ can all vary in space and time (Tschinkel 1988;Hooten & Wikle 2008;Christensen et al. 2013).
Long-distance dispersal events can also be included in this model (Miller Neilan & Lenhart 2011), although the times and locations of any long dispersal events must be determined before the model is run. To do this, Equation (1) is solved forward in time until the time of the first long-distance dispersal. The abundance at that site is adjusted by a set amount. Then, the model is solved forward in time again until either another long-distance dispersal event occurs (and the process is repeated), or the simulation terminates. In this article, only long-distance dispersal events to sites that are unoccupied are considered. Alterations to the optimization would be required to consider the more general case of long-distance dispersal to an already occupied site (Ding 2007).
In general, the goal of invasive species management is usually either to minimize the environmental damage over time or to eradicate the species. In some circumstances, the goal may be to eradicate while also minimizing environmental damage throughout the eradication. Mathematically, all of these goals can be addressed using the objective (Baker & Bode, in press This balances the invasive species population damages against the costs of control through the coefficients ν and ω (see Baker & Bode, in press). The first term is the final invasive species population, multiplied by the weighting ν. The second term is an integral of space and time of the invasive species population, weighted by ω, and the control effort. This assumes a linear relationship between the invasive species population and environmental damage (Parker et al. 1999). Both terms in the integrand have been time discounted at rate δ.
The optimal solution is the effort allocation, E * (x, t), which minimizes the objective, Equation (2), subject to the system dynamics, Equation (1). This is identified by solving a system of two coupled partial differential equations (see Appendix S2): with initial and final conditions The boundary conditions chosen for Equation (3) determine the boundary conditions of Equation (4) as follows. If Equation (3) has Dirichlet boundary conditions, then (6) or if Equation (3) has Robin boundary conditions, then (7) wheren is the outward unit normal vector. The optimal control, E * (x, t), is given by Solving Equations (3), (4), and (8) gives the optimal solution (see Appendix S3).

Results
This model has many coefficients, and it is not possible to give the optimal solution for all plausible scenarios. This is further compounded because each coefficient can be a function of time and space and there are many possible initial states for the system. Hence, these results focus on situations where the abundance of the population is minimized at the terminal time ν > 0, ω = 0 and situations where the invasive species dispersal varies. However, extra results are included in Appendix S5 that explore other values of q, cases when ω > 0, cases where population dynamic coefficients vary in space, and the impact of altering the initial conditions. Figure 1 shows how the optimal solutions change when diffusivity varies. When the diffusivity is large, control effort roughly matches the density of the species. However, for the cases with low diffusivity, control effort is initially focused at the edges of the invasion. As time progresses, the effort allocation spreads out at a similar rate and shifts focus toward the center of the invasion. Mathematically, altering the value of the diffusivity is equivalent to altering the size of the invasion-halving the size of the equation is equivalent to quadrupling the diffusivity (see Appendix S4). Hence, this is a result that is about the ratio of the diffusivity of the species compared to the size of the invasion. For large invasions, or those with slow spread, it is better to focus the edges rather than the center. However, it is important to note that there is control effort over the entire invaded area. Figure 2 shows the optimal solution when there is an advection and when there is a long-distance dispersal event. In the advection case, the spread of the species is biased toward the left. The optimal solution mimics the solution with just diffusion, but it tracks the invasion as it moves to the left. When there is a long-distance dispersal event, the amount of control effort directed to the new outlier population is less than the main population. However, the ratio of control effort to density is greater in the outlier population, quickly reducing it to very low density.
These results can be explained intuitively by examining Equation (4)-the adjoint equation. The adjoint, λ, often known as the shadow value, describes what effect adjusting the invasive species density has on the objectivethat is, λ = d J dn (Lenhart & Workman 2007). Hence, it is most beneficial to remove individuals where λ is large (ignoring costs). While it is not possible to solve Equation (4) by inspection, because there is a complex relationship between n, λ, and E (Equation (8)), we can get some understanding of when the adjoint is small from the numerical results. The value of the adjoint tends to be relatively small in places with high invasive density or control effort, and its value increases at greater distances from these places. Hence, isolated areas with low density have higher values of the adjoint and therefore high control effort relative to the invasive species abundance. However, the control effort is still smaller in these regions compared to the main invasion, as the numerical results show, because of the smaller abundance (see Equation (8)).

Bogong High Plains
So far, we have only considered relatively simple onedimensional cases. However, these methods can be applied across real and complex landscapes, and in this section these methods are applied to the case of Hieracium aurantiacum in Australia. H. aurantiacum, originally from Figure 1 The optimal solution in one dimension with varying diffusivity. The top row has D = 12.5, the middle row has D = 1.25, and the bottom row has D = 0.125. The panels on the left show how the invasive species density changes through time, with time on the vertical axis, and dark red indicates regions where the invasive species is at carrying capacity, and black indicates that the species is at very low density. The right-hand panels show how control effort is distributed through space and time, where dark red indicates high control, and white indicates no control. This is broken up into contoured regions to more clearly show where the regions of highest control effort are and how they move through time. The boundary conditions are no flux at x = −25 and x = 25 and the time horizon is T = 2 (though the plot is cropped to T = 1. Other coefficients are q = 0.2, r = 0.25, K = 1, ω = 0, δ = 1, μ = 9, 000, and ν = 500. Europe, was unwittingly introduced to the Bogong High Planes (within the Victorian Alpine National Park) in the 1980s. Planted in Falls Creek village, it was discovered in other areas of the park in the late 1990s and an eradication program began. Although a spatially explicit optimal search effort allocation has been calculated for the region (Hauser & McCarthy 2009), there has been no optimization of how to allocate control effort for eradication.
The habitat suitability index (HSI) (Williams et al. 2008) provides useful information about how the growth rate and carrying capacity vary across the landscape, but the exact values of these coefficients are unknown. Values of the coefficients are chosen such that the model produces plausible dynamics-starting from H. aurantiacum spreading from a single location in Falls Creek village, the model must produce a distribution after 20 years that is consistent with the 1999 distribution. The HSI, which takes values between 0 and 1, alters the growth rate and carrying capacity, which are set to have maximum values of r = 2 and k = 1, respectively. H. aurantiacum spreads via both runners and seeds, meaning that it will spread in all directions, which is modeled with diffusion, as well as with the prevailing winds on the Bogong High Planes, which is in a south-easterly direction. Along with this slow drift of the invasion, there are long-distance dispersal events, where strong wind can carry seeds kilometers from Falls Creek. The diffusivity is set to 50 to model the expansion in all directions, and advection is set to v = (0.25-0.75) km/year so that the invasion drifts with Figure 2 Time slices showing the invasive species density and control effort for cases with advection and long-distance dispersal. Control effort is scaled relative to the maximum effort over the time window, and the invasive density is not visible initially for the long-distance dispersal case because control effort and the invasive density are almost equal.

Figure 3
The spread of H. aurantiacum through the Bogong High Planes as predicted by the model over a 20-year period. Brown areas show the presence of H. aurantiacum, with the darker shading indicating high density. The green shading is the habitat suitability, where dark green is high suitability and white is unsuitable. The blue region is a lake. H. aurantiacum spread starts from a single site in the north-west of the region. The long-distance dispersal events allow H. aurantiacum to establish colonies on the south-east of the lake. the wind. The advection term is smaller than diffusion so that the model still allows local spread in all directions, as the main effect of the wind is modeled with long distance dispersal events. These occur randomly each year from within Falls Creek. They result in a small increase in the H. aurantiacum density at a random location (the increase at the target location is 15% of the density at the source, multiplied by the HSI of the target location), biased by the predominant wind direction and a distance dispersal kernel (Williams et al. 2008). Figure 3 shows the expansion of H. aurantiacum over 20 years, using this model.
The optimal control strategy is shown over a 5-year period ( Figure 4). Due to the fairly slow spread rate of H. aurantiacum, there is only a small change in the spatial region that receives control, although the intensity of control effort increases greatly over the period. The main populations, north-west of the lake, receive the greatest control allocation, similar to Figure 2. The increasing intensity of control through time reflects the ease of reducing the abundance when from high density as well as the effect of diminishing marginal returns on control effort, and this aspect of the solution is consistent with the previous work on optimizing the temporal allocation of resources (Baker & Bode, in press).

Control principles
These results show that it is important to target the source of invasions. If the species is capable of spreading quickly, relative to the invasion size (e.g., if seeds from a plant at the center of the invasion could spread to uninvaded areas), then control effort should be proportional to invasive density. However, if this is not the case, then control effort can be focused more toward the edges of the invasion, rather than the center. This is because the center can be reinvaded quickly as it is surrounded by the invasion. Hence, propagule pressure comes from every direction. This is in contrast to the edge of the invasion, where propagule pressure only comes from one side. Isolated outlier regions should generally receive less attention than a main invasion, due to smaller abundances and limited ability to spread. However, if there are two comparable sites, one in the main invasion and the other in an outlying area, which have equal density, then more effort should be directed toward the outlying site.
The method of controlling outlier populations permeates weed management thinking. The Bradley method (Bradley 1988) suggests to manage invasion starting from the outlier populations, and, supported with modeling by (Moody & Mack 1988), this method has had a big impact on weed management practice (Bossard et al. 2000), and these methods are implemented worldwide (Credit Valley Conservation 2012;Tzankova & Concilio 2014). For example, Rubus fruticosus agg. (blackberry) causes $100 million damage per year in Australia (Page & Lacey 2006), and management guidelines advocate treating small outlying populations initially (NSW Department of Primary Industries Weed Management Unit 2009). As the results in this article show, managers must shift their focus to main populations.

Figure 4
The optimal effort allocation to control H. aurantiacum over 5 years, starting from the distribution in Figure 3. The top row shows the H. aurantiacum density over the 5 years (from year 20 to 25), and the bottom row shows the distribution of management effort over the same time period. Brown areas show the presence of H. aurantiacum, with the darker shading indicating high density, and pink regions indicate the intensity of control effort. The green shading is the habitat suitability, where dark green is high suitability and white is unsuitable. The blue region is a lake. The final time weighting is ν = 1, 000 and diminishing returns is q = 0.2.

Discussion
Solving for the optimal spatiotemporal resource allocation for invasive species, eradication reveals important control principles. The key principle that has emerged from this analysis is that control efforts should target regions that are sources of propagules. These results can be logically explained, but they contrast with previous work, particularly about how to manage outlier populations. Results from population biology have supported the idea that satellite invasions should be removed first, while metapopulation approaches have suggested that source populations should be targeted (Moody & Mack 1988;Epanchin-Niell & Hastings 2010;Chades et al. 2011). The results in this article support the results that originate from metapopulation analyses. The reason that previous analyses using a population biology approach have favored the control of satellite invasions is that the models generally have not accounted for the establishment phase of satellite invasions and hence overstated their importance. In fact, "target the source" would also say to remove satellite invasions quickly in those models.
There are many ways to model the control of invasive species. A popular alternate method has been metapopulation modeling and dynamic programming approaches (Chades et al. 2011;Epanchin-Niell & Wilen 2012;Walker et al. 2015). These analyses include multiple patches and stochastic growth and spread in the optimization. However, they are usually unable to incorporate diminishing returns on control effort, density-dependent growth, and are restricted to a relatively small number of patches. Walker et al. (2015) include diminishing returns and density dependence, but this further limits the number of patches-they only solve a three-patch problem due to computational difficulties. In contrast, the methods in this article allow for a fairly general characterization of invasive species dynamics and can be used to solve large problems; the H. aurantiacum case has 3,500 cellswell beyond the current capabilities of stochastic dynamic programming.
This article does not take into account detectability; rather, it produces the optimal way to manage a species with a known whereabouts. In invasive species management projects, there would often be a shared budget to cover the costs of both detection and control. Including surveillance costs would alter the optimal solution, as detecting long-distance dispersal events is notoriously difficult. Future research should seek to find the optimal resource allocation, split between surveillance and control, over time and space for managing an invasive species. There has been some progress toward this ), but to date there is no fully dynamic optimal solution.