Optimizing strategies for slowing the spread of invasive species

Invasive species are spreading worldwide, causing damage to ecosystems, biodiversity, agriculture, and human health. A major question is, therefore, how to distribute treatment efforts cost-effectively across space and time to prevent or slow the spread of invasive species. However, finding optimal control strategies for the complex spatial-temporal dynamics of populations is complicated and requires novel methodologies. Here, we develop a novel algorithm that can be applied to various population models. The algorithm finds the optimal spatial distribution of treatment efforts and the optimal propagation speed of the target species. We apply the algorithm to examine how the results depend on the species’ demography and response to the treatment method. In particular, we analyze (1) a generic model and (2) a detailed model for the management of the spongy moth in North America to slow its spread via mating disruption. We show that, when utilizing optimization approaches to contain invasive species, significant improvements can be made in terms of cost-efficiency. The methodology developed here offers a much-needed tool for further examination of optimal strategies for additional cases of interest.


Introduction
Due to global change, many species have become established outside their native habitats [1][2][3].Some of these invasive species are harmful, inflicting severe damage to agriculture, forestry, biodiversity, human health, and infrastructure [1,2,[4][5][6][7].Commonly, such invasive species begin their infestation from a certain area where they become established and spread outward [8][9][10].Biosecurity measures often combine prevention, detection, and early eradication to prevent the establishment of invasive species [11][12][13].If these efforts fail and the invasive species becomes established, it often becomes practically impossible to entirely remove it, but it may still be possible to contain it to prevent or slow down its further spread to nearby areas [9,10,14,15].Accordingly, containment projects are implemented worldwide to slow the propagation speed of various invasive populations.To name but a few examples, this includes populations of various insects, such as the spongy moth (Lymantria dispar) [15,16] and the emerald ash borer (Agrilus planipennis) in the US [16,17]; various marine species, such as the crayfish species Procambarus clarkia [18] and zebra mussels (Dreissena polymorpha) in the US and in Europe [19]; and numerous invasive plants worldwide [20].
In turn, containment often necessitates treatment over large areas and includes surveillance and eradication in regions close to the invasion area, suppression of the population in regions where it has already been established, and limiting the dispersal of the species into uninvaded areas at the borders between these regions to reduce propagule pressure, where all these treatment strategies are interrelated and can be used simultaneously [21].A major question is, therefore, how managers should distribute the treatment efforts over space and time to minimize the net cost due to both treatment and damage.In particular, at what speed should the manager allow the species to propagate in order to minimize the net cost, and what is the optimal spatial allocation of the treatment along the front line where the species propagates?Researchers have integrated knowledge about population dynamics and economic costs into bioeconomic models of invasive species [10,[22][23][24]; but to date, the primary focus of that literature has been on management in relatively small geographic areas.In such cases, the manager needs to either eradicate the population (remove it entirely from the target area) [25][26][27][28][29] or control it at a certain density [30][31][32].Complete eradication is generally feasible only if the species' density is still low and the species is easily detectable even when its density is low [25,33], or if the species is subject to an Allee effect, namely, its population declines or grows more slowly when its density is low [28,29,34,35].
Several bioeconomic studies have explored the optimal containment of an invasive species to stop or slow its spread.Some of these studies have examined where to build fences to contain the species [36], some examined the optimal strategies for the surveillance of small populations of the species [37][38][39], and some examined the optimal distribution of treatment efforts to completely stop the spread of the species [21,[40][41][42][43].In turn, models that examined the optimal propagation speed have considered the cost of slowing the spread to a given speed, without an explicit model of the population's spatial dynamics and treatment [44].
One factor that limits studies that explore the optimal spatial distribution of treatments for slowing population spread is that standard optimization algorithms are limited in their capacity to solve spatially extended dynamical optimization problems over large areas.For example, Pontryagin's maximum principle [32,45] has been used to study certain spatially-extended bioeconomic models [41,46], but its application for solving more complex models requires its implementation specifically tailored to the specific model, and it might not be applicable in practice to all cases and boundary conditions, such as cases that incorporate dispersal kernel instead of diffusion.Also, methods like Dynamic Programming, which work well for finding the optimal treatment of local populations [47], often operate too slowly for solving models that incorporate large landscapes, and therefore, approximation algorithms might be necessary [48].
Here, we consider a population of an invasive species that spreads over space via some dispersal kernel, and we develop a novel algorithm that finds, at least approximately, (a) the optimal distribution of the treatment over space to slow the spread of the population to a certain speed and (b) the optimal target speed at which the species should be allowed to propagate.In particular, the target speed could be (a) zero (stopping the spread), (b) the natural propagation speed of the species (no treatment), (c) some other positive speed (slowing the spread), or (d) some negative speed (reversing the spread or eradicating).
The algorithm is general and enables the examination of a large variety of models.Here we demonstrate how the algorithm is used to conduct a comparative analysis, examining how the response of the species to the treatment method affects the optimal treatment strategy.A particular case study that we analyze is the management of the spongy moth (Lymantria dispar) to slow its population spread in North America [15].In that case, we consider mating disruption as the main treatment method.In addition, we analyze a general and simpler model to further examine how the response of the species to the treatment affects the optimal treatment.One case that we examine using the general model is where the number of individuals abated per dollar invested does not depend on the species' density.This assumption is reasonable in models of large plants that are observable and thus easily detectable.It has been applied in various models of optimal eradication of Spartina alterniflora [25] and hybrid Spartina [27].We also examine cases where the number of individuals abated per dollar invested decreases as the species density declines.This characterizes scenarios in which it is difficult to detect the species' individuals and an aerial treatment is necessary, such as using pesticides for treating insects [26,29,32].

Defining the problem
We study an invasive species population with a density n(x,t) that can vary over a one-dimensional space axis, x, and time, t.The investment in treatment that suppresses the population, A(x, t), can also vary over space and time.We assume that, initially, n has the shape of a population "front," in which the population has reached its carrying capacity in a certain region (where the x is large), and has not yet invaded regions that are further away from it (n is close to zero where x is small).Without any treatment (A = 0), the population propagates to the uninvaded area, such that the population front approaches a shape that propagates leftward at a certain speed, v 0 .However, by applying treatment (A > 0), the manager can slow down the propagation speed to some v < v 0. The manager might even stop the speed entirely (v = 0), or reverse the direction of the front's propagation such that the invasive species would retreat (v < 0).
The goal of the algorithm developed here is to find the treatment function A that slows the propagation speed down to a given, constant speed v, while using the minimal amount of annual investment.In turn, by using the algorithm to examine various values of v, we can obtain a comprehensive analysis of the best strategy, revealing the optimal target speed v and the distribution of treatment over space to achieve that speed.
Before we describe the algorithm, we describe here the bioeconomic model, which comprises: (a) a model of the dynamics of the invasive species, including its response to the treatment; (b) an objective function, defining the goal of the manager (minimizing the total costs due to treatment and damage); and (c) the constraints-restrictions on the set of possible controls.The algorithm is general and could be applied to a wide variety of cases, and accordingly, we define here a general framework defining the dynamics of n.At the same time, we describe in detail how we parametrize the case study of the spongy moth management in North America and we demonstrate how the algorithm is applied to that case.

Population dynamics
We consider a widely-used population model in which the population is spread over space and disperses according to a kernel function [49][50][51].The underlying assumption is that, in addition to the adult population, the population may comprise propagules (juveniles, larvae, or seeds) that disperse over space; however, the development from propagules to adult is relatively fast, allowing us to consider a single variable, n(x, t), characterizing the adult population density.We assume that treatments affect adults or propagules after their dispersal and we distinguish between (a) continuous-time dynamics, characterizing a population with overlapping generations, and (b) discrete-time dynamics, characterizing a population with no overlapping generations, such as that of univoltine insects or annual plants.
Specifically, if the dynamics are time-continuous, they are given by [21,49,50,52] dnðx; tÞ dt where d(n) is the rate of adult natural mortality, R(n, A) is the rate of adult removal due to the treatment, b(n) is the rate of propagule production, and G is the dispersal kernel that determines how propagules disperse over space [51].Note that the rates d, R, and b vary over time and space due to their dependency on n(x, t) and A(x, t).In turn, if the dynamics are time-discrete, they are given by This characterizes a population in which the new generation replaces the old one and dispersal occurs at some particular life stage, commonly as seed or larva, and t is measured in units of generation, commonly a year.In the discrete-time scenario, we assume that b is sufficiently small, such that the local dynamics converge to some steady state and do not exhibit oscillatory or chaotic dynamics.

Functional forms-simple model
The method developed here for solving the model is general, but to demonstrate the results, we consider particular functional forms in the simulations.We begin here with a simple and general model, while in the next sub-section we describe the model that applies to the spongy moth.For the simple model, we consider continuous-time dynamics (Eq (1)) with the following birth and death rates [21]: where γ is the rate of adult natural mortality, r is the per-capita rate of juvenile production when the population is small, and k is the carrying capacity.We consider the following, widely-used form of the dispersal kernel function [51]: where σ is the characteristic dispersal length.
In turn, the response of the population to the treatment, R(n, A), is monotone increasing with n and A and satisfies R(n, 0) = 0. We consider the following, general form, that has been used previously in the literature [41,53]: where βn α is the marginal efficiency of the treatment (number of individuals removed per dollar investment), and 0�α�1 is the species' invisibility, determining how the efficiency of the treatment diminishes as n decreases.Specifically, α = 0 characterizes cases in which the number of individuals abated per dollar invested does not depend on the species' density.This functional form has been applied in various models of optimal eradication of Spartina alterniflora [25] and hybrid Spartina [27].These plants are easily detectable, and consequently, the cost of abating a given area covered by the species does not depend on how much of the total area has already been removed [25,27].On the other hand, α > 0 characterizes cases in which the number of individuals abated per dollar invested decreases as the species density declines [41,53].In particular, α = 1 characterizes cases in which the efficiency is proportional to n, which has been used to model treatment via pesticides that remove a certain fraction of the population in a given treatment [29,32].

Functional forms-spongy moth population model
The spongy moth is a univoltine insect whose new generation replaces the previous one each year.Accordingly, its population can be characterized by the discrete-time dynamics, Eq (2).
To slow the spread of the spongy moth in North America, mating disruption has emerged as a prevalent strategy [26,54].Mating disruption involves the distribution of synthetic pheromone flakes throughout the affected zones during the mating season.The flakes mimic the female moth's natural scent, leading male moths on a futile quest and significantly hindering their ability to locate actual females.Consequently, fewer females are fertilized, which results in lower fecundity.
To incorporate the effect of mating disruption, we need to consider a local birth function, b(n, A), that integrates this effect.To derive the form of b, we first derive an expression for the probability that a given female is fertilized (i.e., found by at least one male), P(n, A), which depends on the population density, n, and the level of investment in mating disruption, A, in a given location.This probability can be calculated by assuming that males search for females independently from one another [26,29].Then, the probability that a female is not found by any male is given by the Poisson coefficient exp(-λn m ), where λ 0 is the probability that a given male finds a given female during the mating season, and n m is the number of males.
In turn, since the moth's sex ratio is approximately 1:1, we assume that n m = n/2.We also assume, as in previous studies [29], that λ declines with the investment in mating disruption and is given by λ = λ 0 /(1+aA), where λ 0 is the probability that a given male finds a given female without interference, and a > 0 is determined by the cost and the efficiency of the false pheromone flakes [29].It follows that In turn, each fertilized female lays one egg mass, which implies that the number of egg masses, n 0 , produced by the n/2 females, is given by At low population densities, each egg mass develops, on average, into r females (2r adults).At higher densities, due to resource limitations at the larval stage, fewer eggs develop into adults.
To incorporate this, we consider the following birth function: The parameter values are estimated from the literature.Previous studies estimated that r is between 2 and 10 [55,56], the carrying capacity is in the range of 1000 -3000 egg masses per hectare [14], and a is roughly 0.08 USD -1 per hectare [26,29].Studies have also shown that the annual dispersal distance of the spongy moth population, σ, ranges from 2.5 to 20 km [57] considering the dispersal kernel described in Eq (4).Lastly, to estimate λ 0 note that the dynamics described in Eqs (6)(7)(8) result in an Allee effect because, if n is below a certain threshold, n a , such that b < 1, the population declines.Empirical evidence shows that the Allee threshold is about 2-4 egg masses per hectare [58], namely, n a = 4r -8r adults per hectare.The Allee threshold appears where b = 1, which implies that λ 0 = (2/n a )ln (r/(r -1)) [29].For values of r between 2 and 10, this implies that λ 0 is between 0.17 and 0.0026.Note that some of these parameters could be scaled out as they could affect the units but not the underlying dynamics.Specifically, σ only determines the units of x; a only determines the units of A; and n can be measured in units of the population's carrying capacity, in which case the single parameter kλ 0 would displace the parameters k and λ 0 .

Cost of treatment, objective function, and constraints
Our goal is to find the treatment function A(x, t) that minimizes the net cost of treatment over time while ensuring that n propagates leftward at a rate that does not exceed v (where v is some target speed that is pre-determined and could itself be subject optimization, as described below).Specifically, A also shifts leftward at a rate v, alongside the front, such that it has the form where ÂðxÞ is stationary and does not depend on t.In turn, the annual cost of treatment (ACT) is defined as the cost of the treatment during one year over all spatial locations: where the cost due to one unit of investment is normalized to one and could reflect the cost of expenditures plus the cost due to environmental side effects from the treatment.If A is given by Eq (9), the ACT is the same in all years (because Â does not depend on t), and n approaches a quasi-equilibrium shape that only moves leftward along the x -axis at a rate v: as t ! 1, where nðxÞ has the shape of a population front.The objective is, therefore, to find a treatment function ÂðxÞ that minimizes the ACT under the constraint that Aðx; tÞ ¼ Âðx þ vtÞ would keep the front propagating at a rate v (i.e., n is given by Eq (11) as t ! 1).
We denote this optimal treatment, ÂðxÞ, as A opt (x), and the population density nðxÞ associated with A opt (x) as n opt (x).On top of that, note that the treatment and the population density must be non-negative (a negative treatment or population size would be biologically meaningless).Therefore, there are additional constraints, implying that for all x.In particular, Eq (12B) implies that ÂðxÞ must not be too large, and the exact restrictions that Eq (12B) imposes on ÂðxÞ depend on the particular form of the dynamics.in the space of all front shapes.To determine whether this local optimum is indeed the global optimum, the algorithm is applied to various different initial front shapes to verify that it always converges to the same solution.Although there is no guarantee that the algorithm will always converge to the global optimum, it at least provides an approximation to that optimum.Finally, since the algorithm gets v as an input, we use it to explore the entire range of possible propagation speeds and find the optimal v.Note that this algorithm extends the algorithm developed in [21] to cases in which v 6 ¼ 0, and also to the study of more general response functions of the species to the treatment, such as those incorporated in Eqs ( 5) and (6)(7)(8).The complete code can be found on Zenodo [59].Results as raw data can be found on Dryad [60].

Algorithm: General framework
The algorithm initiates with some heuristic "guess" of a front, ñðxÞ, which could be any function that is monotonically increasing with x, approaches 0 as x !-1, and approaches the population's carrying capacity as x ! 1.In turn, the algorithm implements a function "find ÂðñÞ" (see below) that calculates Â, the function that minimizes the ACT while guaranteeing that ñ does not travel faster than v in any location (e.g., if Aðx; tÞ ¼ ÂðxÞ and nðx; 0Þ ¼ ñðxÞ, then dnðx; 0Þ=dt � v � dnðx; 0Þ=dx for all x).The algorithm also implements a function "find nð ÂÞ" (see below) that finds the front shape, n, that is approached in the long run if Aðx; tÞ ¼ Âðx þ vtÞ.(Note that n may differ from ñ that is used to find Â because there could be locations x where ñðxÞ propagates slower than v.) In the initial step, the algorithm uses these functions to calculate Â associated with the initial front shape, and then to find n associated with that Â. Next, the algorithm examines several front shapes that are slight modifications of n.Specifically, for numerical purposes, the algorithm iterates through all values of x between x min and x max at some fine resolution Δx, where x min is sufficiently small (nðx min Þ � 0) and x max is sufficiently large (nðx max Þ is close to the population's carrying capacity).At each location x 0 within that range, the algorithm evaluates two modifications of n: (a) n " (x 0 ), which differs from nðxÞ only at x = x 0 and its value there is between nðx 0 Þ and nðx 0 þ DxÞ, and (b) n # (x 0 ), which also differs from n only at x = x 0 , where its value is between nðx 0 Þ and nðx 0 À DxÞ.Specifically, , where p = 1,2,3. . .defines the search's resolution.
In turn, for each modification, the algorithm calculates the associated treatment function Â using "find ÂðñÞ", where ñ is the examined modifications.If the modification results in an improvement ( Â associated with n " (x 0 ) or n # (x 0 ) for some x 0 results in a lower ACT than that associated with the original n), then the algorithm calculates n that is associated with that modification (using "find nð ÂÞ") and continues with that new function as n.The algorithm repeats that step and continues to modify n and Â, until it converges to a value n and Â that cannot be improved upon with local modifications, meaning any further local changes to n would result in a larger ACT.Also, the algorithm begins with a resolution p = 1 and keeps increasing that resolution until it does not find any improvement, even with the finer resolutions.We denote these final n and Â as n opt and A opt , respectively.
Algorithm: Finding Â that keeps a given front shape ñ from propagating faster than v in all locations ("find ÂðñÞ") Slow the spread, v > 0. First, note that the input function, ñðxÞ, can be any function with a shape of a "front": ñðxÞ is monotone increasing with x ðdñ=dx � 0), approaches 0 as x !-1, and approaches the population's carrying capacity as x ! 1.Initially, the algorithm determines the time and space intervals, Δt and Δx, respectively.If the dynamics are time-discrete, then the time intervals are always one generation (Δt = 1), and accordingly, the algorithm determines Δx such that v = Δx/Δt, i.e., Δx = v.If the dynamics are time-continuous, Δt is determined as the time during which the front would propagate a given Δx if it propagates at a rate v (Δt = Δx/v).
In turn, the algorithm simulates the dynamics of n for a period Δt with A(x) = 0 for all x, where initially nðx; 0Þ ¼ ñðxÞ.We denote the resulting function as n 1 (x) (i.e., n1 (x) = n(x, Δt)).The requirement that the speed cannot exceed v implies that, if n 1 ðx 1 Þ > ñðx 1 þ DxÞ at some location x = x 1 , then treatment must be used to suppress the population at x 1 (i.e., Âðx 1 Þ must be > 0).In turn, the algorithm seeks to find ÂðxÞ-the "minimal" investment needed to ensure that the speed does not exceed v in any location.Therefore, if exactly to ñðx 1 þ DxÞ within Δt time units.
Specifically, if the dynamics are time-discrete (Eq (2)), ÂðxÞ is such that, despite the possibility of having n 1 ðxÞ > ñðx þ DxÞ, the birth rate of a population with n 1 (x) individuals is the same as the natural birth rate of a population with ñðx þ DxÞ individuals: if n 1 ðxÞ > ñðx þ DxÞ, and ÂðxÞ ¼ 0 otherwise.In other words, ÂðxÞ is such that it reduces the effective population size from n 1 (x) to ñðx þ DxÞ.In particular, considering the case study of the spongy moth, where b (n,A) is given by Eqs (6-8), Eq (13) implies that ÂðxÞ is given by the following formula: where we denote n 1 = n 1 (x) and ñ ¼ ñðx þ DxÞ to simplify the notations.In turn, if the dynamics are time-continuous (Eq (1)), the rate at which n(x) declines due to the treatment is given by R(n(x),A(x)).Therefore, the condition that ÂðxÞ reduces n(x) from n 1 (x) to ñðx þ DxÞ within Δt time units can be written as: if n 1 ðxÞ > ñðx þ DxÞ, and ÂðxÞ ¼ 0 otherwise, where R -1 is the inverse function of R with respect to n (i.e., R(R -1 (n,A), A) = n).In the general case, ÂðxÞ can be found numerically for any given x from Eq (15) via the Newton-Raphson method [61] since R(n,A) is monotone increasing with A. In some special cases, however, a closed-form solution to Â exists, which could simplify the procedure.Particularly, if dynamics are time-continuous and R(n,A) = βAn α (Eq (5)) with 0�α<1, then R -1 (n, A) = (βA) -1 n -α , and it follows from Eq (14) that Stop or reverse the spread, v�0.Consider the case in which v < 0: The algorithm needs to find Â that reverses the direction of the front's propagation (i.e., the invasive species is being back-propagated at a speed -v).In that case, within Δt time units, n 1 (x) needs to decline within Δt time units to become � ñðx À DxÞ, where Δx is now defined as the distance traveled at a speed -v within Δt time units: -v = Δx/Δt.Therefore, if v < 0, the condition for ÂðxÞ > 0 is that n 1 ðxÞ > ñðx À DxÞ.Otherwise, everything else is the same here as in the algorithm for the case v > 0 except that ñðx À DxÞ replaces ñðx þ DxÞ in Eqs (13)(14)(15)(16).Finally, if v = 0, the condition for Â > 0 is that n 1 ðxÞ > ñðxÞ, and accordingly, ñðxÞ replaces ñðx þ DxÞ in Eqs (13)(14)(15)(16).
Algorithm: Finding n from ñ and Âð"find nð ÂÞ) For a given Â, which has been calculated from a given front ñ using "find ÂðñÞ," it is guaranteed that the propagation speed of nðx; 0Þ ¼ ñðxÞ at t = 0 is �v in all locations.(The speed equals v where ÂðxÞ > 0 but could be lower where ÂðxÞ ¼ 0).Namely, if we allow nðx; 0Þ ¼ ñðxÞ to evolve according to Eq (1) with Aðx; tÞ ¼ Âðx þ vtÞ, its shape may change because some regions of n will move leftward at a speed v while some regions will move slower.
We would like to find nðxÞ -the front shape that evolves in the long run and has a propagation speed of exactly v if treatment is given by Â, i.e., nðx; tÞ ¼ nðx þ vtÞ as t ! 1 if n(x,t) evolves according to Eq (1) with Aðx; tÞ ¼ Âðx þ vtÞ.To find n, the algorithm begins with nðx; 0Þ ¼ ñ and simulates the dynamics of n with A ¼ Âðx þ vtÞ, until n converges to a front shape that moves at a rate v but does not change its shape otherwise, i.e., nðx; tÞ ¼ nðx þ vtÞ.

Algorithm for non-stationary, discrete-time dynamics
According to the spongy moth population model described in Eqs (2,(6)(7)(8), large values of r may lead to chaotic dynamics of the untreated population.Therefore, one cannot expect that the population density in the untreated area-where x is large-will approach its carrying capacity.In such cases, the algorithm sets a location x 1 beyond which the population is not treated, and the natural dynamics of the population are simulated for x > x 1 .This way, the algorithm still converges to some front shape and some treatment shape within the range x < x 1 .In the present study, however, we restrict attention to cases where r is sufficiently small and the population converges to its carrying capacity in the non-treated area, while future studies could use the algorithm to analyze a broader set of problems.

Robustness check
The algorithm finds n opt by performing a local search in the space of all possible front shapes: it examines local modifications of nðxÞ until it reaches a configuration for which no local modification exhibits an improvement.Therefore, the solution n opt and its associated A opt correspond to a local minimum in the space of all front shapes.This means that it is possible that n opt and A opt are only approximations of the global optimum.To further explore whether there exists a better choice of n and Â, we repeat the algorithm with various different choices of the initial front shape ñ, and we verify that it converges to the same n opt and A opt .

Determining the optimal propagation speed
The algorithm we described here finds the optimal treatment, A opt (x), its ACT, and the corresponding long-term front shape, n opt (x), for a given value of v.By using this algorithm to find A opt and the ACT associated with it for various values of v, we can examine the interplay between the propagation speed and the cost of treatment, and find the optimal target speed v. Specifically, consider the net present value (NPV), given by minus the cost of treatment and the cost of damage over time: where the constant C is the annual cost per unit of an infested area due to the establishment of the invasive species there (we measure the infested area relative to the infested area at t = 0).Neglecting the short-term dynamics and assuming that Aðx; tÞ ¼ Âðx þ vtÞ and nðx; tÞ ¼ nðx þ vtÞ, it follows that where ACT*(v) is the ACT if Â ¼ A opt for a given target speed v (see a similar analysis in [44]).
In turn, we would like to find the value of v that maximizes the NPV.Note that −1 < v � v 0 , where v 0 is the speed at which the population front propagates if A = 0. Therefore, the NPV is maximized in one of the following three cases: (a) no treatment, i.e., A = 0 and v = v 0 (optimal if the NPV always increases with v); (b) eradicate as fast as possible, i.e., v !-1 (optimal if the NPV always decreases with v); or (c) control the species at some intermediate value of v (where the NPV has a local maximum).Specifically, it follows from Eq (18) that Therefore, the optimal strategy is to abandon treatment (A = 0 and v = v 0 ) if dACT*(v)/dv < C/δ for all v and eradicate as fast as possible if dACT*(v)/dv > C/δ for all v. Otherwise, the optimal value of v is where and On the other hand, if ACT*(v) in convex for some v = v 1 , it is more efficient to let the species propagate at a speed greater than v 1 for some time and then at a speed lower than v 1 for some time rather than letting it spread at a speed v 1 at all times.Note that Eq (20) has also been suggested also by [44], where here we also show how to calculate ACT*(v) from a more detailed population dynamical model.

Optimizing the treatment improves its cost-efficiency significantly
The algorithm finds the optimal treatment, at least approximately, for a wide variety of birth and death functions, as well as for various functional forms of the invasive species' response to the treatment.Our results for several cases of interest show that optimizing the spatial distribution of the treatment may lead to significant improvement in cost-efficiency.Specifically, the optimal treatment may be significantly more cost-effective than a treatment designed to stop the front while keeping its original shape or keeping a shape with a gradual linear decline.In particular, slowing the front but keeping its shape similar to that of an untreated population would result in an annual cost of treatment (ACT) higher by tens or hundreds of percentages than the ACT of the optimal treatment that slows the population's spread to the same speed.This result appears in the analysis of the general model (compare

The algorithm is robust
Our robustness check reveals that the algorithm is robust: it provides consistent results that are not sensitive to changes in initial conditions.Specifically, the algorithm converges to the same n opt and A opt for various choices of the initial front shape (Figs 2 and 3).This implies that the algorithm is robust and that it converges to either the optimum or to some approximation of the optimum that would be difficult to improve further.In particular, for the general model (Eqs 1, 3-5), the algorithm begins in one simulation with a front shape that emerges when no treatment is applied (  , where b and d are given by Eq (3) with r = 2 year -1 , k =2, and γ = 1 year -1 ; G is given by Eq (4) with σ = 1 km; and R given by Eq (5) with β 1.25 USD -1 year -1 and α = 0.2.
https://doi.org/10.1371/journal.pcbi.1011996.g002also in that case, the algorithm is robust and resulted in the same n opt and A opt for both initial choices of the front shape.

The optimal spatial distribution of the treatment largely depends on the population's responses to the treatment
The analysis of the general model (Eqs (1, 3-5)) shows that, if the marginal cost per removal of one unit of n is independent of n (α = 0), the optimal strategy is to treat only those areas where the treatment reduces the population density to zero (A(x) > 0 only if n(x) = 0; Fig 4 -bottom row).However, if treatment becomes less efficient as n decreases, as occurs in the general model (Eq (5)) with 0<α�1, the optimal treatment dictates treating in certain areas where n(x) > 0 (Figs 2 and 4 -top row).This is because, if α > 0, it is more cost-effective to reduce the population in locations where its density is higher, and treatment in areas where n(x) is larger reduces propagule pressure in areas where n(x) is close to zero.This result is consistent with the result of [21], who examined the special case in which v = 0 and considered a similar (though not identical) model.They showed that, if α = 0, treatment is needed only in locations x where n(x) = 0 (i.e., only a "barrier zone" is needed); but if treatment efficiency is proportional to the population size (α = 1), treatment is needed in some areas where n(x) > 0 (i.e., optimal treatment dictated treating in some "suppression zone").Our result thus extends this result to cases where the manager  ACT* decreases as v increases, until v equals the natural speed at which the species propagate without treatment v = v 0 �50, where ACT* = 0.The ACT* is also higher when α is larger.The sub-panels demonstrate the shape of the optimal treatment profile, A opt , and the population density, n opt , as a function of the location, x.The three sub-panels on the top row show results for α = 0.6, and those on the bottom row for α = 0.Each sub-panel shows the optimal solution for a given target speed v, where both the treatment and the population density continuously move leftward at that speed (Eqs (9,11)).The sub-panels demonstrate that, if α = 0, treatment is applied only where n(x) = 0, whereas if α = 0.6, treatment is distributed over broader areas, including some areas where n(x) > 0. If v < 0 (left sub-plots), a large concentration of treatment is peaked at the border between (a) the region where the population density is large and (b) the region slows the spread (v > 0) or reverses the spread (v > 0) and to cases where 0 < α < 1.In turn, Baker and Bode (2016) also examined optimal distribution of treatment for stopping the spread in a case where treatment efficiency is proportional to n.They found that, similar to our results, treatment should decline gradually in the area that needs to be protected from the species.In turn, the analysis of the spongy moth population case study, in which treatment is implemented via mating disruption (Eqs (2,4,(6)(7)(8)), reveals that the optimal strategy is to use mating disruption over a wide area that includes the regions where the population is below the Allee threshold but also a region where the population is above that threshold (Figs 3 and 5).

To reverse the spread, the manager should apply the treatment within a narrow area
The spatial distribution of the optimal treatment also depends on the target propagation speed, v.If v � 0, our results show that A opt (x) > 0 only if x is smaller than a certain threshold, and A varies smoothly within the range where it is positive (Figs 4 and 5).However, if v < 0, there is an apparent difference between the general model and the spongy moth model.Following the model of the spongy moth with mating disruption (Eqs (2, 4, 6-8)), A opt (x) remains a smooth function when v is negative (Fig 5).But following the general model (Eqs (1, 3-5)), the optimal strategy when v < 0 is to treat very aggressively within a narrow area that creates a border between (a) the area where the species' density is low and (b) the area where the species is established and its density is close to its carrying capacity (Fig 4).No additional treatment is applied in the area where the species' density is high, and some treatment is applied in the area where its density is low in order to reduce it further to zero.This treatment profile moves to the right alongside with the population front at a speed -v (Eqs (9,11)).This implies that the optimal strategy for reversing the spread is to "push" the invasive species backward from the area where its density is low toward the area where its density is high.

Tradeoff exists between the target propagation speed and the annual cost of treatment
The annual cost of treatment depends largely on the target propagation speed, v, and on the response of the species to the treatment.Our results reveal a tradeoff between slowing the spread and the annual cost of the optimal treatment, ACT* (v): In both the general model and the spongy moth model, slowing the spread more aggressively necessitates higher investments (i.e., lower v results in a higher ACT*; Figs 4 and 5).This tradeoff has been suggested by Sharov & Liebhold (1988), who examined specific ACTs that were taken as given (without optimizing a front).In addition, population's response to the treatment also affects the ACT*.In the general model, the larger α is, the more expensive it is to slow the population's spread (Fig 4).In the gypsy moth population model, a larger kλ 0 implies that the Allee effect appears at lower population densities, which makes it more difficult to slow the spread.Accordingly, the larger kλ 0 is, the more expensive it is to slow the population's spread and the faster ATC increases as v decreases (Fig 5).Moreover, in both the spongy moth and the general model with 0 � α < 1, stopping the spread (v = 0) and even reversing the spread (v < 0) can be achieved by investing a finite cost, and this may be optimal if C/δ is sufficiently large.In turn, note that in the general where the rest of the treatment is applied.All the parameters except α and v are the same in all the panels: Eq (1) is considered, where b and d are given by Eq (3) with r = 2 year -1 , k = 2, and γ = 1 year -1 ; G is given by Eq (4) with σ = 25 km -1 ; and R given by Eq (5) with β = 1/(1α) USD -1 year -1 .The raw data with the results of all the simulations for each α and v can be found on Dryad [60].https://doi.org/10.1371/journal.pcbi.1011996.g004ACT* decreases as v increases, and reaches zero when v equals the species' natural propagation speed without treatment (v = v 0 �12 km year −1 ).Higher kλ 0 values, indicating a lower Allee threshold relative to carrying capacity, resulting in an increased ACT*, where ACT* is measured in units of thousands of USD per one-kilometer strip of land.The dashed black lines show, for each line, where the marginal cost of slowing the spread (slope of ACT*) equals the marginal benefit of slowing the spread by one kilometer over a one-kilometer strip, estimated as 16K USD per year by [62].In particular, the optimal speed according to this estimation is v = -2.4km/year if kλ 0 = 100; v = 4.8 km/year if kλ 0 = 300; and v = 9.6 km/year if kλ 0 = 1000.In turn, the sub-panels demonstrate the optimal treatment profile (A opt ), and population density (n opt ) across locations (x in units of σ = 10 km).Each sub-panel shows the optimal solution for a specific target speed v, with both treatment and population density advancing leftward at this speed (Eqs (9,11)).All the parameters except α and v are the same in all the model with α = 1, it becomes impossible to stop the spread because ACT* approaches infinity as α approaches one from below.

The optimal target propagation speed increases as the cost of damage decreases
The optimal choice of v, which maximizes the NPV (Eq (18)), can be found from the function ACT* (v) (Fig 4) together with the ratio between the cost of damage due to the species per unit area, C, and the discount rate, δ (Eq (20)).In particular, if C/δ is below a certain threshold, such that the slope of ACT* (v) is greater than C/δ even near v = v 0 , then abandoning the treatment (A = 0 and v = v 0 ) is optimal.But if C/δ is above that threshold, the optimal choice of v is where ACT* (v) is concave and its slope equals C/δ (Eq (20)) [44].In turn, Figs 4 and 5 demonstrate that ACT* (v) (Eq (21)) is a concave function for a wide range of treatment methods and propagation speeds, which implies that the marginal cost of slowing the spread by additional kilometer per year increases as v decreases.Therefore, the optimal choice of v is smaller if C is larger.
In particular, [62] have conducted a thorough study to valuate the damage incurred by the spongy moth.They found that slowing the spread by one kilometer annually across a one-kilometer strip yields a benefit of 16K USD per year.(Specifically, Leuschner et al. (1996) found that the benefit of slowing the spread by 4 kilometers over a strip of about 800 Kilometers in length is worth 795M USD total or 50.9MUSD per year, and therefore, the annual benefit per kilometer square is 50.9M/(4800) = 16 K USD.)Our results show that, for parameter values used in Fig 5, the optimal target speed is v = -2.4km/year if kλ 0 = 100; v = 4.8 km/year if kλ 0 = 250; and v = 9.6 km/year if kλ 0 = 1000.These results, particularly those where kλ 0 = 100 or kλ 0 = 300, are in line with previous conclusions [63].These results imply that the strong Allee effect is a major reason why it could be cost-effective to slow the spread of the spongy moth.

Limitations & future directions
Our study has several limitations, which could be addressed in future studies.First, the algorithm proposed here is tailored to optimizing the containment of an invasive species, which is only one component in a more comprehensive control program that could address prevention surveillance [64].Second, the algorithm finds the spatial distribution of the population and the treatment that optimizes the annual costs in the long run.Additional research is needed to determine the most cost-effective management strategies during the initial phase-until the population front reaches its asymptotic shape.Third, the algorithm searches for the optimal front shape locally, and therefore, there is no guarantee that it always finds the global optimum.Fourth, the model does not incorporate the construction of physical barriers that prevent dispersal, such as fences and dams that restrict the movement of mammals [65] and aquatic species [18], and future research is needed to examine how to incorporate such strategies [36].Fifth, we restricted attention to non-chaotic parameter regimes of the discrete-time dynamics.The algorithm can work in the chaotic regimes, but future studies are needed to analyze these cases.Finally, the dynamics in our model are deterministic, whereas real-world invasions exhibit stochastic dynamics and uncertainty.In particular, invasion could appear in more distant location and the locations where new populations appear might not be predictable.Future studies could panels: Eq (2) is considered, where b is given by Eqs (6-8) with r = 2 and a = 0.08 USD -1 , and G is given by Eq (4).The raw data with the results of all the simulations for each kλ 0 and v can be found Dryad [60].https://doi.org/10.1371/journal.pcbi.1011996.g005

FindingFig 1 .
Fig 1. Flow diagram summarizing the algorithm for finding the optimal treatment.https://doi.org/10.1371/journal.pcbi.1011996.g001 Fig 2A where ACT = 28.58 and Fig 2D where ACT = 8.58) as well as in the analysis of the spongy moth case study (compare Fig 3A in which ACT = $278.2Kand Fig 3D in which ACT = $56.29$,where ACT is measured in units of thousands of USD per one-kilometer strip of land over which the front propagates).A similar improvement in cost-efficiency is obtained if one uses the optimal treatment instead of a treatment that maintains a front with a shape of a straight line.This efficiency gain is evident in both the general model (compare Fig 2E where ACT = 25.82 and Fig 2H where ACT = 8.52) and the spongy moth case study (compare Fig 3A where ACT = $197.3Kand Fig 3D where ACT = $56.15K).
Fig 2A), and in another simulation with a piecewise-linear front shape (Fig 2E), and it resulted in the same n opt and A opt in both cases (Fig 2D and 2H).We performed the same analysis with the case study of the spongy moth (Eqs 2, 4, 6-8; Fig 3), and

Fig 2 .
Fig 2.The algorithm finds the optimal treatment, which slows the population's propagation to a given target speed, v, while minimizing the annual cost of treatment (ACT) (simple model).(A) The algorithm begins with a population front, ñ, that has evolved naturally when the species has propagated at a speed v 0 (v 0 > v; blue line).The algorithm then finds the treatment function Â for which the speed of the leftward movement of the front does not exceed v in any location.(B-D) The algorithm finds new front shapes, as well as the treatment functions that hold these fronts propagating leftward at a speed v.In each iteration, the algorithm changes ñ and Â to those for which the ACT is lower.(E-H) The same algorithm as in (A-D) is executed, only here it begins with a piecewise-linear population front (E).In both simulation no. 1 (A-D) and no. 2 (E-H) of the algorithm, ñ and Â converge to a similar shape ((D) is similar to (H)), and we denote ñ and Â of this final outcome of the algorithm (shown in (D) and (H)) as n opt and A opt , respectively.Parameters are the same in all panels: The target speed is v = 10 km/year, and the dynamics of n follow Eq(1), where b and d are given by Eq (3) with r = 2 year -1 , k =2, and γ = 1 year -1 ; G is given by Eq (4) with σ = 1 km; and R given by Eq (5) with β 1.25 USD -1 year -1 and α = 0.2.

Fig 3 .
Fig 3.The algorithm finds the optimal treatment, which slows the population's propagation to a given target speed, v, while minimizing the annual cost of treatment (ACT) (spongy moth population model).The description of the panels is similar to that in Fig 2. (A-D) The algorithm begins with a population front ñ that has evolved naturally (A); it finds front shapes that could be slowed with lower annual costs (B-D), until reaching a population front ñ ¼ n opt for which ACT is minimized.(E-H) The same algorithm as in (A-D) is executed, only here it begins with a piecewise-linear population front (E).As in Fig 2, in both simulations no. 1 (A-D) and no. 2 (E-H) of the algorithm, the population front and the corresponding treatment function converge to a similar shape ((D) is similar to (H)).Units: Distance (x) is shown in units of σ = 10 km: population size (ñ) is in units of its carrying capacity; treatment cost ( Â) is in USD per hectare per year; and ACT is given in thousands of USD per onekilometer strip of land.Parameters are the same in all panels: The target speed is v = 220 m/year, and the dynamics of n follow Eq (2), where b is given by Eqs (6-8) with r = 2, kλ 0 = 100, and a = 0.08 USD -1 .https://doi.org/10.1371/journal.pcbi.1011996.g003

Fig 4 .
Fig 4. The annual cost of treatment decreases with the target speed v and increases with α (simple model).The main (middle) panel shows the annual cost of treatment associated with the optimal treatment, ACT*, as a function of the target propagation speed of the population, v, for various choices of α (Eq (5)).ACT* decreases as v increases, until v equals the natural speed at which the species propagate without treatment v = v 0 �50, where ACT* = 0.The ACT* is also higher when α is larger.The sub-panels demonstrate the shape of the optimal treatment profile, A opt , and the population density, n opt , as a function of the location, x.The three sub-panels on the top row show results for α = 0.6, and those on the bottom row for α = 0.Each sub-panel shows the optimal solution for a given target speed v, where both the treatment and the population density continuously move leftward at that speed (Eqs(9,11)).The sub-panels demonstrate that, if α = 0, treatment is applied only where n(x) = 0, whereas if α = 0.6, treatment is distributed over broader areas, including some areas where n(x) > 0. If v < 0 (left sub-plots), a large concentration of treatment is peaked at the border between (a) the region where the population density is large and (b) the region

Fig 5 .
Fig 5.The annual cost of treatment decreases with the target speed v and increases with kλ 0 (spongy moth population model).As in Fig 4, the main (middle) panel shows the annual cost of treatment associated with the optimal treatment, ACT*, as a function of the target propagation speed v for three choices of kλ 0 .ACT* decreases as v increases, and reaches zero when v equals the species' natural propagation speed without treatment (v = v 0 �12 km year −1 ).Higher kλ 0 values, indicating a lower Allee threshold relative to carrying capacity, resulting in an increased ACT*, where ACT* is measured in units of thousands of USD per one-kilometer strip of land.The dashed black lines show, for each line, where the marginal cost of slowing the spread (slope of ACT*) equals the marginal benefit of slowing the spread by one kilometer over a one-kilometer strip, estimated as 16K USD per year by[62].In particular, the optimal speed according to this estimation is v = -2.4km/year if kλ 0 = 100; v = 4.8 km/year if kλ 0 = 300; and v = 9.6 km/year if kλ 0 = 1000.In turn, the sub-panels demonstrate the optimal treatment profile (A opt ), and population density (n opt ) across locations (x in units of σ = 10 km).Each sub-panel shows the optimal solution for a specific target speed v, with both treatment and population density advancing leftward at this speed (Eqs(9,11)).All the parameters except α and v are the same in all the