Optimal control of an invasive species using a reaction-diffusion model and linear programming

,


INTRODUCTION
Invasive species management is important worldwide, as invasive species are one of the top threats to global biodiversity (Parker et al. 1999, Mack et al. 2000, Pimentel et al. 2005, McCleery et al. 2015).Management of an invasive species is a complicated task for several reasons.First, if the species is newly introduced, it is generally difficult to characterize its population dynamics in its new habitat.Consequently, it is difficult to predict the possible impact on native species or to predict its current and future spatial distribution.Second, management practices may not be efficient due to factors such as low detectability or underdeveloped removal methods, both caused by a lack of training and experience facing the new invasion.This uncertainty ultimately creates challenging decisions for natural resource managers: Given actual funding and other constraints, how should resources be invested without knowing the potential ecological impact, future spatial distribution of the species, or efficacy of that management?
Researchers have provided several insights for invasive species management, and not surprisingly it appears that the rule is the earlier the better.Indeed, studies have shown that investment in prevention is more efficient than investment in post-detection control (Leung et al. 2002(Leung et al. , 2005)).If an invasive species is newly detected, it is then important to invest in search effort to control the population as soon as possible.Mehta et al. (2007) proposed computing the optimal search effort as a function of the species' characteristics and the cost of proposed management.An effective search effort can enable detection of the species sufficiently early, and therefore, the population might be small enough that eradication might still be a viable option (Sharov andLiebhold 1998, Olson andRoy 2002).However, in most cases the management of an invasive species starts only because the damages caused by the species are visible.To determine whether containment, eradication, or non-action is the best management action, one can use a decision analysis framework as Moore et al. (2011) did for Acacia paradoxa in South Africa.
The optimal management decision can also be constrained by other factors.For example, some political jurisdictions require that action be taken to mitigate the impacts of an invader (Simberloff et al. 1997, Carlton 2003) and non-action is not an option: The invasive population has to be reduced as much as possible with the available resources.In this case, reducing the invader's population within the management unit is an obligation, and resources must be deployed over several years, even in the face of uncertainty concerning, for example, the population growth rate, dispersal behavior, range, diet, predators, and current spatial distribution of the species.
As with most environmental management problems, the problem of Dynamic (i.e., over several years) Spatial Allocation of Resources (DSAR) should be approached using an adaptive management framework (Chad es et al. 2015).In this case, the optimal strategy is expected to perform well on average over all possible invasion scenarios, and the species' population is monitored every year in order to reduce the uncertainty in key parameters.But in the early stage of an invasion, the available resources, the coordination between management agencies, and the research capacity might not be sufficient to implement this approach.From a modeling point of view, adaptive management generally relies on the framework of stochastic optimal control or Markov decision processes (Marescot et al. 2013), two mathematical frameworks allowing to compute an optimal strategy.Unfortunately, these frameworks generally limit the size of the decision problem, that is, the number of management units, actions, and the time frame.For large problems, a different framework must be used to determine the best management option.Hauser and McCarthy (2009) developed a framework to compute the optimal spatial allocation of effort as a function of the searching efficacy and the species' probability of presence.In their work, the spatial dynamics of the species is not explicit but is included in the cost of having the species undetected at a site.In Blackwood et al. (2010), the spatial dynamics of the species and its density are made explicit.Their formulation is particularly convenient because it allows explicit computation of the optimal spatial resource allocation over time without requiring intensive computational resources.More recently, Baker (2016) proposed an analytic formulation of the optimal control effort over time, which exhibits an important principle: It is optimal to focus the control effort on the source of ❖ www.esajournals.org 2 October 2017 ❖ Volume 8(10) ❖ Article e01979 the infestation, where the invader is generally most abundant.This work is particularly interesting in its generality; the spatial dynamics of the species is described using a reaction-diffusion (RD) model, with density-dependent growth rate and advection, which allows consideration of different dispersal speeds depending on the dispersal directions.
All of these approaches attempt to minimize an objective function based on the sum of the ecological and economic costs resulting from the invasive species.In other words, the optimal management option minimizes the management cost, as well as the ecological and societal damage, over the long term.Using such objective functions allows the identification of the optimal strategy that represents a trade-off between management cost and species damage.Without considering ecological and societal costs, the optimal strategy is simple to compute: Do nothing.Without considering management cost, it is again simple to compute: Remove all the individuals from the landscape.Both strategies are generally unrealistic, which underlines the importance of defining each cost properly.However, although management costs remain relatively easy to estimate from data, it is particularly difficult to quantify the ecological and societal damages (but see Kaiser and Burnett 2010 for an example).A common practice to address this difficulty is to instead define management constraints, for example, defining a maximum management budget or defining a maximum invasive species population level.Adding these further constraints in a convenient way also requires another framework.
Epanchin-Niell and Wilen (2012) proposed a practical solution, in which the management objective and constraints are defined with linear equations.Then, the DSAR problem can be defined as a linear programming optimization problem (Hof and Bevers 2002), which allows time-efficient computation of the optimal control strategy for large problems.The work of Epanchin-Niell and Wilen (2012) allows us to draw general conclusions about optimal spatial allocation of effort in stylized cases, but they only considered the presence/absence of the species and a simple dispersal model where adjacent cells get infected in the next time period.Here, we propose a broader look at the possibilities that the linear programming framework offers, and illustrate its use on the management of the invasive Argentine black and white tegus (Salvator merianae, hereafter "tegus") in southern Florida.We considered the problem of spatial allocation of resources over time, when the invasive population has to be minimized under several constraints and uncertainty.
We first paired linear programming with a discrete RD model (Levin 1974).Discrete RD models are relatively general (Holmes et al. 1994), and we proposed a simulation approach to parameterize the dispersal coefficients, based on any dispersal kernel.The pattern of invasion proposed by RD models is simple: In an infinite and suitable landscape, invasion is described by a growing circle centered on the point of first introduction.Discrete RD models are particularly useful as they define the DSAR with a very general framework (Hof and Bevers 2002), and the optimal spatial strategies can be computed time-efficiently with common linear programming solvers.We use this framework to study several questions: Where to allocate traps in order to minimize the invader's population?What are the minimal required resources to implement containment?How much and where to allocate resources to protect a specific area from invasion?Because there is uncertainty about the tegus' population dynamics, we considered 25 different spatial dynamic models.These three questions were then studied for each of the 25 models, and the results were generalized depending on the risk attitude of the decision-maker.

Problem definition
Tegus are large lizards native to southeastern Brazil, Paraguay, Uruguay, and northern Argentina.Individuals were likely introduced via the pet trade around 2002 in Florida City, FL (Engeman et al. 2011).Managers are concerned that these habitat and dietary generalists will affect native species, particularly those of conservation concern such as the American crocodile (Crocodylus acutus, Mazzotti et al. 2015), federally recognized as threatened.The Everglades National Park (ENP) and the Turkey Point Power Plant (TPPP; see Fig. 1 for the locations) are of primary concern as far as they host the American crocodile and several other native species.In South Florida, tegus are typically found in disturbed ❖ www.esajournals.orghabitats, including urban areas, agricultural fields, and roads and levees, that extend into natural areas.They are active during the warmest period of the year and overwinter the coldest months (i.e., October to February/March) by brumating in burrows.More details on the invasion of tegus in Florida can be found in Klug et al. (2015) and Johnson et al. (2017).

Discrete reaction-diffusion model
We propose to model tegu population dynamics using a difference equation model, which allows a description of the population over discrete space and time (Levin 1974, Holmes et al. 1994).More precisely, we propose to use the following model, obtained by using a finite difference scheme of the classic RD model (Allen 1987): where n(i, t) is the population density at site i 2 1; . . .; s f gat time t and n(j, t À 1) is the density at site j 2 1; . . .; s f gand time t À 1. s is the total number of sites in the landscape, f nðj; t À 1Þ ð Þis the per-capita rate of recruitment, and g ji is called the dispersal coefficient, the proportion of the population in site j moving to site i.We first propose to use a density-dependent per-capita rate of growth, for any site j 2 1; . . .; s f g : The studied areas are colored in green and beige.Green areas are the roads and levees with a 1-km buffer on each side.Beige areas are the agricultural and urban areas.We assumed that these areas represent the only suitable habitats for tegus (Salvator merianae) such that they are not present anywhere else.On the contrary, control is available only in the green areas and the Turkey Point Power Plant (TPPP).The initial point of introduction is represented with a star.Red lines show the boundaries of the Everglades National Park as well as those of the TPPP.Finally, the black circles are the 25, 50, 75, and 99% envelope, centered on the first point of introduction.For example, the 25% envelope is represented with the smallest disk and contained at least 25% of the population for each of the model.(b) Spatial distribution of the EDDMaps data (https://www.eddmaps.org/)used to estimate the dispersal coefficient of the reaction-diffusion model.(2) where e is the population's intrinsic rate of growth and K j is the carrying capacity of site j.
Þ nðj; t À 1Þ thus describes the natural birth and death phenomenon in site j and it can be viewed as the population in site j, just after the birth and death pulse, but just before any dispersion occurs.
g ji is the post birth and death population that leaves site j to reach site i.Eq. 1 simply indicates that the population at any site i is the sum of all individuals that are immigrating (i.e., all the terms in the sum for j ¼ 1; . ..; s; j 6 ¼ i), plus the individuals staying in the site (i.e., only the term j = i).Other formulations of discrete RD model can be found in Hof and Bevers (2002).
Finally, we denote the population vector of length s, ...; s , where each component is computed according to Eq. 1 and represent site's density at time t.

Tegu spatial distribution
Modeling tegu spatial distribution using a discrete RD model requires an estimate of the intrinsic rate of growth e, the initial population value N 0 , the carrying capacity K i ð Þ i¼1; ...; s , and the dispersal coefficient g ji À Á j;i¼1;...; s .Intrinsic rate of growth.-For the rate of growth, we rely on Johnson et al. (2017), where a probability distribution of e was estimated from expert elicitation.Johnson et al. (2017) estimated that there is a probability of 0.3 that the population's growth rate e is negative.In this case, the population is decreasing by itself and will go extinct after several years even if the population is not controlled.For negative values of the growth rate, other questions might be investigated, such as minimizing the time to extinction, but we choose to focus only on the case where the growth rate is positive.Details on the estimation of the growth rate for a RD model can be found in Bonneau et al. (2016).
Initial population.-Wealso relied on Johnson et al. (2017) for the initial state of the population.We used the mean of the estimated size of the population in 2009, initializing the population with 770 individuals.We also relied on experts to determine the location of the first point of introduction (25°25 0 44.2884″N, 80°28 0 51.64″ W, see Fig. 1).The 770 initial individuals were spread around the first point of introduction such that the carrying capacity is respected.
Carrying capacity.-For carrying capacity, little information is available but it has been observed that tegus are mostly found in disturbed habitats, such as roads and levees on public land.We used a carrying capacity of 0 in marsh and mangrove, which represents undisturbed wetland habitat.We used the best estimate of the experts, that is, K = 240 per km 2 , for all the disturbed/suitable habitats, including marsh habitats within a 1-km buffer from each side of roads and levees.We considered s = 3465 sites with suitable habitat for tegus and each site is a square of 500 m length.See Fig. 1 for the location of the suitable habitats.
Dispersal coefficient.-Thedispersal coefficient g ji represents the proportion of the population in site j that emigrates to site i.We assume that the dispersal kernel follows a normal distribution with mean zero and variance r 2 .Little is known about the dispersal of tegus, so we used a common symmetric distribution, although we are lacking data to support this choice.However, this choice is particularly convenient, as it allows using EDDMaps presence-only observations available at https:// www.eddmaps.org/to estimate the variance r.As for the growth rate e, we propose to fit a probability distribution (a normal distribution here) on r in order to account for uncertainty.We then used a simulation approach to estimate the dispersal coefficients ðg ji Þ j;i¼1; ...; s .The basic idea of our approach consists of first using the dispersal kernel to simulate dispersion of a large toy population and then to estimate the proportion of individuals that emigrates in each site of this large toy landscape.These proportions are finally used to compute g ji , depending on the distance between site j and site i.This approach is general and can be used with any dispersal kernel.For clarity, details are proposed in Appendix S1.
In theory, it is possible to simulate long-distance dispersal, but one possible consequence is all g ji > 0, even for very distant pair of sites (j, i).If long-distance dispersal needs to be included, we recommend considering instead some scenarios (see the following section) where long-distance dispersal events have happened.Finally, models with a different dispersal kernel can also be considered to reflect structural uncertainty on ❖ www.esajournals.orgdispersal.For example, an asymmetric dispersal kernel can also be considered in the set of possible models.

Parameter uncertainty
The growth rate e and the standard deviation r of the dispersal kernel are described by continuous distributions.In order to propose a finite set of possible values for these parameters, we used a Gaussian quadrature (Miller and Rice 1983).Gaussian quadrature can be used to represent a continuous probability distribution on finite support points, while preserving the moments of the distribution.We approximated each continuous distribution of e and r using a discrete probability distribution with five support points.We thus considered five possible values of the growth rate and, for each of them, we computed five possible values for r.In total, we considered 25 different pairs of growth rate and standard deviation of the normal distribution kernel.The models are denoted M 1 , . .., M 25 , associated with the coupled parameters e 1 ; r 1 ð Þ, . . ., e 25 ; r 25 ð Þ and model probabilities p 1 , . . ., p 25 .The possible values of the growth rate and the variance of the Gaussian dispersal kernel are provided in Fig. 2. We sorted the models by ascending probability such that: Note that it is possible to add other types of uncertainty by adding more models.For example, one may wish to analyze the effect of a secondary point of introduction.In this case, models incorporating this event can be added to the set of possible models, and the optimal spatial control can be computed to be robust to this extra source of uncertainty.
In the following, we show how this spatial model can be used to answer practical questions.

Tegus spatial distribution
In order to minimize the tegu population without budget constraints, it is optimal to control the area of highest density (Baker 2016).However, in our case the spatial distribution of tegus is unknown and we can only compute a prediction of the spatial distribution for each model.We propose to summarize the prediction of all models by computing the minimal disk or envelope, centered on the first point of introduction which contains at least 25%, 50%, 75%, and 99% of the tegu population.For example, the 25% envelope contains 25% of the tegus population with probability 1.It is difficult to assess the quality of these predictions as little is currently known on the spatial distribution of tegus.However, we propose to compare these envelopes to the presence-only observations available at https://www.eddmaps.org/.

Optimal spatial control
Let us define the control vector Y t , where each component Y t (i) represents the amount of control effort used in site i.And let us denote a i the catch per unit effort (CPUE) of the control method in site i.When control is implemented, Eq. 1 becomes: where a i 9 Y t (i) represents the number of individuals that have been removed from site i.For tegus, control is implemented using traps and Y t (i) is then the number of days a trap has been placed in site i.The CPUE was estimated from the current control program as a i = a = 0.035, for all i = 1, . . ., s. Trapping is only relevant when tegus are active and thus, the amount of control Y t has to be bounded.The active period of tegus is approximately 238 d in South Florida, and we assume here that a maximum of five traps can be used in each site, which represents traps every 100 m on a linear feature.Thus, we have a maximum of 0 ≤ Y t (i) ≤ 5 9 238 = 1,190 trapping days per site i.It is also necessary to determine the per-unit effort management cost.The active season trapping cost is estimated to be US$1,000, and the per-unit management cost is then c u ¼ ðUS$1; 000=238Þ.It is likely that a nonhomogeneous management cost would be more realistic.For example, the distance between the trap position and the management office, as well as the type of roads between these two positions, can significantly impact the cost.Theoretically, non-homogeneous cost can be modeled by defining c u as a vector of size s, where each component is the per-unit management cost for a given site.Unfortunately, we lacked relevant data to parameterize such a model.Finally, we assumed that traps can be used on roads and levees on public lands, but not on private land.We defined C as the set of sites where control can be implemented.It is composed of roads and levees, with a 1-km buffer on each side of them.See Fig. 1 for the area available for control in green.

Minimize tegu abundance with a fixed number of traps
General method.-We first start with the practical situation where a fixed number of traps are available and one wants to compute their optimal locations in order to minimize the number of tegus after control.Note that this also corresponds to the problem of maximizing the number of trapped animals.In the following, we show how this problem can be defined as a linear integer programming (LIP) problem.
We first suppose that traps are not moved during the season, and we let p be the vector of the number of traps that are used in the sites.Our problem consists in finding the optimal trapping strategy, p*, or in other words, the optimal allocation of trapping effort minimizing the remaining number of tegus.The number of removed tegus depends on the strategy, but also the number of individuals that are present in the site, and thus on the population model.Let us define R t M m p ð Þ as the vector of the number of removed tegus in the sites at time t when model M m is true and strategy p is used.The number of remaining tegus after control at time t when model M m is true, denoted , is then computed as follows: where Formally, the objective of our problem is to find: To define a realistic problem, all the variables should be bounded and we now define the constraints of the LIP problem.
It is easy to show that the number of removed tegus should not exceed p i Â a Â 238 ð Þ , which is the maximum number of removed tegus in any site i when p i traps are used in this site during the entire active period.On the other hand, it should as well not exceed n M m ði; tÞ, that is, the number of tegus that are present in the site at time t when model M m is true.The number of removed tegus in site i, R t M m p i ð Þ, should then verify: for all i 2 f1; . . .; sg: The number of traps should also be bounded and we defined n Trap Max as the maximal number of traps that are used such that: for all i 6 2 C; p i ¼ 0: (9) Eq. 6 through Eq. 9 represents the constraints of the LIP problem and together with Eq. 5, they make sure that the optimal trapping strategy will respect the characteristic of our problem.Our LIP problem can be summarized as follows: 8i2f1;ÁÁÁ;sg; 0 p i 5; 8i6 2 C; p i ¼0: Note that for this problem, tegu population size has to be minimized for only one time step, and thus, t = t 1 .
Coping with uncertainty.-Problem(10) can be solved for each possible model M 1 , . . ., M 25 , providing 25 different strategies p Ã M 1 ; . . .p Ã M 25 .When choosing among them, we must consider the efficiency of the strategy under every model, all representing a possible invasion scenario.For example, strategy p Ã M 2 is optimal when M m 2 is the underlying model, but it can exhibit weak performance if the underlying model is not model M m 2 .Then, a decision-maker has several choices based on their risk attitude: (1) A risk-averse attitude would choose the strategy that performs best in the worst-case scenario, (2) a risk-neutral attitude would choose the strategy that performs best, on average, over all models, or (3) a risk-seeking attitude would choose a strategy that performs best in the best-case scenario.
Among p Ã M 1 ; . ..; p Ã M 25 , the MiniMax strategy p MiniMax is a robust strategy that minimizes the maximal possible remaining number of tegus, over all possible models: In other words, the value of each strategy p Ã M m is computed in case any of the possible models is the remaining number of tegus when strategy is the worst possible performance of the strategy p Ã M m over all possible models.As a consequence, p MiniMax is the strategy that minimizes the remaining number of tegus in the worst-case scenario.This is the strategy that should be selected by a risk-averse decision-maker.
Second, the MiniMean strategy, p MiniMean , is the strategy that achieves the minimal expected number of remaining tegus: In other words, p MiniMean is the strategy that minimizes the remaining number of tegus, on average, over all possible models.This is the strategy that should be selected by a risk-neutral decision-maker.
Third, the MiniMin strategy, p MiniMin , is the strategy that minimizes the minimal possible remaining number of tegus over all possible models: In other words, p MiniMin is the strategy that minimizes the remaining number of tegus in the bestcase scenario.This is the opposite of p MiniMax and it should be selected by a risk-seeking decisionmaker.Finally, we also computed the optimal control effort using an average model, M av .In this case, the growth rate and standard deviation of the dispersal kernel were fixed to their average value.For the average model, the standard deviation of the dispersal kernel is r = 3 and the growth rate e = 0.28.Note that to compute the average growth rate, we used both negative and positive values from Johnson et al. (2017) in order to use an average model that includes all possibilities.

Robust containment strategy
The term containment is used to describe a strategy that aims to confine the population to a small part of the landscape.In practice, it is difficult to determine the size of the containment area that is achievable and the budget that it might require.Containment is generally only achievable in the early stage of the infestation and when ❖ www.esajournals.org8 October 2017 ❖ Volume 8(10) ❖ Article e01979 abundance is relatively low.We assume that containment is still achievable for tegus and that containment will help protect incursion by tegus into the Everglades National Park (ENP) and the Turkey Point Power Plant (TPPP).General problem.-Weseek to implement a diskshaped containment area of radius d, centered on the point of first introduction.The objective is to minimize the control cost to attain a containment area over a time frame of four years.We keep this time frame relatively short for computational reasons but also because budgets are generally planned for a short time horizon and the credibility of different tegu population models is likely to be updated in subsequent years.
Solving such a dynamic optimization problem requires us to reject logistic growth of the population and consider exponential growth instead: Using exponential growth allows us to express the density in one site as a linear combination of the density in other sites, which is required to use the LIP framework.On the contrary, quadratic terms appear when logistic growth is used and the population dynamic is computed over a time frame >1.On the previous example, using logistic growth was possible because the population dynamics were not explicitly accounted for, but only the population at a given time step t.One consequence of the exponential growth is that the tegu population size and the control cost will be overestimated.For models where the predicted density remains low, the difference with a logistic growth will be minimal.In our case study, it may not be the case that tegu population density is low everywhere.Interpretation of results from the following sections for such high densities must be done with caution, until such time as exponential growth is verifiable or methods allowing logistic population growth are developed.
The linear programming problem is as follows.If model M m is assumed to be true, minimizing the control effort to design a containment area of radius d remains finding the control vector p c t M m ; . . .; p c tþ4 M m that minimize the control cost P tþ4 t 0 ¼t P n i¼1 c u Â p c t 0 M m .Once again, for the problem to be fully described, several constraints have to be defined: (1) The density of tegus is computed according to Eq. 3 with exponential growth, (2) n M m i; t À 1 ð Þ is initialized according to model M m (3) the density of tegus is 0 anywhere outside the containment area during the four years, (4) the number of trapping days per site does not exceed the maximal number of trapping day using five traps, that is, p c t 0 M m 238 Â 5 for all time t 0 , and (5) the number of removed tegus does not exceed the number of tegus in the site, that is, a Â p c t 0 M m n M m ðt 0 ; iÞ for all time t 0 .We suppose that control starts in 2016 and thus t = 2016.Note that in this case, we supposed that traps can be moved during the season, and thus, p c t 0 M m can be lower than 238 d.We also assumed that trapping is available everywhere in the landscape; otherwise containment is impossible.Note that constraint (3) is drastic, and one can be more flexible in authorizing the tegu's density to be below a given threshold outside of the containment area instead of fixed to zeros.
Coping with uncertainty.-In this case, it is hard to follow the same method as in the previous example with a fixed budget.Indeed, an optimal strategy p cÃ M m is designed to guarantee that no tegus will be outside of the containment area for a minimal cost, only if model M m is true.But for other models, some tegus might be present outside of containment.It is then hazardous to compare two strategies only in terms of cost, as far as one strategy can be cheap, but only allows containment for few models.Instead, we compared the strategies in terms of the number of tegus that are still present outside of the containment area when all the models are accounted for.We denote the number of remaining tegus outside of the containment area when using strategy p cÃ m 2 but model M m 1 is the true model.Note that we have five traps per site allow removal of all individuals.Indeed, we recall that a maximum of five traps can be used in any site and when the density is too high, more than five traps would be needed to remove all the individuals.Similar to the problem with a fixed budget, we considered a MiniMax strategy, p c MiniMax , which minimizes the number of tegus outside of the containment area in the worst-case scenario: ❖ www.esajournals.orgIn addition to these three strategies, we considered a Max strategy p c Max , which on any time step and any site, use the largest control effort among all the possible strategies: where p cÃ t 0 M m ðiÞ is the amount of effort in site i during year t 0 when strategy p cÃ M m is used.Finally, we also computed the optimal control effort for the containment of tegus using the average model M av , as described before.
Note that to compute the remaining number of tegus outside of the containment area, V c M m : ð Þ, we computed the population dynamics using Eq. 3 with exponential growth; thus, the carrying density K = 240 is not necessarily guaranteed.But to compute the associate cost of the different strategy, we only consider a maximal density of K = 240, preventing comparison of strategies with an irrelevantly high number of tegus and cost.

Protection of Everglades National Park and the Turkey Point Power Plant
We consider a final case, where the problem is to compute the minimal effort that is needed to protect specific parts of the landscape (the ENP and the TPPP in our case), for all scenarios.Unlike the previous case study, containment is not applied on all the landscape but just for specific parts.In addition, we consider that having tegus is these zones is not an option, in any of the possible scenarios.Then, we ask what is the minimal CPUE value, and where do we apply control in order to protect these two areas in the worstcase scenario?We define the worst-case scenario as the model that predicts the highest density of tegus in the study area in 2015; it is model M 2 in our case and we let M Worst = M 2 .
The linear programming problem consists of minimizing the management cost, such that there are no tegus in the ENP and the TPPP for the next four years.The constraints on the tegu population dynamics and on the maximal number of tegus are unchanged.

RESULTS
Where are the tegus?
The boundary of each envelope is provided in Fig. 1.The radius of each envelope is 6.8, 9.7, 12.1, and 17.4 km.In other words, 25% of the tegus population are predicted to be found within a distance of 6.8 km of (x 0 , y 0 ) with probability 1.The average density in each envelope is 136, 121, 105, and 76 tegus per km 2 .As expected under the assumptions of a RD model, the density of tegus decreases with increasing distance from the point of first introduction.Note that a simple "rule of thumb" control strategy could be used by randomly selecting cells within the 25% envelope, a core area where the tegus density is highest.
The predicted density for the least likely model M 1 , the average model M av , and the most likely model M 25 are available in Fig. 3.The predicted density for all models are available in Appendix S1: Fig. S2.The highest infestation levels are predicted by the model with the lowest likelihood, and the density of tegus is a decreasing function of the distance from the point of first introduction.Nearly 86% of the observations from EDDMaps are located in the 25% envelope, 95% of the observations in the 50% envelope, 97% of the observations in the 75% envelope, and 99.9% of the observations in the 99% envelope.
Where to control with a fixed number of traps?
The MiniMax, MiniMean, MiniMin strategies, and the optimal strategy for the average model are displayed in Appendix S1: Fig. S3.For all strategies, all of the traps are located within the 25% envelope.Note that all the 25 strategies ❖ www.esajournals.org10 October 2017 ❖ Volume 8(10) ❖ Article e01979 minimized the maximal values over all models; that is, there are 25 other possible r MiniMax strategies.As explained earlier, the r MiniMax strategy is designed to minimize the worst-case scenario and it occurs for the first possible model, with the lowest likelihood.In this case, the CPUE and/or the number of traps are too low to really impact the tegu population.Any strategy that randomly selects locations in the 50% envelope is expected to perform the same (only in the worst-case scenarios).Indeed, a maximum of a 9 5 9 238 ≃ 42 tegus can be trapped using five traps during 238 d in the same cell.And in the worst-case scenario, every cell in the 50% envelope has a minimum of 56 individuals.For a maximum of 200 traps, the risk-neutral and risk-seeking strategies r MiniMean and r MiniMin consist of spreading the traps in the 25% envelope, as close as possible to the point of first introduction.With 300 traps, the risk-averse (r MiniMax ) and risk-neutral (r MiniMean ) strategies use up to two traps around the point of first introduction.This is because the density of tegus is predicted to be high in this area for most of the models.In contrast, the risk-seeking strategy (r MiniMin ) focuses on situations where it is possible to completely eradicate tegus from the controlled cells with only one trap.As a consequence, this strategy consists of spreading the traps inside the 25% envelope as much as possible, which helps reduce the risk of using more traps than needed in any given cell.
When using an average model, there is no uncertainty around the tegu density, and it is therefore easier to determine if using more traps in a given cell is wasteful or not.Then, two traps are used in six (or 16) cells around the point of first introduction when 200 (or 300) traps are available.The remaining traps are then spread in the 25% envelope.

Is containment a feasible option?
We estimated that the probability that a containment area can be established is equal to 0.66 when the width of the containment area is ≤4 km and is equal to 0.93 for a containment width from 5 to 8 km.For the larger containment width ≥9 km, containment is achievable with a probability of 0.99.Fig. 4a shows the distribution of the expected number of tegus outside of the containment area.The optimal containment strategy for a given model is not supposed to achieve containment for any other models, and thus, the number of individuals outside of the containment area might not be zero for every model in the model set.It is apparent from Fig. 4a that more tegus are expected to be outside the containment area for the smallest containment widths, as it requires removal of tegus from sites with high predicted density (close to the point of first introduction).Fig. 4b shows the distribution of the cumulative containment cost, for all the containment strategies r cÃ M m , for models M m where containment is feasible.A wider radius d allows containment for more models, which causes the jump in containment costs (e.g., between 4 and 5 km and between 8 and 9 km).The new models for which containment is feasible describe scenarios with high density of tegus and thus need a higher containment cost.Thus, even if in general cost decreases with an increase of containment radius, in our case there is situation where containment for d 2 is more costly than for d 1 , even if d 1 < d 2 , simply because the estimated cost includes more models for d 2 .But when containment is possible for the same set of models, cost for d 1 is higher than for d 2 .Indeed as it is noted previously, density decreases with distance to the point of first introduction; thus, more tegus will have to be removed for small radius.

Establishment of a containment area
The expected number of remaining tegus outside of the containment area associated with each robust strategy is presented in Fig. 5a.Not surprisingly, the risk-seeking strategy (r c MiniMin ) left the highest expected number of tegus outside of the containment area, while the risk-averse strategy r c Max is the one that minimizes the number of tegus outside of the containment area.Using the strategy r Ã av computed with the average model provided a near optimal result for a containment radius lower than 4 km.For the largest radius, the amount of control is not sufficient and the strategy performs poorly, similar to the risk-seeking strategy r c MiniMin .Indeed, the average model only predicts low tegu density outside of the 50% envelope (see Fig. 3), which means that a low level of removal is used to remove the population within the higher containment radius.One downside of using the average strategy is that it ignores the possibility that tegus can be found outside the 50% envelope, and this strategy is thus not appropriate for the possible models that predict a high density of tegus outside the 50% envelope.In addition, with a containment radius of 1 km, the containment area cannot be implemented given the prediction of the average model.Finally, the risk-averse strategy (r c MiniMax ) and risk-neutral strategy (r c MiniMean ) show similar results, with a slightly better performance for the risk-seeking strategy (r c MiniMean ).These two strategies exhibit similar performance compared to r c Max and are indistinguishable in Fig. 5a.
The performances have to be considered in the context of the associated control cost because the strategies that minimize tegu density outside of the containment area are also the most expensive strategies (see Fig. 5b).Note that the costs of r c MiniMax , r c MiniMean , and r c Max greatly increase for the containment radius ≥5 km.
This increase is explained by the fact that containment becomes achievable for more models, with high predicted density as described earlier.
MiniMin and r Ã av do not take into account these models, and thus, they are not affected and their cost smoothly decreases with the containment radius.

Protection of Everglades National Park and the Turkey Point Power Plant
The minimal CPUE values to protect the ENP and TPPP from tegus over four years is displayed in Appendix S1: Fig. S4.A smaller CPUE seems to be needed to protect ENP.In this region, the CPUE should not have to exceed 8.2 9 10 À4 , 1.48 9 10 À3 , 6.46 9 10 À3 , and 3.06 9 10 À2 during the first, second, third, and fourth years.Given that the current CPUE is estimated to be 0.035, this objective is achievable at least for the next four years with only one trap per cell.For the TPPP, a maximal CPUE of 9.76 9 10 À2 , 4.32 9 10 À2 , 0.15 and 0.608 is needed at time t = 1, . . ., 4. This is achievable for the first three years with the current constraints.But then for the next year, more than five traps per cell is needed.We note that trapping is not only required inside the TPPP, but also around, especially during the first two years, such that the number of animals that will disperse in this area will remain small.

DISCUSSION
This article presents a framework for studying the optimal allocation of control effort with a constrained budget and/or constrained location of controlled cells.This framework is based on linear programming, used to minimize a linear objective function while accounting for linear constraints.Linear programming is particularly useful as it allows us to solve problem that accounts for practical constraints.First, linear programming allows to solve problem of reasonably large size (11,000 square cells of 500 m by 500 m in our case), which is interesting as the area to be managed can be large, such that an entire national park.Second, it allows accounting for constraints and objectives that can be met in practical situations, such that a fixed budget per year, maintaining the population below a given threshold or considering a discrete control variable such that a number of traps.Third, it also allows considering several models of population dynamics, which reflects the uncertainty on the species' biological parameters.Nonetheless, for dynamic problems our framework cannot be used with a logistic growth function, which can be an important limitation.This work builds upon Hof and Bevers (2002) and Epanchin-Niell and Wilen (2012), using a more general population dynamics model and taking advantage of the constraints to model realistic problems.This framework is sufficiently general and flexible to be used in other applications involving control of an invasive fauna or flora.
The use of linear programming to optimize spatial control is made possible by the discrete RD model, allowing us to describe the spatial dynamic of the invasive using linear equations.Reaction-diffusion models rely on simple and general assumptions of the species' population dynamics, which seems to be suitable to our situation, where only little information is available on the species.Of course, if the problem of modeling the spatial distribution of the species was considered alone, there might be no advantage to select RD models instead of, for example, occupancy modeling tools (Bailey et al. 2014) or enhanced spatial distribution models (Sullivan et al. 2012).But here we focus on a coupled problem of spatial control and high uncertainty on the species characteristics, which is why a compromise has to be found.In our case, control has to be implemented, with or without sufficient information on the species, and our framework allows us to use available information in order to estimate the best allocation of control effort.We view the RD model as a null hypothesis, which supposes that the invasion front is a moving circle, centered on the point of first introduction, inside which the density is decreasing with the distance from the center.Unfortunately, it is hard to formulate more hypotheses or to validate them with data, which suggests that our results should be interpreted with caution.For example, Klug et al. (2015) showed that the activity of tegus is different depending on the habitat and the relative position in the invasion area (core vs. periphery).Activity can change the potential of being trapped and our analysis can be improved by defining the CPUE as a function of distance from the point of introduction and habitat.Their study also showed that male dispersal behavior is different in the periphery compared to the core of the population, which might be hard to incorporate in a RD model as it will require a densitydependent dispersal coefficient.On the other hand, one can imagine using a different dispersal coefficient, growth rate, and carrying capacity in the different habitat types.Although Klug et al. (2015) described habitat preferences, there is high uncertainty about current spatial distribution and the only available data are from EDDMaps.When comparing the EDDMaps data to the model predictions (see Fig. 1), there are two major discrepancies: (1) Few observations were recorded east of the ENP, while no observations were recorded in the TPPP, and (2) observations were mostly spread south of the point of first introduction.The EDD-Maps data are issued from a non-regular sampling scheme and some places in the study area are not frequently visited, which might explain the absence of observations in some regions.In addition, mostly private lands can be found north of the presumed point of introduction, which might explain the small number of observations in this zone.Gathering more information on the species could certainly allow us to consider more realistic models and prescribe a better use of the available resources.This highlights an important question, how much should we learn about a species to control it?This is a complicated question because gathering information about the species uses resources, which will not be used to control the population.On the other hand, learning more about the species, for example, the spatial distribution, diet, reproductive and dispersal behavior, could help produce a more efficient control strategy in the future.When the uncertainty can be structured, the Value of Information (Canessa et al. 2015, Williams andJohnson 2015) can be used to quantify the benefit of decreasing the uncertainty on the species' characteristics.Also, control and information gathering are not necessarily contradictory.There is certainly information to be extracted from the control action.For example, if no individuals are captured in the trap, this might change the likelihood of the species being present at this location.In principle, this can be accounted for by including directly some uncertainty measure in the objective function.In this work, we propose considering a large range of population dynamics (i.e., 25 different invasion scenarios) that are likely to include the true dynamics, and then design the control strategy accordingly.
A potentially major limitation of our framework is the hypothesis of an exponential growth when considering population dynamics in the optimization problem.Using an exponential growth when logistic growth is a better assumption tends to overestimate the population.One consequence is that control effort can be wasted in some places where the population is predicted ❖ www.esajournals.orgto be higher than what it is in reality.Thus, the optimal strategy for a population that grows exponentially might not be optimal when used for a population with logistic growth.Another consequence is that control cost will be overestimated.However, the framework that we propose to design a robust containment strategy or to protect of particular areas of interest could be used in an iterative way.In other words, one can use the framework with a time horizon of 1 to compute the optimal strategy for the current year, predict the population density after control using logistic growth, compute the optimal strategy for the next year with the new predicted density as a starting point, and so on.For example, the robust containment strategy can be computed first for 2016 only.The predicted density for 2017 can be computed using Eq. 3, with logistic growth for f and the robust containment strategy for Y.The robust containment strategy for 2017 can then be computed using the predicted density as an input and so on.Finally, the hypothesis of exponential growth is disputable for our case study.We considered 25 different prediction models and for those where the population is predicted to be low, the exponential growth hypothesis might not be a problem.On the contrary, when the population is predicted to be high in 2016, the optimal strategies that we computed might overestimate the control effort and be inefficient in reality.As a consequence, the MiniMax strategy is particularly affected by the exponential growth hypothesis while the impact on the MiniMin strategy is certainly lower.More generally, results of the two last sections (Robust containment strategy and Protection of ENP and the TPPP) should be interpreted with caution.
When applying this framework to Argentine black and white tegus in South Florida, we first showed that at least 25% of the population can be found within approximately a distance of 7 km from the point of first introduction.When at most 300 traps can be used, we showed that they are better to be used within the 25% envelope in order to maximize the number of captures.These results compare favorably with the recent work of Baker (2016), which shows that control should be implemented near the source of the infestation.How traps are allocated in the 25% envelope then depends on the risk attitude: (1) Gather traps close to the point of first introduction (risk-averse and risk-neutral) or (2) distribute traps evenly around the first point of introduction (risk-seeking).Second, it appears that containment might not be possible for some scenarios, but generally the containment success probability increases with the containment radius d.On the contrary, containment cost decreases with an increase of containment radius d.Finally, we discuss the minimal required CPUE that is necessary in the TPPP and in the ENP in order to protect these areas from the invasion of tegus.The situation for ENP is unique, as far as there are few tegus that are currently predicted by some models.In those cases, only a small and achievable CPUE is required for the next four years.For the TPPP, the situation is more critical because some models predict that this area is already invaded by tegus.To be robust to all of the model predictions, the current CPUE is predicted to be sufficient at least for the next three years.For the two following years, more than five traps would need to be placed in the cells inside this area, violating one of our assumptions.Models predict more tegus in the TPPP compared to the ENP because it is accessible via a larger number of cells.Indeed, because we restricted the suitable habitat for tegus, the ENP is only accessible via one path of suitable habitat when the TPPP shared a long boundary with suitable habitat.We also discussed using an average model to inform the management decision.A major downfall of this approach is that it does not allow us to consider extreme, but still possible, cases.For the tegu example, the extreme cases are the models that predict a very high or very low tegu density.As we have shown for the largest width of the containment area, using an average model can result in a decision of poor quality if a pessimistic model is true.Actually, the performance of the strategy computed with the average model is similar to the risk-seeking strategy (see Fig. 5).
It is interesting to note that current control design, resulting from the cooperative trapping of multiple agencies, represents a mix of the different strategies presented in this article: (1) Most of the traps are located in the 25% envelope, which could satisfy the objective of maximizing the number of trapped animals; (2) some traps are located around the TPPP and the ENP, which could satisfy the objective of protecting these ❖ www.esajournals.org15 October 2017 ❖ Volume 8(10) ❖ Article e01979 zones from invasion; and (3) some traps are located near the edge of the 50% envelope, which could satisfy the objective of containment.These results can help managers in deciding the best trap density in the landscape, depending on their prioritization of these objectives and given that our framework is adapted to their situation.More generally, this work allows managers to understand how trap locations and control objectives are linked.To design a control strategy, it is important to first define the objectives (e.g., to maximize capture [M], to contain the species [C], or to protect a delimited zone from invasion [P]).Second, trap locations can be determined.For objective (M), traps should be located in the area of highest density and the number of traps per cell should be decided in order to be able to remove all tegus from the cell.For objective (C), trap should be located inside and near the edge of the containment area, with a decreasing trap intensity from the first point of introduction.The choice of an adapted containment radius should realize a compromise between likelihood of success, cost, and size of the area that should be protected.For objective (P), traps should be located inside the protected zones, but also around when the density of invaders is already high, thus preventing future invasion.In any case, a better prediction of the invader's spatial distribution can highly increase control efficiency and minimize waste of control effort.

ACKNOWLEDGMENTS
Funding was provided by the Invasive Species and Greater Everglades Priority Ecosystem Science programs of the U.S. Geological Survey.We are grateful to anonymous referees for their help in improving this article.Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.The EDDMaps data are publicly available at no charge at https://www.eddmaps.org/ (managed by the University of Georgia).

Fig. 1 .
Fig. 1. (a)The studied areas are colored in green and beige.Green areas are the roads and levees with a 1-km buffer on each side.Beige areas are the agricultural and urban areas.We assumed that these areas represent the only suitable habitats for tegus (Salvator merianae) such that they are not present anywhere else.On the contrary, control is available only in the green areas and the Turkey Point Power Plant (TPPP).The initial point of introduction is represented with a star.Red lines show the boundaries of the Everglades National Park as well as those of the TPPP.Finally, the black circles are the 25, 50, 75, and 99% envelope, centered on the first point of introduction.For example, the 25% envelope is represented with the smallest disk and contained at least 25% of the population for each of the model.(b) Spatial distribution of the EDDMaps data (https://www.eddmaps.org/)used to estimate the dispersal coefficient of the reaction-diffusion model.

Fig. 2 .
Fig. 2. Possible values of the growth rate e (a) and of the standard deviation r of the Gaussian dispersal kernel (b).Five possible values of the growth rate e and their probabilities were first determined.Then for each growth rate, five possible values of r and their probabilities were computed.

Fig. 3 .
Fig. 3. Predicted density of tegus (Salvator merianae) in 2016 in the study area.The prediction of tegus density is provided in km 2 for three different models.The first prediction (a) is obtained by using the least likely combination of the growth rate e and the standard deviation r of the Gaussian dispersal kernel, while the third prediction (c) is obtained by using the most likely combination.The second prediction (b) is obtained by using the average value of the growth rate e and the estimated value of the standard deviation r using the EDDMaps data (https://www.eddmaps.org/).

Fig. 4 .
Fig. 4. (a) Distribution of the expected number of tegus outside of the containment area after 4 yr as a function of the containment width.All of the 25 models are considered.(b) Distribution of the cumulative control cost in order to establish a containment area of different width (containment radius).The distribution of the cumulative of all strategies r cÃ Mm is reported only for the model where containment is possible.The cost is summed over four years.

Fig. 5 .
Fig. 5. (a) Expected number of tegus (Salvator merianae) outside of the containment area, that is, P 25 m¼1 p m Â V c Mm r c ð Þ, for the different control strategies r c and different containment radius after 4 yr of management.The Max, MiniMax, and MiniMean strategies exhibit similar performances and are indistinguishable from the figure.(b) Cumulative cost over four years, in US dollars, of implementing each of the containment strategy.