Spatially explicit control of invasive species using a reaction–diffusion model

Ecological Modelling 337 (2016) 15–24

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Spatially explicit control of invasive species using a reaction–diffusion model Mathieu Bonneau a,∗ , Fred A. Johnson b , Christina M. Romagosa a a b

Department of Wildlife Ecology and Conservation, University of Florida, 110 Newins-Ziegler Hall, P.O. Box 110430, Gainesville, FL 32611-0430, USA Wetland & Aquatic Research Center, U.S. Geological Survey, 7920 NW 71 Street, Gainesville, FL 32653, USA

a r t i c l e

i n f o

Article history: Received 28 August 2015 Received in revised form 21 May 2016 Accepted 24 May 2016 Keywords: Allocation Burmese pythons Control Invasive species Reaction–diffusion model Simulation Spatial distribution Ecological modeling

a b s t r a c t Invasive species, which can be responsible for severe economic and environmental damages, must often be managed over a wide area with limited resources, and the optimal allocation of effort in space and time can be challenging. If the spatial range of the invasive species is large, control actions might be applied only on some parcels of land, for example because of property type, accessibility, or limited human resources. Selecting the locations for control is critical and can significantly impact management efficiency. To help make decisions concerning the spatial allocation of control actions, we propose a simulation based approach, where the spatial distribution of the invader is approximated by a reaction–diffusion model. We extend the classic Fisher equation to incorporate the effect of control both in the diffusion and local growth of the invader. The modified reaction–diffusion model that we propose accounts for the effect of control, not only on the controlled locations, but on neighboring locations, which are based on the theoretical speed of the invasion front. Based on simulated examples, we show the superiority of our model compared to the state-of-the-art approach. We illustrate the use of this model for the management of Burmese pythons in the Everglades (Florida, USA). Thanks to the generality of the modified reaction–diffusion model, this framework is potentially suitable for a wide class of management problems and provides a tool for managers to predict the effects of different management strategies. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Invasive species (IS) are responsible for severe ecological damages and economic losses (Olson, 2006; Dehnen-Schmutz et al., 2007; Pysek and Richardson, 2010), and their management (e.g. control and/or eradication) is a challenging problem due to uncertainty about the species biological characteristics and limited resources with which to undertake control. Although the most cost-effective method for minimizing damage is one that prevents the arrival of non-indigenous species (Lodge et al., 2006; Keller et al., 2007), the implementation of such a strategy is not obvious when the potential invasiveness of a species is unknown and the species has commercial value (Keller and Springborn, 2013). As a consequence, the management of IS often starts long after their introduction, when ecological damages are visible and the species is already well established. In such a situation, the area occupied by the species can be large, and managers have to act with

∗ Corresponding author. Tel.: +1 3522643504. E-mail addresses: mbonneau@ufl.edu (M. Bonneau), [email protected] (F.A. Johnson), cmromagosa@ufl.edu (C.M. Romagosa). http://dx.doi.org/10.1016/j.ecolmodel.2016.05.013 0304-3800/© 2016 Elsevier B.V. All rights reserved.

limited resources in order to slow the invasion and/or eradicate the IS. Managers typically face the decision of where, when, how much and/or how often control has to be implemented. The optimal management of IS, reviewed in Epanchin-Niell and Hastings (2010), is discussed in a growing number of articles, in which biological models are coupled with efficient optimization routines (Mehta et al., 2007; Bogich et al., 2008; Hauser et al., 2009; Carrasco et al., 2010; Giljohann et al., 2011). But in many practical situations, the decision problem is spatially explicit by nature (Meier et al., 2014). Therefore, the problem is to decide where in the entire management area to allocate control effort. As the complexity of this optimization problem grows with the number of control units in the management area, finding the optimal management strategy among all possible alternatives becomes computationally intractable. Another important issue arises when there must be an accounting of the temporal dynamics of the invader population. Then, managers have to optimize their choices, not only for the current time step, but for possible scenarios in the future. One approach consists in using heuristic methods to find a suboptimal management strategy (Schapaugh and Tyre, 2012; Nicol and Chadès, 2011). Although these methods should be favored first, heuristic methods can be hard to use in the most general case

16

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

because of their computational complexity. And in some cases, even if the management area remains large, the areas where control can be implemented can be small. In this case, the challenge is not so much to compute the optimal control strategy, but to predict and compare the efficiency of the different options available to the manager. Then a natural approach is simulation based (Higgins et al., 2000; Grevstad, 2005; Provencher et al., 2007; Frid et al., 2013), which requires a dynamic, spatially explicit model of the invader population that takes into account the effects of control. Reaction–diffusion (RD) models are general, spatially explicit models that use a partial differential equation to describe changes in population density over time and space. Although some limitations of these models are known (Svenning et al., 2014), they have been widely applied in ecology (Holmes et al., 1994; Cantrell and Cosner, 1996, 2003; Shigesada and Kawasaki, 1997; Bogich et al., 2008; Kaiser and Burnett, 2010; Acevedo et al., 2012; Leroux et al., 2013) and provide a good basis to describe the pattern of diffusion (Andow et al., 1990; Hastings et al., 2005). The basis of this model is described by Fisher’s equation with logistic growth (Fisher, 1937; Skellam, 1951): 2

2

n ∂(n) ∂ (n) ∂ (n) + ) + n(1 − ), = D( 2 2 K ∂x ∂y ∂t   







(1)

Local growth

Spread

where n is the species density at a given location (i.e. number of individuals per square kilometer), D is the diffusion coefficient (i.e. the rate at which individuals disperse),  is the intrinsic rate of growth (i.e. the birth rate minus the death rate) and K is the socalled carrying capacity (i.e. the maximal density per location). The left term is responsible for the spread or diffusion of the species, and the right term for local population growth (here a logistic growth function). Classic RD models make the assumption that the invader is introduced at a given point in space and then diffuses in a homogeneous landscape in every direction with the same probability, such that movement can be characterized as a random walk. In an infinitely large and suitable landscape, the pattern of invasion is described by a growing circle centered on the point of first introduction. Although simple, this colonization pattern is quite general and can be used to model invasions in practice with reasonable computational complexity. But to be used to explore the control of IS, the Fisher equation (1) has to be modified in a way that accounts for the effect of control, not only on the local growth of the population, but on the diffusion process as well. To account for the effect of control on population growth, a useful approach is based on the Schaefer model (Schaefer, 1957), extensively used in the field of fisheries economics. See for example (Clark, 1990) for a detailed analysis of the Schaefer model. This model is defined without any spatial component: if ˇ is the control mortality (or fishing mortality or harvest rate) of the control method (or any removal method), the change in the population density is defined as follows: ∂(n)/∂t = n(1 − n/K) − ˇn. Several authors (Oruganti et al., 2002; Neubert, 2003; Kurata and Shi, 2008) extended this last equation for two dimensions by adding the diffusion term, as it appears in the original Fisher’s equation:

∂(n) =D ∂t





2

2

∂ (n) ∂ (n) + ∂x2 ∂ y2



Spread



 

n +n  1− K

 







−ˇ .



Local growth

This modification is known as the Constant–Effort Harvesting model (CEH). For example, Neubert (2003) used this model to discuss the optimal locations of marine reserves, as well as the optimal level of harvesting. As one can see, the CEH model makes the assumption that the diffusion term is not affected by control (i.e. the spread term is unchanged), or equivalently that individuals cannot

be removed while dispersing. For invasive species this assumption is debatable, as far as they are generally controlled regularly over the year, thus both in-situ and during dispersal. We propose to modify the diffusion term to explicitly account for the fact that individuals can be removed when they are diffusing. This results in a CEH model with an additional linear control mortality, where the density at any point is explicitly affected by the probability that immigrants have been removed. We first compared the prediction of the CEH model and the modified RD model using a simulated data set. We then illustrated the use of this model for management of the Burmese python (Python molurus bivitattus) in the Everglades, Florida. 2. Materials and methods 2.1. Accounting for the effect of control actions This section present the modified CEH model with a total linear removal rate. A complete proof of this new model is available in Appendix S1. We assume that J control actions are available to the manager. The control actions are associated with the yearly control mortality ˛T1 , . . ., ˛TJ and weights rT1 (x), . . ., rTJ (x). And let’s denote CTj ⊂ R2 , j = 1. . .J, the spatial domain where control action Tj is applied. The different control actions can be for example, physical, chemical or biological. Then, the yearly control mortality is the added mortality created by the control action over a year (i.e. a reproduction cycle). Without loss of generality, we suppose that any point in space can be controlled by at most one control J action, which implies that ∩j=1 CTj = ∅. In practice, a control action Tj can result from different control actions and the control mortality has to be computed accordingly. But in the model, Tj always appears as a single control action. As we consider the spatial problem of management, x denotes a spatial coordinate: x = (x, y) ∈  ⊆ R2 , where  is the entire spatial domain or management area. Then x ∈ CTj means that location x is controlled with control action Tj . Formally, rTj (x) is the probability that an immigrant to location x crossed an area controlled by action Tj before arriving in x. And thus, ˛Tj × rTj (x) is the probability that the individual will be removed before reaching location x. In practice, obtaining an analytic value of rTj (x) is not straightforward and an approximation is needed. We propose an intuitive way to compute rTj (x) in the next section. Finally, let ˇ(x) denote the yearly control mortality of the control action applied at location x. If x is not controlled, then ˇ(x) = 0 and if x is controlled with action Tj , ˇ(x) = ˛Tj . In such a situation the population dynamic at location x is described as follows:

∂(n(x, t)) ∂ (n(x, t)) ∂ (n(x, t)) + ) − n(x, t) ˛Tj rTj (x) = D( ∂t ∂x2 ∂ y2 2

2

J

j=1

+ n(x, t)((1 −

n(x, t) ) − ˇ(x)) K

(2)

This modified RD model corresponds to the CEH model with an additional linear control mortality term and will be further denoted the LCM model. The only difference from the CEH model is the addiJ tional term −n(x, t) j=1 ˛Tj rTj (x), which appears in the reaction term. As expected, not only the growth of the population is influenced by the control method at this location but also by control in its neighborhood, as control of dispersing individuals naturally decreases the number of individuals that are immigrating. One can

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

note that the CEH and LCM model are equivalent when no control actions are implemented. In terms of management this has several implications, the most obvious being that the estimation of population density will be lower when using this modified equation compared to that with the CEH model. Another implication is that the locations of the control actions are critical to overall control efficiency. 2.2. Implementation 2.2.1. Approximation of the probability that an immigrant crossed a controlled area We first present the approximation of the probabilities rTj (x). For any location xEnd ∈ , rTj (xEnd ) is the probability that an immigrant from xStart to location xEnd will cross a controlled location before arriving in xEnd . Using Bayes’ rule we have: rTj (xEnd ) =



P (xStart → xEnd )

xStart ∈ 



× P xStart → xEnd ⊂ CTj | xStart → xEnd . where P (xStart → xEnd in xStart ) is the probability that individuals

move to xEnd and P xStart → xEnd ⊂ CTj | xStart → xEnd is the probability that individuals crossed a controlled location while moving to xEnd . The later probability is related to all the possible paths that can be used by the individual to reach the two locations, denoted PxStart →xEnd , their probabilities of occurrence P PxStart →xEnd and either or not the path crosses a control area. The previous equation can be rewritten as follows: rTj (xEnd ) =



P (xStart → xEnd )

xStart ∈ 

×





P PxStart →xEnd ıPx

Start →xEnd

,

PxStart →x

End

where ıPx →x is equal to one if the path crosses a controlled area Start End and zero otherwise. Unfortunately, computing the probability of moving from xStart to xEnd , as well as counting all the possible paths between xStart and xEnd , is not straightforward and is beyond the scope of the article. We thus propose to use the following approximation. Let’s denote dm the average yearly movement rate and thus during a year, individuals can only migrate to locations not further than dm kilometers. Then, let’s denote N(x) the disc centered in x of radius dm . In other words, N(x) is the domain containing all the possible locations where an immigrant to x can come from. It follows that N(xStart ) contains all the possible migration locations from xStart . Let further assume that individuals from xStart can immigrate to any points inside N(xStart ) with the same probability, which implies that: P (xStart

1 → xEnd ) = . Area (N(xStart ))

This is equivalent to assume that an individual can be located anywhere from its last year location in a range of dm kilometers with equal probability. This approximation seems to be reasonable when the habitat is supposed to be homogeneous over the landscape. Finally, note that the area of the region N(xStart ) does not depend on the position of xStart and thus we can use Area (N(xEnd )) instead.  Next,  = ıPx →x has to be approximated. For P xStart →xEnd

Start

End

this we propose to simplify the computation in assuming that a path crosses a controlled location if and only if the starting point is itself controlled and thus  is equal to one if xStart is controlled (i.e.

17

xStart ∈ CTj ) and zeros otherwise. Finally, we obtain the following approximation for rTj (x):



rTj (xEnd ) =

Area CTj ∩ N(xEnd ) Area (N(xStart ))

,

One can note that the simplification that we propose clearly underestimates the value of  as far as all the paths that crossed a control region but not start on a controlled location are not included. As a result rTj (xEnd ) is also underestimated, as well as the effect of control on dispersion. A more rigorous approximation might be found using the theory of Brownian motion. However, we will show later in this article that this first approximation already enables us to improve the prediction over the CEH model. Fig. S1 provides an example of the value of the approximation of rTj (x).

2.2.2. Boundary condition and numerical solution In practice, the domain is obviously finite and delimited by frontiers, and specific conditions at these frontiers have to be defined. We consider that the frontiers of the domain are the natural barriers present in the landscape. More precisely, let’s consider that the frontiers of the domain are the boundaries between suitable and unsuitable habitat. For example, if we study the spread of a terrestrial invasive species on an island, the natural barrier is the island coast. Mathematically, these frontiers are called boundaries and the invader density is assumed to follow a particular equation at these locations, which is called the boundary condition. Here we describe two types of boundary conditions: (i) the Dirichlet boundary condition and (ii) the Neumann boundary condition. The Dirichlet boundary condition is used for porous barriers like unsuitable habitat. It is appropriate when the species is allowed to cross the barrier, but can never return. In contrast, the Neumann boundary condition is used for non-porous habitat, in which crossing the boundary is not possible. The Dirichlet boundary condition is expressed as follows:

∀x ∈ ı, t ≥ 0,

n(x, t) = bD (x),

where ı is the domain boundary and bD is a real value function that returns the value of the density on this boundary. When bD is a constant fixed to 0, the boundary acts like a sink: the invader crosses the boundary and never returns. The Neumann boundary condition is expressed as follows:

∀x ∈ ı,

∂(n) (x) = bN (x), ∂v

where bN is a real value function that returns the value of the derivative of the density function on the boundary and v is the outward pointing unit normal vector to the boundary ı. When bN is also constant and fixed to 0, the boundary acts like a wall: the invader never crosses the boundary and chooses another direction. When the domain is bounded, the estimation of the population density has to be performed using numerical methods. Common and straightforward techniques are finite difference methods (Strikwerda, 2004), which aim to replace the derivatives in the differential equation by a finite difference term. As in Acevedo et al. (2012), we use the Crank–Nicholson method, known to be

18

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

unconditionally stable. Computing the density function of the population involves solving the following problem:

⎧ J 2 2 ⎪

⎪ ∂(n(x, t)) n(x, t) ∂ (n(x, t)) ∂ (n(x, t)) ⎪ ⎪ D( + ) − n(x, t) ˛Tj rTj (x) + n(x, t)((1 − ) − ˇ(x)), x ∈ , t > 0, = ⎨ K ∂t ∂x2 ∂ y2 j=1

⎪ n(x, t) = bD (x), ⎪ ⎪ ⎪ ⎩ 0

x ∈ ı, t ≥ 0,

n(x, 0) = n (x),

x ∈ .

The first two rows explain how the population density changes in the domain . The third row of the system explains the boundary condition. Here we use the Dirichlet condition but one could use the Neumann instead. One can also mix different boundary conditions on different parts of the domain boundary. The last row of the system is the initial condition. It gives the initial value of the population density at the beginning of the study or at the beginning of the invasion. More details on the implementation of the Crank–Nicholson method are available in Appendix S2 and a hypothetical example is provided in Appendix S3. We also provide Matlab codes of the Crank–Nicholson method in supplementary material. The simplicity of the RD model is obviously one of its advantages, but it can also be an important drawback for practical applications. Indeed, the invasion process is often complicated by the possibility of long-distance dispersal and secondary points of introduction. For the application that is considered in this article, these two features are not necessarily relevant, but may be in many other applications. To extend the use of this model, we propose in Appendix S4 a more general framework that allows these two extensions. For the rest of the article, we consider that both long-distance dispersal and secondary points of introduction do not apply. 2.3. Illustrations 2.3.1. Comparison with the CEH model on two hypothetical examples We compared the prediction of the CEH and the LCM model on the simulated spatial distribution of hypothetical species. We recall that a fundamental assumption of the RD model is that the species movement follows a random walk (Skellam, 1951). We thus used this framework to simulate a data set of the spatial distribution of a hypothetical species. We then used this data set to compare with the prediction of each RD model. Simulations of the data set. We simulated the spatial distribution of the species using a cellular automaton model on a square lattice composed of 100 × 100 cells of 1 km2 . We initialized the population with K individuals in 9 cells located in the center of the domain (see Fig. 1a). We used a constant carrying capacity of K = 100. We supposed that the species has a yearly reproduction cycle and considered the general case where it is diffusing during the year and then reproduces, which determines the end of the year. The diffusion of each individual is first simulated in allowing them to make nM moves according to a non-biased random walk. As we considered a two dimensional lattice, the possible moves are: moving to the cell located on the right, left, up or down of the individual current position. After a movement, if the individual is in a controlled cell, it is removed with a probability ˛/nM , where ˛ is the yearly control mortality. Finally, when the nM moves are performed, growth is in turn simulated using the following equation:

 

Nit+1 = Nit + Nit

 1−

Nit K



the species density is then used to create a data set of simulated species spatial distributions. We used this process to simulate 500 species spatial distributions. Finally, the predictions of the CEH and LCM models are compared to this data set. To compare each model in different situations, we used the cellular automaton model with various species and control characteristics. In order to examine variation in the diffusion and growth rate of the population, we considered two different set of parameters. The first set corresponds to a low-dispersal scenario: (, nM ) = (0.3, 5) and the second to a high-dispersal scenario: (, nM ) = (0.6, 15). To be able to compare the prediction when control is applied on a large or fine scale, we defined a single square controlled area composed of 121 cells (large scale) and a straight line composed of 19 cells starting at the bottom right corner of the square controlled area (fine scale, see Fig. 1b). For both of the low and high dispersal scenarios we considered the following yearly control mortality: ˛ = 0.2%, ˛ = 0.4% and ˛ = 0.6%. For low-dispersal scenarios, we simulated the population dynamics over 20 years and for high-dispersal scenarios we simulated the population dynamics over 10 years to keep the population away from the boundaries. Note that we implemented Dirichlet boundary conditions but the study area is large enough so that the population is not able to reach the edges of the domain. Thus, the boundary conditions do not impact the prediction. Parameters of the RD models and comparison criterion. For both RD models the control mortality, growth rate and carrying capacity are equal to the ones that we used for the cellular automaton model. For the estimation of the diffusion coefficient we proceed as follows. We first simulated the population dynamics using the cellular automaton model without control over 20 years. Each year, we computed the square root √ of the occupied area, denoted at and √ used the relation at = 2t D (see Eq. (7)) to estimate the diffusion coefficient D using a linear regression (Skellam, 1951). Once



− ˇi

,

(3)

where Nit is the density of individuals in cell i in year t and ˇi is equal to ˛ if the cell is controlled and 0 otherwise. This process is repeated until the time horizon is reached. The final spatial distribution of

Fig. 1. Configuration of the initial state of the population matrix and controlled cells for the hypothetical examples. The population is initialized with K individuals in 9 cells located in the center of the domain. The occupied (unoccupied) cells are colored in black (white). The controlled cells are colored in gray. Control is applied on a square area composed of 121 cells and on a straight line composed of 19 cells. Note that some initially occupied cells are also controlled (dark gray).

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

19

Fig. 2. (a) In black is the quarantine area where pythons are controlled. Here is an example with a quarantine area of h = 40 km, which is the maximal possible width of the quarantine area. (b) Predicted response of the pythons population in the Arthur R. Marshall Loxahatchee National Wildlife Refuge as a function of the time horizon T, control mortality ˛ and quarantine width h. The results are presented as a regression tree and a response of S means that the density in the refuge is predicted to be lower than 1 individual per km2 . The response is F otherwise.

the RD models are parameterized, each of them is used to predict the spatial distribution of the hypothetical species under control in 20 or 10 years depending on the dispersal scenario. Finally, we simply compared the prediction in terms of the average sum of the absolute difference between the map created by each of the RD models and the 500 simulated maps in the data set:

Error (CEH or LCM) 1

= 500 500

s=1

n c



abs N CEH or LCM (xi ) − NsSimu (xi )

 ,

(4)

i=1

where nc = 10,000 is the number of cells, NCEH (xi ), NLCM (xi ) and NsSimu (xi ) are the predicted densities in cell xi at the final step by the CEH model, the LCM model, and in the simulated map number s, respectively. Then, we compared the prediction of each model using the percentage of increased error produced by the CEH model, denoted pctLCM : pct LCM =

100 ∗ (Error(CEH) − Error(LCM)) Error(LCM)

In other words, it means that the prediction error of the CEH model is equal to Error(LCM)(1 + (pctLCM /100)). 2.3.2. Management of Burmese pythons Context. Burmese pythons are particularly well adapted to the available habitats in the Everglades where they can find a variety of native prey species, such as small mammals, birds, and alligators. Pythons are now recognized as widespread in the southern Everglades and are believed responsible for significant ecological damage (Dorcas et al., 2012). There are several lines of evidence that show pythons are progressively moving north and thus will become a threat to ecologically sensitive areas like the Arthur R. Marshall Loxahatchee National Wildlife Refuge (here after the refuge). Control of the Burmese python in the Everglades is particularly challenging because the habitat is highly unfavorable for control: a large proportion of the Everglades offers suitable habitat for pythons, but only a small portion of the Everglades is easily accessible by humans (e.g., roads and levees). In addition, the efficacy of control actions is limited by the detection probability (i.e. the probability that a human observer detects a python given

its presence), which is suspected to be around 0.01 (Dorcas and Willson, 2013). To date, control of pythons through direct removal from the wild is mostly conducted on roads and levees. These locations are likely to remain the only feasible locations for systematic control in the foreseeable future. To illustrate the model, we are looking here for the effect of increasing the control mortality on roads and levees only. We were interested if increasing the control mortality on these locations would be sufficient to delay or prevent the invasion of the refuge. See Fig. S2 for the location of the refuge. More details on the initial invasion of pythons in southern Florida can be found in (Willson et al., 2011). Increasing the control mortality on roads and levees. We assumed here that control can be undertaken only on public lands, located to the south of the refuge (see Fig. S2a). We assumed no improvement in the existing control method in areas other than roads and levees. Thus, we assumed a constant control mortality of ˛T1 = 0.01 in these areas. On the roads and levees bordering the refuge, we consider differing control mortalities: ˛T2 =0.01, 0.1, 0.3 and 0.5. The scenario ˛T1 = ˛T2 = 0.01 is close to the current situation and can be viewed as the baseline. The scenarios with ˛T2 > 0.01 represent possible improvements in control. Finally, we assumed that control starts in 2015, is undertaken every year, and then estimated the python population in the refuge until 2030. Quarantine area over the refuge We also provide an illustrative example in order to investigate potential scenarios that can decrease the python population. In this case, we considered that control is implemented around the refuge only, on a strip of various widths, denoted h. The strip can be viewed as a quarantine area to prevent the refuge invasion. We considered that control is implemented only south of the northern boundary of the refuge (see Fig. 2). We investigated several quarantine widths: h = 5, 10, 15, 20, 30 or 40 km; control mortalities: ˛ = 0.01, 0.1, 0.3; and time horizon: T = 2020, 2025, 2030.

2.3.3. RD model for Burmese pythons Grid. The domain of interest starts in extreme South Florida and ends just north of Lake Okeechobee (see Fig. S2). The south, east, and west boundaries are the Florida coasts, which are natural barriers to the spread of pythons. As a consequence, we implemented the Neumann boundary condition with bN ≡ 0 on these boundaries. The North boundary is determined more arbitrarily.

20

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

It was chosen far away so that the predicted number of pythons in this area remained small even for the longest time horizon (i.e. 2030). Thus either boundary condition will not have much effect on the population dynamic in the refuge. We implemented a Dirichlet boundary condition on the northern boundary with bD ≡ 0, which seems most appropriate. The domain is divided into square cells of 1km by 1km. As in the hypothetical example, x = y = 1. The density of pythons within the cell is calculated with the new RD model as the density in the center of the cell. As the density is expressed in number of pythons per square kilometer, it can be directly interpreted as the number of pythons within the cell. Then, the estimated number of pythons in the refuge is obtained as the sum of the number of pythons in every cell that is located in the refuge. We considered that a cell is controlled with control mortality ˛T2 either if its center is located at a maximal distance of 200 m from a road or levee, or if it is located on the border of the refuge. Then, all the remaining cells in the public areas are controlled with the control mortality ˛T1 = 0.01. Model parameters. To implement our model for pythons, we required estimates of the intrinsic rate of population growth, the diffusion coefficient, carrying capacity, and initial conditions. We relied on published information and the online database (EDDMaps) of presence-only python observations (EDDMaps, 2014) to estimate these parameters. The intrinsic rate of growth was estimated in two ways: (a) using the exponential rate of growth in python observations in EDDMaps; and (b) using published information on survival and reproductive rates (Willson et al., 2011). First, we restricted EDDMaps observation records to the period 1995–2010: only one observation occurred prior to this in 1979 and population growth may have been limited by an unusually cold weather in January 2010. We also limited observations to those occurring south of the northern border of the refuge to eliminate the possibility of independent introductions of pythons in Florida distant from the original introduction(s) in the Everglades. We used the model: log(ct ) = 0 + (t),

(5)

where t is the year and c is the number of observations. For this model  = 0.429 (SD=0.030) (F114 = 210.7, P < 00001). Second, using the estimates or assumed values of age-specific reproduction and survival provided by Willson et al. (2011), we constructed the following projection model for female pythons aged 1, 2, and older:

⎛

Y1,t+1

⎜ ⎝ Y2,t+1 A3+,t+1

⎞

⎞ ⎛ Y1,t ⎞ ⎟ ⎝ ⎟ ⎜ ⎠ = 0.5 0.0 0.0 ⎠ ⎝ Y2,t ⎠ ⎛

0.0

0.0

3.5

0.0

0.7

0.9

(6)

A3+,t

The dominant eigenvalue of the projection matrix provides the asymptotic growth rate  = 0.468, which compares favorably with that estimated from the growth in python observations. We used an estimate of  = 0.468 for the purposes of simulating our model. The diffusion coefficient was estimated using the method described by Holmes (1993), in which the square root of the area occupied is plotted over time. According to Skellam (1951): √ √ at = 2t D, (7) where a is the occupied area. To estimate at we determined the area of the convex hull encompassing all observations from EDDMaps south of the refuge during 1995–2010. The estimated slope of this model was 2(D)1/2 = 9.753 (SD=0.454) (F115 = 461.4, P < 0.001). Using  = 0.468 we derived an estimate of D = 16.2 km2 /yr. We chose this method over that based on radio telemetry because of the inherent difficulties in differentiating exploratory and dispersal movements, in specifying an appropriate time step

Table 1 Parameters of the LCM model used for pythons.

 = 0.468 D = 16.2 (km2 /yr) K = 100 pythons. dm = 39.4 (km/yr) ˛T1 = 0.01 ˛T2 = 0.01, 0.1, 0.3 and 0.5 t0 = 1995 xFlamingo =(508790 E,2821227 N) (UTM coordinates, zone 17, northern hemisphere) 0 nt (x) = 0 if x = / xFlamingo 0 nt (xFlamingo ) = K

between observations, and in ensuring an adequate sample of long-distance dispersal (Turchin, 1998). The method that we used, described by Holmes (1993), assumes exponential growth of the population and a random walk of individuals. To help ensure exponential growth, we limited data to the period of 1995-2010, during which the population appeared to be growing exponentially. We also tested whether six radio-tagged pythons released at their capture site (K.M. Hart, personal communication) appeared to follow a random walk using the method described by Swihart and Slade (1985). Based on successive locations of pythons (n = 184) during 2006–2008, we could not reject the null hypothesis of a random walk or overdispersion (as opposed to aggregation) for any of the six pythons (mean (t2 /r2 ) = 5.10, experimentwise error rate P > 0.26). The same data were used to provide an estimate of the average yearly movement rate, dm = 39.4. Information on carrying capacity, population size, date, and location of the python introduction is much more speculative. We could find no empirical information on carrying capacity of pythons, although we assume it is lower than the 200 brown tree snakes per km2 provided by Kaiser and Burnett (2010). We admit that the carrying capacity should ideally vary as a function of the habitat (i.e. mangrove, marsh, urban and agricultural areas), but without the necessary information we used an average estimate. For illustrative purposes we set K = 100, although this is merely a guess on our part. The best information about the date, location, and size of the python introduction in southern Florida is provided by Willson et al. (2011), although there is considerable uncertainty in these values. For our purposes, we initialized the model in 1995, with K individuals located at the center of the spatial distribution of the observations from EDDMaps (near Flamingo, Florida: xFlamingo =UTM 508790E, 2821227N). Thus n0 (xFlamingo ) = K and n0 (x) = 0, ∀ x ∈  \ xFlamingo . The model parameters are summarized in Table 1. Model validation. We tried to assess the validity of the model with the presence-only data from EDDMaps in 2013. We used our model to predict the distribution of pythons in 2013 without any control action. Then, we compared the predicted infested zone from our model with that suggested by the presence-only data. For the RD model, we considered a cell to be infested if the population density of pythons was greater than 0.1. We recall that because we considered no control actions were implemented, predictions made by the CEH and LCM models are equal. 3. Results 3.1. Comparison with the CEH model on a hypothetical example As it is shown in Fig. 3, for every simulated example, the CEH model produced a much higher prediction error than the LCM model. The percentage of error with the CEH model increases with the percentage of removed individuals. This is simply due to the fact that the highest control mortalities are associated with the highest effect of control on dispersal, which is not accounted for in the CEH

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

Fig. 3. Percentage of increased error produced by the CEH model as a function of the percentage of removed individuals at the control locations, nR . There is one line for the low dispersal scenario and for the high dispersal scenario. We recall that the error of the CEH model is equal to Error(LCM)(1 + pctLCM /100).

model. The error of the CEH model is higher in the low dispersal scenario. On every simulated example, the error of the CEH model is at least 40% higher than with the LCM model (i.e. pctLCM > 40 %). As illustrated by the examples provided in Fig. 4, most of the prediction error occurs around the squared controlled area, but not inside, which is explained by the fact that either a larger (or a fewer number depending on the situation), of individuals escape the controlled area. Finally, the prediction error is the highest on the straight line controlled area, especially for the CEH model. As opposed to the square control area, there is not enough controlled cells for the CEH model to reduce the population as it is in the simulation. 3.2. Management of Burmese pythons No control action. For comparison, we first used the RD model to predict the density of pythons in the refuge when no control action is implemented (Fig. 5). We first note that their number increases exponentially. Second, the model predicts that pythons

21

Fig. 5. Predicted number of pythons in the refuge when public areas are controlled with a control mortality of ˛T1 = 0.01 and increasing control mortality on roads and around the refuge ˛T2 .

will be present in the refuge in 2015 (approximately 4 pythons) and thus, pythons are already in the refuge when we start the simulation of the control method. Finally, at the end of the simulation, the refuge is still not saturated. Increasing control mortality on roads. Fig. 5 shows the predicted change in the number of pythons in the refuge when public areas are controlled with a control mortality of ˛T1 = 0.01 and increasing control mortality on roads and levees around the refuge ˛T2 = 0.01, 0.1, 0.3, 0.5. We first note that none of the control mortalities are sufficient to eradicate pythons from the refuge. As is shown in Fig. 6, the proportion of all cells with increased control mortality is low. Thus it is still possible for pythons to reach the refuge using some uncontrolled path (e.g. from the east) or to travel via the public areas where the control mortality is very low. In addition, once pythons cross the border of the refuge, the control mortality is too low to control the population inside the refuge. We also provide in Fig. 7, the percentage decrease in the number of pythons compared to the case where no control action is applied. For each time horizon, the percentage decrease is linear in the control mortality ˛T2 and the percentage of decrease is higher for the

Fig. 4. Examples of the per cell average absolute difference between the simulated and the predicted maps for each model. The average absolute error is computed using Eq. (4) for each cell, without the second summation. The error values of the CEH and LCM model are provided on figure (a) and (b) respectively. Both examples are for the low-dispersion scenario and a removal rate of nR = 60%. Yellow cells indicate a maximal error of 60. In other words, the difference between the predicted model and the simulation in this cell is, on average, equal to 60. Dark blue cells indicate that prediction and simulation are equal. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)

22

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

Fig. 6. (a) In black are the positions of the cells near roads and around the refuge, where the control mortality ˛T2 is used. ˛T1 is used elsewhere inside the red zones. (b) Predicted density of pythons in the domain in 2030 and with ˛T2 =0.5 . K = 100 is the carrying capacity, i.e. the maximal number of pythons that can live in a 1 km2 cell.

Fig. 7. Decrease in the number of pythons in the refuge when public areas are controlled with a control mortality of ˛T1 = 0.01 and increasing control mortalities ˛T2 on roads and around the refuge. The decrease of pythons is expressed as a percentage compared to the case where no control action is implemented. There is one line per each time horizon.

two last times periods(i.e. 2025 and 2030). With the lowest control mortality, the population of pythons in the refuge is only decreased by 14% in 2030, where it is decreased by 43% for the highest control mortality. Quarantine area over the refuge. We provide the results of the analysis as a regression tree (Fig. 2b.). For the different values of the time horizon T, control mortality ˛, and quarantine width h, the results of the simulation can either be S for success, i.e. the density of python in the Arthur R. Marshall Loxahatchee National Wildlife Refuge is lower than 1 individual per km2 or F for failure, if the population in the refuge is greater than 1 ind/km2 . As is shown earlier, the population in the refuge remains low in 2020 and any control mortality and quarantine width are sufficient to maintain the density under 1 ind/km2 . For the other time horizon, i.e. t=2025, 2030, the population in the refuge is under 1 ind/km2 only when the control mortality is higher than 0.2 and the quarantine width is higher than 15 km. Model validation. While validation of our model for pythons is not possible, its realism may be judged by comparing the predicted range in 2013 with that indicated by the presence-only records in EDDMaps (Fig. 8). The greatest discrepancy is in the northwest section of the domain, suggesting that a model incorporating

Fig. 8. Spatial range of the infestation determined by the RD model and as represented by the EDDMaps data. The black dots are the presence-only observations from the data base from 1995 to 2013. In grey is the predicted occupied area using the RD model.

directional spread (advection), a correlated random walk, or heterogeneous environment might be more appropriate (Holmes et al., 1994). A secondary point of introduction could also explain this discrepancy. 4. Discussion In this article we proposed a simulation based framework, which allows predictions of the effect of different management options (i.e. varying controlled locations and control mortality) on the spatial distribution and density of an invasive species. This framework is an extension of the CEH model and allows consideration of situations where control is applied both in-situ and during dispersal. The extension of the CEH model is justified by mathematical arguments. Using a simulated dataset, we showed that the LCM model actually provides better estimates when control is also applied during dispersal. RD models have been extensively used in the literature to characterize the invasion of an invasive species. However, only the

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

theoretical characteristics of the model have been used, such as the theoretical speed of the invasion front (Carrasco et al., 2010) or the critical length of suitable habitat (Neubert, 2003). To the best of our knowledge, only (Kaiser and Burnett, 2010) used a RD model to predict the spatial distribution of an invasive species under control. They explored early detection and rapid response for Brown treesnakes (Boiga irregularis) on the Oahu Island, Hawaii. Unfortunately, the effect of control on the invasive population lacks any mathematical or biological justification. Indeed, they considered that control in one cell affects the population in all other cells, no matter the distance between cells. For example, if only one cell is controlled, the effect on adjacent cells is the same as the effect on cells located possibly hundreds of kilometers away. In our approach, the effect of control during dispersal is explicitly added in the proof of the Fisher’s equation and the LCM model is obtained by mathematical argument. First, we proposed to demonstrate that the LCM model provides a better prediction than the CEH model when control is also applied during dispersal by using a simulated dataset. We simulated the dataset using a cellular automaton model, which is particularly suitable as it assumes that individuals follow a random walk, the basic assumption of the RD model. We compared the prediction of the CEH and the LCM model to the one provided by the cellular automaton model and found that for various diffusion and growth rates, and control scales, the predictions of the LCM model were closer to the simulated one. Second, we used the simulation framework to discuss the management of Burmese pythons in the Everglades and its implications for the Arthur R. Marshall Loxahatchee National Wildlife Refuge. We used a realistic situation, where control is available only on public lands located south of the refuge and where the control mortality can be increased only on roads and on the border of the refuge. Our most efficient management scenario is based on a control mortality of 0.5 on roads and around the refuge, which is far from the likely present value of about 0.01. Using our simulation model, none of the management scenarios were sufficient to control the population within the refuge. The maps of predicted python density seem to indicate that the proportion of control cells with increasing control mortality remains too low to significantly impact the python population. Thus, the highest priority in the effort to control pythons in southern Florida appears to be finding methods that can increase the detection rate, and thus capture mortality, of this cryptic predator (Dorcas and Willson, 2013). We proposed a second scenario to understand how big the controlled area should be to maintain a low density of pythons in the refuge. We found that a control mortality greater than 0.2 and a quarantine area of 15 km were needed to achieve this objective. Hopefully, this analysis shows that the invasion of pythons in the refuge can be avoided when sufficient area is controlled. But once again, this remains dependent on greatly increasing the control mortality, which will be difficult in practice. Indeed, a control mortality ≤0.2 is predicted to be insufficient, even for the wider quarantine area. RD models make the implicit hypothesis that dispersion of the individuals is random, but it is not clear to what extent this hypothesis can be violated and the RD model still be an adequate approximation. Indeed, even if at the individual scale the diffusion is random, at the population scale the density of individuals is smooth and the invasion front is regular. RD models is a way to represent a regularly growing and diffusing population and appear to be sufficiently general to fit many situations. Nonetheless, RD models propose to model the short-diffusion of individuals as a regular and deterministic spatial process. But in practice, some stochastic events can modify this regular dispersion. We show in an appendix that one can easily extend RD models in order to account for long-distance dispersal events (LDD) and secondary points of

23

introduction (SPI). The resulting RD model would be stochastic. Many other extensions of the RD models are available in the literature. For example, one can use varying growth rates and diffusion coefficients to mimic heterogeneous landscapes. Another extension would add an Allee effect. Another example is the Lotka–Volterra predator–prey model (Murray, 1975). In this case, the RD model is now constituted of two correlated partial differential equations, one for the predator population dynamics and the other for the prey. When the predator is an invasive species and the prey an endangered species, one can predict the effect on the endangered species population of controlling its non-native predator. Or one can also compute the minimal control effort on the invasive species that guarantees a fixed minimal abundance of the endangered species. Other extensions are possible and in general are well documented. See for example (Cantrell and Cosner, 2003) for a detailed presentation of RD models and their possible extensions. RD models can also be used in an iterative fashion to consider non-fixed management strategies (commonly called dynamic or adaptive strategies), where the control locations and/or the control mortality is not constant over time. Stochastic RD models can be considered with the simulation of LDD and SPI events and then can be useful to define the transition probabilities of any Markov decision process. Indeed, the transition probabilities can be obtained by using Monte-Carlo simulations, based on our model. Finally, the approximation of the probability that an immigrant crossed a controlled location, rTj , can certainly be improved by more complex computation. Unfortunately, we don’t know at this time the strength of the approximation that has been proposed in this article. For the purpose of simulating the effect of a management strategy, other models can be found in the literature (Wikle and Hooten, 2010). In particular, one can prefer models that are both discrete in time and space. But their implementation is generally more difficult and their relative improvement compared to RD models is not clear. We believe that RD models make a good trade-off between model complexity and flexibility for use in different situations. In addition, due to the generality of RD models and the relative ease with which they can be parameterized, the framework proposed in this article might be suitable for the study of different management problems and invasive species. And more generally, this framework can be used in other contexts; for example in the control of the spread of disease. In addition, the high level of detail we provide on the implementation of the RD model should facilitate its use for other applications. Data accessibility Pythons observations: Pythons observations are available online (EDDMaps, accessed June 2014). Road locations: Derived from the TRIGER/Line shapefiles (Triger, accessed 21 August 2014). Public area locations: Uploaded as online supporting information. Matlab scripts: Uploaded as online supporting information. Acknowledgements We would like to thank Miguel Acevedo and Mariano Marcano for fruitful discussions about the use of the Crank–Nicholson method, and Sergei Pilyugin for help with reaction–diffusion models in general. We also thank Ann Foster for providing the necessary GIS data and Kristen Hart for providing telemetry data for the python example. This study was funded by the U.S. Geological Survey’s (USGS) Invasive Species Program and Greater Everglades Priority Ecosystem Studies. We thank the U.S. Fish and Wildlife Service for hosting a workshop concerning python control at the

24

M. Bonneau et al. / Ecological Modelling 337 (2016) 15–24

National Conservation Training Center. Any use of trade, product, or firm names in this article is for descriptive purposes only and does not imply endorsement by the U.S. Government. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ecolmodel.2016. 05.013. References Acevedo, M., Marcano, M., Fletcher Jr., R., 2012. A diffusive logistic growth model to describe forest recovery. Ecol. Model. 244, 13–19. Andow, D., Kareiva, P., Levin, S., Okubo, A., 1990. Spread of invading organisms. Landsc. Ecol. 4 (2), 177–188. Bogich, T., Liebhold, A., Shea, K., 2008. To sample or eradicate? A cost minimization model for monitoring and managing an invasive species. J. Appl. Ecol. 45 (4), 1134–1142. Cantrell, R., Cosner, C., 1996. Spatially explicit models for the population dynamics of a species colonizing an island. Math. Biosci. 136, 65–107. Cantrell, R., Cosner, C., 2003. Spatial Ecology via Reaction–Diffusion Equations. John Wiley & Sons, Ltd. Carrasco, L., Baker, R., MacLeod, A., Knight, J., Mumford, J., 2010. Optimal and robust control of invasive alien species spreading in homogeneous landscapes. J. R. Soc. Interface 7, 529–540. Clark, C., 1990. Mathematical Bioeconomics: The Optimal Management of Renewable Resources, Pure and Applied Mathematics, a Wiley-Interscience Series of Texts, Monographs and Tracts. Dehnen-Schmutz, K., Touza, K., Perrings, C., Williamson, M., 2007. The horticultural trade and ornamental plant invasions in Britain. Conserv. Biol. 21, 224–231. Dorcas, M., Willson, J., 2013. Hidden giants: problems associated with studying secretive invasive pythons. In: Lutterschmidt, W.I. (Ed.), Reptiles in Research. Nova Biomedical, New-York, pp. 367–385. Dorcas, M., Willson, J., Reed, R., Snow, R., Rochford, M., Miller, M., Meshaka, W., Andreadis, P., Mazzotti, F., Romagosa, C., H.K.M., 2012. Severe mammal declines coincide with proliferation of invasive Burmese pythons in everglades national park. Proc. Natl. Acad. Sci. U. S. A. 109 (7), 2418–2422. EDDMaps, 2014. Early Detection & Distribution Mapping System. The University of Georgia – Center for Invasive Species and Ecosystem Health, Available online at: http://www.eddmaps.org/florida/species/subject.cfm?sub=20461 (accessed 05.05.14). Epanchin-Niell, R., Hastings, A., 2010. Controlling established invaders: integrating economics and spread dynamics to determine optimal management. Ecol. Lett. 13 (4), 528–541. Fisher, R., 1937. The wave of advance of advantageous genes. Ann. Eugen. 7, 355–369. Frid, L., Hanna, D., Korb, N., Bauer, B., Bryan, K., Martin, B., Holzer, B., 2013. Evaluating alternative weed management strategies for three Montana landscapes. Invasive Plant Sci. Manage. 6 (1), 48–59. Giljohann, K., Hauser, C., Williams, N., Moore, J., 2011. Optimizing invasive species control across space: willow invasion management in the Australian Alps. J. Appl. Ecol. 48, 1286–1294. Grevstad, F., 2005. Simulating control strategies for a spatially structured weed invasion: Spartina alterniora (loisel) in pacific coast estuaries. Biol. Invasions 7 (4), 665–677. Hastings, A., Cuddington, K., Davies, K., Dugaw, C., Elmendorf, S., Freestone, A., Harrison, S., Holland, M., Lambrinos, J., Malvadkar, U., Melbourne, B., Moore, K., Taylor, C., Thomson, D., 2005. The spatial spread of invasions: new developments in theory and evidence. Ecol. Lett. 8, 91–101. Hauser, C., McCarthy, M., 2009. Streamling ‘search and destroy’: cost-effective surveillance for invasive species management. Ecol. Lett. 12, 683–692. Higgins, S., Richardson, D., Crowling, R., 2000. A dynamic spatial model for man aging alien plant invasions at the landscape extent. Ecol. Appl. 10, 1833–1848.

Holmes, E., Lewis, M., Banks, J., Veit, R., 1994. Partial differential equations in ecology: spatial interactions and population dynamics. Ecology 75, 17–29. Holmes, E., 1993. Are diffusion model too simple? A comparison with telegraph models of invasion. Am. Nat. 142, 779–795. Kaiser, B., Burnett, K., 2010. Spatial economic analysis of early detection and rapid response strategies for an invasive species. Resour. Energy Econ. 32, 566–585. Keller, R., Springborn, M., 2013. Closing the screen door to new invasions. Conserv. Lett. 00, 1–8. Keller, R., Lodge, D., Finnoff, D., 2007. Risk assessment for invasive species produces net bioeconomic benefits. Proc. Natl. Acad. Sci. U. S. A. 104, 203–207. Kurata, K., Shi, J., 2008. Optimal spatial harvesting strategy and symmetry breaking. Appl. Math. Optim. 58, 89–110. Leroux, S.J., Larrivée, M., Boucher-Lalonde, V., Hurford, A., Zuloaga, J., Kerr, J.T., Lutscher, F., 2013. Mechanistic models for the spatial spread of species under climate change. Ecol. Appl. 23 (4), 815–828. Lodge, D., Williams, S., MacIsaac, H., Leung, B., Reichard, S., Mack, R., Moyle, P., Smith, M., Andow, D., Carlton, J., McMichaell, A., 2006. Biological invasions: recommendations for the U.S. policy and management. Ecol. Appl. 16, 2035–2054. Mehta, S., Haight, R., Homans, F., Polasky, S., Venette, R., 2007. Optimal detection and control strategies for invasive species management. Ecol. Econ. 61, 237–245. Meier, E., Dullinger, S., Zimmermann, N., Baumgartner, D., Gattringer, A., Hülber, K., 2014. Space matters when defining effective management for invasive plants. Biodivers. Res. 20, 1029–1043. Murray, J., 1975. Non-existence of wave solution for the class of reaction–diffusion equations given by the Volterra interacting-population equations with diffusion. J. Theor. Biol. 52, 459–469. Neubert, M., 2003. Marine reserves and optimal harvesting. Ecol. Lett. 6, 843–849. Nicol, S., Chadès, I., 2011. Beyond stochastic dynamic programming: a heuristic sampling method for optimizing conservation decisions in very large state spaces. Methods Ecol. Evol. 2 (2), 221–228. Olson, L., 2006. The economics of terrestrial invasive species: a review of the literature. Agric. Resour. Econ. Rev. 35, 178–194. Oruganti, S., Shi, J., Shivaji, R., 2002. Diffusive logistic equation with constant yield harvesting, I: Steady states. Trans. Am. Math. Soc. 354, 3601–3619. Provencher, L., Forbis, T., Frid, L., Medlyn, G., 2007. Comparing alternative management strategies of fire, grazing, and weed control using spatial modeling. Ecol. Model. 209, 249–263. Pysek, P., Richardson, D., 2010. Invasive species, environmental change and management, and health. Annu. Rev. Environ. Resour. 35, 25–55. Schaefer, M., 1957. Some considerations of population dynamics and economics in relation to the management of the commercial marine fisheries. J. Fish. Board Canada 14, 669–681. Schapaugh, A., Tyre, A., 2012. A simple method for dealing with large state spaces. Methods Ecol. Evol. 3 (6), 949–957. Shigesada, N., Kawasaki, K., 1997. Biological invasions: theory and practice. Oxf. Ser. Ecol. Evol. Skellam, J., 1951. Random dispersal in theoretical populations. Biometrika 38, 196–218. Strikwerda, J., 2004. Finite Difference Schemes and Partial Differential Equations. Siam. Svenning, J., Gravel, D., Holt, R., Schurr, F., Thuiller, W., Münkemller, T., Schiffers, K., Dullinger, S., Edwards, T., Hickler, T., Higgins, S., Nabel, J., Pagel, J., Normand, S., 2014. The influence of interspecific interactions on species range expansion rates. Ecography, 1–12. Swihart, R., Slade, N., 1985. Testing for independence of observations in animal movements. Ecology 66, 1176–1184. Turchin, P., 1998. Quantitative Analysis of Movement: Measuring and Modeling Population Redistribution in Animals and Plants. Sinauer Associates, Inc., Sunderland, MA. Wikle, C.K., Hooten, M.B., 2010. A general science-based framework for dynamical spatio-temporal models. Test 19 (3), 417–451. Willson, J., Dorcas, M., Snow, R., 2011. Identifying plausible scenarios for the establishment of invasive Burmese pythons (Python molurus) in Southern Florida. Biol. Invasions 13, 1493–1504.