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

Mathieu Bonneau, Fred Johnson A., Brian J. Smith, Christina M. Romagosa, Julien Martin, Frank J. Mazzotti

► To cite this version:

Mathieu Bonneau, Fred Johnson A., Brian J. Smith, Christina M. Romagosa, Julien Martin, et al.. Optimal control of an invasive species using a reaction-diffusion model and linear programming. Ecosphere, 2017, 8 (10), pp.Non Paginé. 10.1002/ecs2.1979 . hal-02627999

HAL Id: hal-02627999

https://hal.inrae.fr/hal-02627999

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Optimal control of an invasive species using a reaction-diffusion model and linear programming

MATHIEU BONNEAU,<fnref n="1" />,⁴,† FRED A. JOHNSON,<fnref n="2" /> BRIAN J. SMITH,<fnref n="3" /> CHRISTINA M. ROMAGOSA,<fnref n="1" /> JULIEN MARTIN,<fnref n="2" /> AND FRANK J. MAZZOTTI<fnref n="3" />

**Citation:** Bonneau, M., F. A. Johnson, B. J. Smith, C. M. Romagosa, J. Martin, and F. J. Mazzotti. 2017. Optimal control of an invasive species using a reaction-diffusion model and linear programming. Ecosphere 8(10):e01979. 10.1002/ecs2.1979

**Abstract.** Managing an invasive species is particularly challenging as little is generally known about the species' biological characteristics in its new habitat. In practice, removal of individuals often starts before the species is studied to provide the information that will later improve control. Therefore, the locations and the amount of control have to be determined in the face of great uncertainty about the species characteristics and with a limited amount of resources. We propose framing spatial control as a linear programming optimization problem. This formulation, paired with a discrete reaction-diffusion model, permits calculation of an optimal control strategy that minimizes the remaining number of invaders for a fixed cost or that minimizes the control cost for containment or protecting specific areas from invasion. We propose computing the optimal strategy for a range of possible model parameters, representing current uncertainty on the possible invasion scenarios. Then, a best strategy can be identified depending on the risk attitude of the decision-maker. We use this framework to study the spatial control of the Argentine black and white tegus (*Salvator merianae*) in South Florida. There is uncertainty about tegu demography and we considered several combinations of model parameters, exhibiting various dynamics of invasion. For a fixed one-year budget, we show that the risk-averse strategy, which optimizes the worst-case scenario of tegus' dynamics, and the risk-neutral strategy, which optimizes the expected scenario, both concentrated control close to the point of introduction. A risk-seeking strategy, which optimizes the best-case scenario, focuses more on models where eradication of the species in a cell is possible and consists of spreading control as much as possible. For the establishment of a containment area, assuming an exponential growth we show that with current control methods it might not be possible to implement such a strategy for some of the models that we considered. Including different possible models allows an examination of how the strategy is expected to perform in different scenarios. Then, a strategy that accounts for the risk attitude of the decision-maker can be designed.

**Key words:** decision under uncertainty; invasive species; linear programming; reaction-diffusion model; *Salvator merianae*; spatial control.

**Received** 30 August 2016; revised 3 August 2017; accepted 10 August 2017; final version received 8 September 2017.
**Corresponding Editor:** Michael P. Perring.

**Copyright:** © 2017 Bonneau et al. This is an open access article under the terms of the Creative Commons Attribution
License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

⁴Present address: URZ UR143, INRA, 97170, Petit-Bourg (Guadeloupe), France.

† **E-mail:** mathieu.bonneau@inra.fr

<footnote n="1">¹Department of Wildlife Ecology and Conservation, University of Florida, 110 Newins-Ziegler Hall, P.O. Box 110430, Gainesville, Florida 32611 USA</footnote>

<footnote n="2">²U.S. Geological Survey, Wetland and Aquatic Research Center, 7920 NW 71 Street, Gainesville, Florida 32653 USA</footnote>

<footnote n="3">³Ft Lauderdale Research and Education Center, University of Florida, College Avenue, Ft Lauderdale, Florida 33314 USA</footnote>
## 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, 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 and Liebhold 1998, Olson and Roy 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ès 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
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 advect-
tion, which allows consideration of different disper-
sal speeds depending on the dispersal directions.

All of these approaches attempt to minimize an
objective function based on the sum of the ecolog-
ical and economic costs resulting from the inva-
sive species. In other words, the optimal
management option minimizes the management
cost, as well as the ecological and societal dam-
age, over the long term. Using such objective
functions allows the identification of the optimal
strategy that represents a trade-off between man-
agement cost and species damage. Without con-
sidering ecological and societal costs, the optimal
strategy is simple to compute: Do nothing. With-
out considering management cost, it is again sim-
ple 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 esti-
mate from data, it is particularly difficult to quan-
tify the ecological and societal damages (but see
Kaiser and Burnett 2010 for an example). A com-
mon 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 conve-
nient 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 Epan-
chin-Niell and Wilen (2012) allows us to draw
general conclusions about optimal spatial alloca-
tion of effort in stylized cases, but they only con-
sidered the presence/absence of the species and a
simple dispersal model where adjacent cells get
infected in the next time period. Here, we pro-
pose a broader look at the possibilities that the
linear programming framework offers, and illus-
trate its use on the management of the invasive

Argentine black and white tegus (*Salvator meri-
anae*, hereafter “tegus”) in southern Florida. We
considered the problem of spatial allocation of
resources over time, when the invasive popula-
tion has to be minimized under several con-
straints and uncertainty.

We first paired linear programming with a dis-
crete RD model (Levin 1974). Discrete RD models
are relatively general (Holmes et al. 1994), and
we proposed a simulation approach to parame-
terize the dispersal coefficients, based on any dis-
persal kernel. The pattern of invasion proposed
by RD models is simple: In an infinite and suit-
able landscape, invasion is described by a grow-
ing circle centered on the point of first
introduction. Discrete RD models are particu-
larly 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 program-
ming 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.

MATERIALS AND METHODS

Problem definition

Tegus are large lizards native to southeastern
Brazil, Paraguay, Uruguay, and northern Argen-
tina. Individuals were likely introduced via the
pet trade around 2002 in Florida City, FL (Enge-
man 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 (*Crocody-
lus acutus*, Mazzotti et al. 2015), federally recog-
nized 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 croco-
dile and several other native species. In South
Florida, tegus are typically found in disturbed
<figure><img src="image_2.png" /><figcaption>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.</figcaption></figure>

habitats, 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 bru-mating 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):

$$n(i, t) = \sum_{j=1}^{s} [(1 + f(n(j, t - 1)))n(j, t - 1)]g_{ji} \quad (1)$$

where $n(i, t)$ is the population density at site $i \in \{1, \dots, s\}$ at time $t$ and $n(j, t - 1)$ is the density at site $j \in \{1, \dots, s\}$ and 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 \in \{1, \dots, s\}$:
$$f(n(j, t - 1)) = \epsilon \left( 1 - \frac{n(j, t - 1)}{K_j} \right), \quad (2)$$

where $\epsilon$ is the population's intrinsic rate of growth and $K_j$ is the carrying capacity of site j. The term $(1 + f(n(j, t - 1)))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. Then, $[(1 + f(n(j, t - 1)))n(j, t - 1)]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 \neq 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$, $N^t = (n(i, t))_{i=1,...,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 $\epsilon$, the initial population value $N^0$, the carrying capacity $(K_i)_{i=1,...,s}$, and the dispersal coefficient $(g_{ji})_{i=1,...,s}$.

**Intrinsic rate of growth.**—For the rate of growth, we rely on Johnson et al. (2017), where a probability distribution of $\epsilon$ was estimated from expert elicitation. Johnson et al. (2017) estimated that there is a probability of 0.3 that the population's growth rate $\epsilon$ 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.**—We also 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'44.2884" N, 80°28'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.**—The dispersal 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 $\sigma^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 $\sigma$. As for the growth rate $\epsilon$, we propose to fit a probability distribution (a normal distribution here) on $\sigma$ in order to account for uncertainty. We then used a simulation approach to estimate the dispersal coefficients $(g_{ji})_{j=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
dispersal. For example, an asymmetric dispersal kernel can also be considered in the set of possible models.

## Parameter uncertainty

The growth rate $\epsilon$ and the standard deviation $\sigma$ 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 $\epsilon$ and $\sigma$ 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 $\sigma$. In total, we considered 25 different pairs of growth rate and standard deviation of the normal distribution kernel. The models are denoted $M_1, \dots, M_{25}$, associated with the coupled parameters $(\epsilon_1, \sigma_1), \dots, (\epsilon_{25}, \sigma_{25})$ and model probabilities $p_1, \dots, 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:

$$p_1 < p_2 < \dots < p_{25}$$

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%

<figure><img src="image_3.png" /><figcaption>Fig. 2. Possible values of the growth rate $\epsilon$ (a) and of the standard deviation $\sigma$ of the Gaussian dispersal kernel (b). Five possible values of the growth rate $\epsilon$ and their probabilities were first determined. Then for each growth rate, five possible values of $\sigma$ and their probabilities were computed.</figcaption></figure>

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 $\mathbf{Y}^t$, where each component $\mathbf{Y}^t(i)$ represents the amount of control effort used in site $i$. And let us denote $\alpha_i$ the catch per unit effort (CPUE) of the control method in site $i$. When control is implemented, Eq. 1 becomes:

$$n(i, t) = \sum_{j=1}^{s} [(1 + f(n(j, t-1))) \times n(j, t-1)]g_{ji} - \alpha_i \times \mathbf{Y}^t(i). \quad (3)$$

where $\alpha_i \times \mathbf{Y}^t(i)$ represents the number of individuals that have been removed from site $i$.
For tegus, control is implemented using traps and $\mathbf{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 $\alpha_i = \alpha = 0.035$, for all $i = 1, \dots, s$. Trapping is only relevant when tegus are active and thus, the amount of control $\mathbf{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 \le \mathbf{Y}^t(i) \le 5 \times 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 = (\text{US}\$1,000/238)$. It is likely that a non-homogeneous 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 $\pi$ be the vector of the number of traps that are used in the sites. Our problem consists in finding the optimal trapping strategy, $\pi^*$, 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 $\mathbf{R}_{M_m}^t(\pi)$ as the vector of the number of removed tegus in the sites at time $t$ when model $M_m$ is true and strategy $\pi$ is used. The number of remaining tegus after control at time $t$ when model $M_m$ is true, denoted $V_{M_m}(\mathbf{R}_{M_m}^t(\pi))$, is then computed as follows:

$$V_{M_m}(\mathbf{R}_{M_m}^t(\pi)) = \sum_{i=1}^{s} n_{M_m}(i, t), \quad (4)$$

where $n_{M_m}(i, t)$ is obtained by using Eq. 3 with $\alpha \times \mathbf{Y}^t(i) = \mathbf{R}_{M_m}^t(\pi_i)$. Formally, the objective of our problem is to find:

$$\pi^* = \underset{\pi}{\operatorname{argmin}} V_{M_m}(\mathbf{R}_{M_m}^t(\pi)). \quad (5)$$

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 $(\pi_i \times \alpha \times 238)$, which is the maximum number of removed tegus in any site $i$ when $\pi_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$, $\mathbf{R}_{M_m}^t(\pi_i)$, should then verify:

$$\begin{gathered} \mathbf{R}_{M_m}^t(\pi_i) = \min(n_{M_m}(i, t), \pi_i \times \alpha \times 238), \\ \text{for all } i \in \{1, \dots, s\}. \end{gathered} \quad (6)$$

The number of traps should also be bounded and we defined $n_{\text{Max}}^{\text{Trap}}$ as the maximal number of traps that are used such that:

$$\left(\sum_{i=1}^{s} \pi_i\right) \le n_{\text{Max}}^{\text{Trap}}. \quad (7)$$

Currently, nearly $n_{\text{Max}}^{\text{Trap}} = 100$ traps are used, but we considered two other scenarios where $n_{\text{Max}}^{\text{Trap}} = 200$ or $n_{\text{Max}}^{\text{Trap}} = 300$ traps are available.

As explained previously, no more than five traps should be used in a site, such that we have:

$$\text{for all } i \in \{1, \dots, s\}, \quad 0 \le \pi_i \le 5. \quad (8)$$

Finally, traps can be used only on sites where control is allowed:
$$ \text{for all } i \notin C, \quad \pi_i = 0. \tag{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:

$$ \left\{ \begin{array}{ll} \pi^* = \min_{\pi} & V_{M_m} (\mathbf{R}_{M_m}^t (\pi)) \\ \text{such that} & \mathbf{R}_{M_m}^t (\pi_i) = \min(n_{M_m}(i,t), \pi_i \times \alpha \times 238), \\ & \sum_{i=1}^{s} \pi_i \le n_{\text{Max}}^{\text{Trap}}, \\ & \forall i \in \{1, \dots, s\}, \ 0 \le \pi_i \le 5, \\ & \forall i \notin C, \pi_i = 0. \end{array} \right. \tag{10} $$

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, \dots, M_{25}$, providing 25 different strategies $\pi_{M_1}^*, \dots, \pi_{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 $\pi_{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 $\pi_{M_1}^*, \dots, \pi_{M_{25}}^*$, the MiniMax strategy $\pi_{\text{MiniMax}}$ is a robust strategy that minimizes the maximal possible remaining number of tegus, over all possible models:

$$ \pi_{\text{MiniMax}} = \underset{\pi_{M_m}}{\operatorname{argmin}} \left\{ \max_{M_{m'}} \left[ V_{M_{m'}} \left( \mathbf{R}_{M_{m'}}^t \left( \pi_{M_m}^* \right) \right) \right] \right\}. \tag{11} $$

In other words, the value of each strategy $\pi_{M_m}^*$ is computed in case any of the possible models $M_{m'}$ is the true model. For example, $V_{M_{m'}}(\pi_{M_m}^*)$ is the remaining number of tegus when strategy $\pi_{M_m}^*$ is used but model $M_{m'}$ is the true model.

Then, $\max_{M_{m'}}[V_{M_{m'}}(\pi_{M_m}^*)]$ is the worst possible performance of the strategy $\pi_{M_m}^*$ over all possible models. As a consequence, $\pi_{\text{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, $\pi_{\text{MiniMean}}$, is the strategy that achieves the minimal expected number of remaining tegus:

$$ \pi_{\text{MiniMean}} = \underset{\pi_{M_m}^*}{\operatorname{argmin}} \left\{ \sum_{M_{m'}} p_{M_{m'}} \times V_{M_{m'}} (\pi_{M_m}^*) \right\}. \tag{12} $$

In other words, $\pi_{\text{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, $\pi_{\text{MiniMin}}$, is the strategy that minimizes the minimal possible remaining number of tegus over all possible models:

$$ \pi_{\text{MiniMin}} = \underset{\pi_{M_m}^*}{\operatorname{argmin}} \left\{ \max_{M_{m'}} \left[ V_{M_{m'}} \left( \mathbf{R}_{M_{m'}}^t \left( \pi_{M_m}^* \right) \right) \right] \right\}. \tag{13} $$

In other words, $\pi_{\text{MiniMin}}$ is the strategy that minimizes the remaining number of tegus in the best-case scenario. This is the opposite of $\pi_{\text{MiniMax}}$ and it should be selected by a risk-seeking decision-maker. Finally, we also computed the optimal control effort using an average model, $M_{\text{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 $\sigma = 3$ and the growth rate $\epsilon = 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
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.*—We seek to implement a disk-shaped 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:

$$f(n(j, t - 1)) = \epsilon n(j, t - 1). \quad (14)$$

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 $\pi_{M_m}^c, \dots, \pi_{M_m}^{c+4}$ that minimize the control cost $\sum_{t'=t}^{t+4} (\sum_{i=1}^n c_u \times \pi_{M_m}^{c'})$. 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, $\pi_{M_m}^{c'} \le 238 \times 5$ for all time $t'$, and (5) the number of removed tegus does not exceed the number of tegus in the site, that is, $\alpha \times \pi_{M_m}^{c'} \le n_{M_m}(t', i)$ for all time $t'$. 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, $\pi_{M_m}^{c'}$ 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 $\pi_{M_m}^c$ 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 $V_{M_{m1}}^c(\pi_{M_{m2}}^c)$ the number of remaining tegus outside of the containment area when using strategy $\pi_{m2}^c$ but model $M_{m1}$ is the true model. Note that we have $V_{M_{m1}}^c(\pi_{M_{m1}}^c) = 0$ when 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, $\pi_{MiniMax}^c$, which minimizes the number of tegus outside of the containment area in the worst-case scenario:
$$ \pi_{\text{MiniMax}}^{c} = \arg \min_{\pi_{M_{m}}^{c*}} \left\{ \max_{M_{m'}} \left[ V_{M_{m'}}^{c} \left( \pi_{M_{m}}^{c*} \right) \right] \right\}. \quad (15) $$

The MiniMean strategy $\pi_{\text{MiniMean}}^c$ which minimizes the expected number of tegus outside of the containment area:

$$ \pi_{\text{MinMean}}^{c} = \arg \min_{\pi_{M_{m}}^{c}} \left\{ \sum_{M_{m'}} p_{M_{m'}} \times V_{M_{m'}}^{c} \left( \pi_{M_{m}}^{c*} \right) \right\}. \quad (16) $$

We also considered the MiniMin strategy, $\pi_{\text{MiniMin}}^c$, which minimizes the minimal possible remaining number of tegus over all possible models:

$$ \pi_{\text{MiniMin}}^{c} = \arg \min_{\pi_{M_{m}}^{c}} \left\{ \min_{M_{m'}} \left[ V_{M_{m'}}^{c} \left( \pi_{M_{m}}^{c*} \right) \right] \right\}. \quad (17) $$

In addition to these three strategies, we considered a Max strategy $\pi_{\text{Max}}^c$ which on any time step and any site, use the largest control effort among all the possible strategies:

$$ \forall t' = t \dots t + 4, \quad i = 1, \dots, s \\ \pi_{\text{Max}}^{c'}(i) = \max_{m=1 \dots 25} \pi_{M_m}^{c*'}(i), \quad (18) $$

where $\pi_{M_{m}}^{c*'}(i)$ is the amount of effort in site i during year $t'$ when strategy $\pi_{M_{m}}^{c*}$ is used. Finally, we also computed the optimal control effort for the containment of tegus using the average model $M_{\text{av}}$ as described before.

Note that to compute the remaining number of tegus outside of the containment area, $V_{M_{m}}^{c}(\cdot)$, 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 worst-case 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_{\text{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². 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_{\text{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
<figure><img src="image_4.png" /><figcaption>Fig. 3. Predicted density of tegus (*Salvator merianae*) in 2016 in the study area. The prediction of tegus density is provided in km² for three different models. The first prediction (a) is obtained by using the least likely combination of the growth rate ε and the standard deviation σ 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 ε and the estimated value of the standard deviation σ using the EDDMaps data (https://www.eddmaps.org/).</figcaption></figure>

minimized the maximal values over all models; that is, there are 25 other possible σ<sub>MiniMax</sub> strategies. As explained earlier, the σ<sub>MiniMax</sub> 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 α × 5 × 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 σ<sub>MiniMean</sub> and σ<sub>MiniMin</sub> 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 (σ<sub>MiniMax</sub>) and risk-neutral (σ<sub>MiniMean</sub>) 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 (σ<sub>MiniMin</sub>) 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
<figure><img src="image_5.png" /><figcaption>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 $\sigma_{M_m}^*$ is reported only for the model where containment is possible. The cost is summed over four years.</figcaption></figure>

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 $\sigma_{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 containment 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 ($\sigma_{MiniMin}^c$) left the highest expected number of tegus outside of the containment area, while the risk-averse strategy $\sigma_{Max}^c$ is the one that minimizes the number of tegus outside of the containment area. Using the strategy $\sigma_{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 $\sigma_{MiniMin}^c$. 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 ($\sigma_{MiniMax}^c$) and risk-neutral strategy ($\sigma_{MiniMean}^c$) show similar results, with a slightly better performance for the risk-seeking strategy ($\sigma_{MiniMean}^c$). These two strategies exhibit similar performance compared to $\sigma_{Max}^c$ and are indistinguishable in Fig. 5a.

The performances have to be considered in the context of the associated control cost because the
<figure><img src="image_6.png" /><figcaption>Fig. 5. (a) Expected number of tegus (*Salvator merianae*) outside of the containment area, that is, $\sum_{m=1}^{25} p_m \times V_{M_m}^c (\sigma^c)$, for the different control strategies $\sigma^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.</figcaption></figure>

strategies that minimize tegu density outside of the containment area are also the most expensive strategies (see Fig. 5b). Note that the costs of $\sigma_{\text{MiniMax}}^c$, $\sigma_{\text{MiniMean}}^c$, and $\sigma_{\text{Max}}^c$ greatly increase for the containment radius $\ge 5$ km.

This increase is explained by the fact that containment becomes achievable for more models, with high predicted density as described earlier. $\sigma_{\text{MiniMin}}^c$ and $\sigma_{\text{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 \times 10^{-4}$, $1.48 \times 10^{-3}$, $6.46 \times 10^{-3}$, and $3.06 \times 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 \times 10^{-2}$, $4.32 \times 10^{-2}$, 0.15 and 0.608 is needed at time $t = 1, \dots, 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 situ-
ation, where only little information is available
on the species. Of course, if the problem of mod-
eling the spatial distribution of the species was
considered alone, there might be no advantage to
select RD models instead of, for example, occu-
pancy modeling tools (Bailey et al. 2014) or
enhanced spatial distribution models (Sullivan
et al. 2012). But here we focus on a coupled prob-
lem of spatial control and high uncertainty on
the species characteristics, which is why a com-
promise 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 cir-
cle, 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 incor-
porate in a RD model as it will require a density-
dependent 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 pre-
dictions (see Fig. 1), there are two major discrepan-
cies: (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 addi-
tion, 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 con-
trol 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 distribu-
tion, diet, reproductive and dispersal behavior,
could help produce a more efficient control strat-
egy in the future. When the uncertainty can be
structured, the Value of Information (Canessa et al.
2015, Williams and Johnson 2015) can be used to
quantify the benefit of decreasing the uncertainty
on the species' characteristics. Also, control and
information gathering are not necessarily contra-
dictory. 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 mea-
sure in the objective function. In this work, we pro-
pose 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 frame-
work 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 assump-
tion tends to overestimate the population. One
consequence is that control effort can be wasted
in some places where the population is predicted
to 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 overesti-
mated. 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 strat-
egy for the next year with the new predicted den-
sity as a starting point, and so on. For example,
the robust containment strategy can be com-
puted first for 2016 only. The predicted density
for 2017 can be computed using Eq. 3, with logis-
tic growth for *f* and the robust containment strat-
egy 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 popu-
lation 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 sec-
tions (*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 imple-
mented 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-seek-
ing). 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 min-
imal 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 acces-
sible 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 deci-
sion. 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 strat-
egy (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 dif-
ferent 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
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.eddmapsc.org/ (managed by the University of Georgia).

## LITERATURE CITED

* Allen, L. J. 1987. Persistence, extinction, and critical patch number for island populations. Journal of Mathematical Biology 24:617–625.

* Bailey, L. L., D. I. MacKenzie, and J. D. Nichols. 2014. Advances and applications of occupancy models. Methods in Ecology and Evolution 5:1269–1279.

* Baker, C. 2016. Target the source: optimal spatiotemporal resource allocation for invasive species control. Conservation Letters. https://doi.org/10.1111/conl.12236

* Blackwood, J., A. Hastings, and C. Costello. 2010. Cost-effective management of invasive species using linear-quadratic control. Ecological Economics 69: 519–527.

* Bonneau, M., F. Johnson, and C. Romagosa. 2016. Spatially explicit control of invasive species using a reaction-diffusion model. Ecological Modelling 337:15–24.

* Canessa, S., G. Guillera-Arroita, J. J. Lahoz-Monfort, D. M. Southwell, D. P. Armstrong, I. Chadès, R. C. Lacy, and S. J. Converse. 2015. When do we need more data? A primer on calculating the value of information for applied ecologists. Methods in Ecology and Evolution 6:1219–1228.

* Carlton, J. 2003. Invasive species: vectors and management strategies. Island Press, Washington, D.C., USA.

* Chadès, I., T. Tarnpolskaya, S. Dunstall, J. Rhodes, and A. Tulloch. 2015. A comparison of adaptive management and real options approaches for environmental decisions under uncertainty. http://espace.library.uq.edu.au/view/UQ:383967

* Engeman, R., E. Jacobson, M. L. Avery, and W. E. Meshaka. 2011. The aggressive invasion of exotic reptiles in Florida with a focus on prominent species: a review. Current Zoology 57:599–612.

* Epanchin-Niell, R. S., and J. E. Wilen. 2012. Optimal spatial control of biological invasions. Journal of Environmental Economics and Management 63:260–270.

* Hauser, C., and A. McCarthy. 2009. Streamlining ‘search and destroy’: cost-effective surveillance for invasive species management. Ecology Letters 12:683–692.

* Hof, J. G., and M. Bevers. 2002. Spatial optimization in ecological applications. Columbia University Press, University of Columbia, New York, New York, USA.

* Holmes, E. E., M. A. Lewis, J. Banks, and R. Veit. 1994. Partial differential equations in ecology: spatial interactions and population dynamics. Ecology 75:17–29. http://www.jstor.org/stable/1939378

* Johnson, F. A., B. J. Smith, M. Bonneau, J. Martin, C. Romagosa, F. Mazzotti, H. Waddle, R. N. Reed, J. K. Eckles, and L. J. Vitt. 2017. Expert elicitation, uncertainty, and the value of information in controlling invasive species. Ecological Economics 137: 83–90.

* Kaiser, B., and K. Burnett. 2010. Spatial economic analysis of early detection and rapid response strategies for an invasive species. Resource and Energy Economics 32:566–585.

* Klug, P. E., R. N. Reed, F. J. Mazzotti, M. A. McEachern, J. J. Vinci, K. K. Craven, and A. A. Y. Adams.
2015. The influence of disturbed habitat on the spatial ecology of Argentine black and white tegu (Tupinambis merianae), a recent invader in the Everglades ecosystem (Florida, USA). Biological Invasions 17:1785-1797.

Leung, B., D. Finnoff, J. F. Shogren, and D. Lodge. 2005. Managing invasive species: rules of thumb for rapid assessment. *Ecological Economics* 55: 24-36.

Leung, B., D. M. Lodge, D. Finnoff, J. F. Shogren, M. A. Lewis, and G. Lamberti. 2002. An ounce of prevention or a pound of cure: bioeconomic risk analysis of invasive species. Proceedings of the Royal Society of London B: Biological Sciences 269:2407-2413.

Levin, S. A. 1974. Dispersion and population interactions. American Naturalist 108:207-228. http://www.jstor.org/stable/2459851

Mack, R. N., D. Simberloff, W. Mark Lonsdale, H. Evans, M. Clout, and F. A. Bazzaz. 2000. Biotic invasions: causes, epidemiology, global consequences, and control. Ecological Applications 10:689-710.

Marescot, L., G. Chapron, I. Chadès, P. L. Fackler, C. Duchamp, E. Marboutin, and O. Gimenez. 2013. Complex decisions made simple: a primer on stochastic dynamic programming. Methods in Ecology and Evolution 4:872-884.

Mazzotti, F. J., M. McEachern, M. Rochford, R. N. Reed, J. K. Eckles, J. Vinci, J. Edwards, and J. Wasi-lewski. 2015. *Tupinambis merianae* as nest predators of crocodilians and turtles in Florida, USA. Biological Invasions 17:47-50.

McCleery, R. A., A. Sovie, R. N. Reed, M. W. Cunningham, M. E. Hunter, and K. M. Hart. 2015. Marsh rabbit mortalities tie pythons to the precipitous decline of mammals in the Everglades. Proceedings of the Royal Society of London B: Biological Sciences. http://rspb.royalsocietypublishing.org/content/282/1805/20150120. https://doi.org/10.1098/rspb.2015.0120

Mehta, S., R. Haight, F. Homans, S. Polasky, and R. Venette. 2007. Optimal detection and control strategies for invasive species management. *Ecological Economics* 61:237-245.

Miller III, A. C., and T. R. Rice. 1983. Discrete approximations of probability distributions. Management Science 29:352-362.

Moore, J. L., M. C. Runge, B. L. Webber, and J. R. Wilson. 2011. Contain or eradicate? Optimizing the management goal for Australian acacia invasions in the face of uncertainty. Diversity and Distributions 17:1047-1059.

Olson, L. J., and S. Roy. 2002. The economics of controlling a stochastic biological invasion. American Journal of Agricultural Economics 84:1311-1316. http://www.jstor.org/stable/1245064

Parker, I., et al. 1999. Impact: toward a framework for understanding the ecological effects of invaders. Biological Invasions 1:3-19.

Pimentel, D., R. Zuniga, and D. Morrison. 2005. Update on the environmental and economic costs associated with alien-invasive species in the United States. Ecological Economics 52:273-288.

Sharov, A. A., and A. M. Liebhold. 1998. Bioeconomics of managing the spread of exotic pest species with barrier zones. Ecological Applications 8:833-845.

Simberloff, D., D. C. Schmitz, and T. C. Brown. 1997. Strangers in paradise: impact and management of nonindigenous species in Florida. Island Press, Washington, D.C., USA.

Sullivan, M. J., R. G. Davies, L. Reino, and A. Franco. 2012. Using dispersal information to model the species-environment relationship of spreading non-native species. Methods in Ecology and Evolution 3:870-879.

Williams, B. K., and F. A. Johnson. 2015. Value of information and natural resources decision-making. Wildlife Society Bulletin 39:488-496.

## SUPPORTING INFORMATION

Additional Supporting Information may be found online at: http://onlinelibrary.wiley.com/doi/10.1002/ecs2.1979/full