Speeding up the simulation of population spread models

Simulating spatially explicit population models to predict population spread allows environmental managers to make better‐informed decisions. Accurate simulation requires high spatial resolution, which, using existing techniques, can require prohibitively large amounts of computational resources (RAM, CPU, etc). We developed and implemented a novel algorithm for the simulation of integro‐difference equations (IDEs) modelling population spread, including stage structure, which uses adaptive mesh refinement. We measured the accuracy of the adaptive algorithm by comparing the results of simulations using the adaptive and a standard non‐adaptive algorithm. The relative error of the population's spatial extent was low (<0·05) for a range of parameter values. Comparing efficiency, we found that our algorithm used up to 10 times less CPU time and RAM than the non‐adaptive algorithm. Our approach provides large improvements in efficiency without significant loss of accuracy, so it enables faster simulation of IDEs and simulation at scales and at resolutions that have not been previously feasible. As an example, we simulate the spread of a hypothetical species over the UK at a resolution of 25 m. We provide our implementation of the algorithm as a user‐friendly executable application.


Introduction
Modelling changes in species' distributions enable environmental managers to make better-informed habitat management (Hulme 2006) and conservation (Guisan, Tingley & Baumgartner 2013) decisions. Predicting spread has two major applications: to the problem of invasive species (Shigesada & Kawasaki 1997;Williamson 1999) and to understanding population responses to the shifting of habitat under climate change (Zhou & Kot 2011;Bullock et al. 2012;Bennie et al. 2013), habitat loss (Fahrig 1997) and habitat fragmentation (Fahrig 2002). This provides motivation for spatially explicit population models which represent the relevant biological processes mathematically and enable the prediction of population dynamics. As population spread happens over large (Continental/National) spatial scales, large-scale models are required. One way of doing this is via species distribution models (e.g. Elith & Leathwick 2009), but these are static and predict only potential ranges which may not be realised (Zhu, Woodall & Clark 2012). They do not generally include the biological mechanisms of demography or dispersal which are likely to influence the rate and extent of spread (Schurr, Pagel & Cabral 2012) and have a comparatively weak underpinning in ecological theory (Thuiller et al. 2013). SDMs are trained only on data from locations within the range of the species. Predictions outside of these ranges may be unrealistic as they are beyond the scope of the data, making predictions at range limits or scenario testing (e.g. adding another species or biological control) unreasonable. These problems can be overcome by using mechanistic models that explicitly include demography and dispersal. Mechanistic models take several forms, including partial differential equations (PDEs), integro-difference equations (IDEs), individual-based models (IBMs) and some metapopulation models, with the choice of model depending on the temporal and spatial scales underpinning the system and the form of the biological processes (e.g. whether dispersal can be characterised by diffusion or is more complex). Another option is hybrid models that use components of SDMs and mechanistic models, but these are vulnerable to over-fitting, caused by needing to estimate or assume large numbers of parameter values (Wisz, Pottier & Kissling 2013).
In many mechanistic population models, landscapes are assumed to be homogeneous, with the demographic and dispersal parameters spatially constant. This is because homogeneity (i) reduces the number of independent parameters in the model, making model specification easier and (ii) reduces the amount of computational resources required to determine the spread rate. In homogeneous landscapes, population growth can be characterised by a single scalar growth rate. Similarly, in spreading populations, one can calculate the spreading speed (the annual increment in the population's range) and how sensitive this is to life-history and dispersal parameters (Shigesada, Kawasaki & Teramoto 1986;Kot, Lewis & den Driessche 1996;Neubert & Caswell 2000). For landscapes where spatial variation is very slight (Gilbert et al. 2014a) or highly localised (Dewhirst & Lutscher 2009), the results for homogeneous landscapes can be used as an approximation, with analogous results extending to some periodic (e.g. regular, repeating) landscapes (Gilbert et al. 2014b). However, landscape homogeneity is often a poor assumption and in many landscapes the dispersal and demographic parameters exhibit strong variation across large spatial scales, rendering this approximation inaccurate (Miller & Tenhumberg 2010;Svenning, Gravel & Holt 2014). Instead, to make predictions on complex landscapes, numerical simulations of the models are required. For these predictions to be accurate, both accurate ecological data and accurate simulation of the model are required. Accurate simulations require a large number of calculations (Scheffer et al. 1995). This means that intense computational power is needed to simulate spread in large two-dimensional landscapes. This is a particular problem for IBMs, such as RangeShifter (Bocedi et al. 2014), where a large number of individual mechanisms can be modelled, but where each organism must be explicitly represented in the numerical framework.
As predictions of species spread rely on computationally expensive simulations, the ability to predict spread is often limited by the model's scale and the available computational resources. This limitation inhibits the types of scientific questions and problems that can be investigated. For example, fitting spread models to data requires numerous model runs, so if each model run is computationally intensive, then the time to fit the model to data will be prohibitively long.
To be simulated numerically, most mechanistic population models must be discretised, with space and time represented as a mesh of points. The simulation's accuracy depends on the number of points, with the accuracy increasing with the number of points (Bocedi et al. 2012). Simulating spatially explicit models involves a large number of calculations at each time step. The computational complexity and computational resources (RAM, CPU time, etc) required to compute these calculations increase linearly or faster with the total number of points or individuals (Cooley & Tukey 1965;Solodovnikov 1985). If simulations over large landscapes are to be accurate, they need a high density of points and a very large total number of points. In addition, most models will require numerous iterations if the model is simulated over a long time, making simulations extremely computationally expensive. Large numbers of runs are also often required for parameter estimation and sensitivity studies, thus further motivating the need for efficient solvers (Csill ery et al. 2010). One way to improve efficiency has been developed for PDEs, using adaptive mesh refinement to allow the density of points (the resolution) to change in time and space (e.g. Davis & Flaherty 1982;Berger & Oliger 1984). The density of points is usually highest where the most accuracy is required (e.g. where the relative population change is greatest, such as in the wave-front). With only small reductions in accuracy (i.e. the accuracy of the simulations compared to the mathematical model), this approach allows limited computational resources to be focused on the spatial locations which have the greatest effect on the distribution as it changes. Adaptive mesh refinement has so far been restricted to PDEs. However, in cases where long distance dispersal or annual cycles play an important role in the model, IDEs may be a more appropriate choice of modelling framework than PDEs, as they can incorporate these processes (Kot, Lewis & den Driessche 1996) and are relatively easy to parameterise, e.g. by using data from the COMPADRE database (Salguero-G omez, Jones & Archer 2015) or population spread data (e.g. Gilbert et al. 2005), and thus, IDEs are extensively used in species modelling (Neubert & Parker 2004;Le Corff & Horvitz 2005;Bullock, Pywell & Coulson-Phillips 2008;Miller & Tenhumberg 2010;Zhou & Kot 2011;Bullock et al. 2012).
In this study we provide a solution to large-scale spread problems by developing an adaptive mesh refinement algorithm to simulate IDEs. At each iteration, the algorithm divides the landscape up into square regions of equal area that are either high resolution (containing a large number of mesh points) or low resolution (containing only one mesh point). We present the algorithm for stage-structured populations with homogeneous dispersal and heterogeneous demography in Section 'Materials and methods'. In Section 'Results', we demonstrate that the algorithm is both highly accurate and uses fewer computational resources (RAM, CPU time) than non-adaptive algorithms for a given problem, as well as enabling us to tackle problems that were previously impossible due to insufficient computational resources. Finally, we discuss extending the algorithm to landscapes with heterogeneous dispersal in the Discussion (Section 'Discussion'). We provide an executable application with a worked example via GitHub (Gilbert et al. 2016) so that others may use this algorithm to simulate the spread of species under a variety of scenarios.

Materials and methods
In this section, we introduce single species stage-structured IDEs with constant dispersal and show how to simulate them with adaptive meshes and achieve improvements in computational efficiency compared to non-adaptive meshes. We restrict ourselves to IDEs on twodimensional landscapes, with the method easily restricted to one spatial dimension. The technical details of the algorithm are presented in full in Appendix S2.
Integro-difference equations are discrete time, spatially continuous population models that have been used to model many different organisms including annual (Bullock, Pywell & Coulson-Phillips 2008) and perennial herbs (Le Corff & Horvitz 2005), shrubs (Neubert & Parker 2004), trees , weevils (Miller & Tenhumberg 2010) and butterflies (Zhou & Kot 2011). Their strength lies in the trade-off between simplicity, flexibility and ease of use in terms of parametrisation and implementation; they can incorporate a range of dispersal mechanisms and demographic structures, including stage-structured populations which represent each life-history stage explicitly (Neubert & Caswell 2000), but with few parameters. The general stage-structured two-dimensional IDE on a landscape Ω relates the stage-structured vector population density u tþ1 ðxÞ at each location x for generation t + 1 with the population for generation t via where denotes the Hadamard (element-wise) product of two matrices.
Here, the population undergoes two consecutive phases, growth and dispersal. The growth phase incorporates the processes of birth, maturation and death that occur during a generation. Growth is local and depends only on the location x and on the population u t ðxÞ at x at time t. Locations may differ in terms of habitat, meteorological conditions, soil quality, etc., impacting growth rates and carrying capacities. The current population may affect demographic rates through intraspecific competition or via Allee effects. Population growth is determined by the population projection matrix Bðu t ðyÞ; yÞ, where the (i,j)-th element is the ratio of the number of individuals in stage j after the growth phase and the number of individuals in stage i at time t (at location y), e.g. the survival probability between juveniles and adults. The dispersal phase is the process by which individuals (or propagules) move between locations. In stage-structured IDEs, it is characterised by the matrix of dispersal kernels K(x À y,y), where the (i,j)-th element K i;j ðx À y; yÞ gives the density of population that has transitioned from stage i to j and has transferred between locations y and x, e.g. dispersal of seeds from adult plants. Each pair of demographic stages is considered separately as individuals' dispersal depends on both their current stage and their stage prior to the growth phase.
In general, the dispersal kernel depends on the disperser's starting location y and the displacement x À y (as in Fig. 1e). However, in this study, we consider only models where the dispersal kernel is spatially constant and has no explicit y dependence (see Fig. 1d), which is a useful simplifying assumption. Prospects for relaxing this assumption are reviewed in the Discussion.
With the exception of homogeneous landscapes (see section 'Analytical Results' in Appendix S1, Supporting Information), quantifying spread requires numerical simulation, with the landscape typically discretised into a square mesh (Powell, White & McMillen 1998). If dispersal is homogeneous, then the integral in (1) can be computed efficiently using fast Fourier transforms (FFTs) (Andersen 1991). If dispersal varies spatially, then the discretised spatial distributions at t and t + 1 are related by a matrix multiplication. Considering either case, the computational resources required increase linearly or faster with the number of mesh points, leading to very large computational costs for simulating the model accurately over large and complex landscapes (see Appendix S1). To reduce these computational requirements, we propose a simulation algorithm with adaptive mesh refinement.
In spreading populations, there is usually little change in the population density in the wave-back (or within range), where the population has reached carrying capacity, or in the far-field, where the population is very small (see Fig. 1a). Behaviour in these regions has only a small effect on the wave-front's behaviour, so computational resources can be focused on the wave-front (or leading edge) with a loss of resolution in the wave-back and far-field. For our novel algorithm, we subdivide the landscape into non-overlapping square regions of equal width Dx, which are further divided into an equal number of cells of width dx. Regions can be either High Resolution (HR) or Low Resolution (LR have a mesh point for each cell, but low-resolution regions are represented by a single mesh point. The algorithm requires the following inputs: maps (provided as raster, i.e. gridded data at the desired resolution) of the initial population distribution, demographic parameters and corresponding growth functions over the entire landscape and the (spatially homogeneous) dispersal kernel. Initially, any region in which any mesh point is in the wavefront is high resolution (see Fig. 1a). The wave-front is any location where the juvenile (or lowest) demographic stage, u t 1 ðxÞ, has population density between the thresholds T LR and T HR . The upper threshold T LR is taken as a proportion of the carrying capacity and T HR is set as small as required to allow small populations to be resolved (e.g. T HR ¼ 10 À13 in the examples in Sections 'Materials and methods' and 'Results'), with suitable values provided for these parameters in the worked example.
In standard simulations of IDEs there are two sequential processes during each generation: population growth and dispersal. However, the adaptive algorithm also incorporates the additional process of mesh refinement. 1. Population Growth at each cell depends on the cell's intrinsic demographic rates and carrying capacity, which are supplied by the user. For high-resolution regions, population growth is calculated by applying each cell's population projection matrix B(, y) to its stage-structured population density u t ðyÞ. For low-resolution regions, each demographic stage's population density is represented by a single scalar, so the growth phase for the whole region is calculated using the mean of the population projection matrix over the region (see Appendix S2, equation (B2)).

Dispersal.
The assumption of homogeneous dispersal reduces the integral for each pair of demographic stages in (1) to a convolution and allows dispersal within a regular square mesh to be calculated using FFTs (Andersen 1991). The key idea with our method is to divide dispersal into two categories of dispersal resolution.
(a) High-resolution dispersal is used to denote the dispersal processes which require more computational resources for simulations to be accurate. High resolution is required for dispersal that (i) is relatively short distance and (ii) has its destination in the wave-front. In practice, this is dispersal to a high-resolution region (a region in the wave-front) from any neighbouring region (short distance). For a region containing high-resolution information, we take the pre-dispersal population density of each cell in the region and its neighbours (for low-resolution neighbouring regions, all cells are assumed to have population density equal to the region's population density) and take the discrete convolution of this distribution with the HR dispersal kernel (a high-resolution discretisation of the dispersal kernel) to get the post-dispersal population density distribution. (b) Low-resolution dispersal is used to denote the dispersal processes which can be simulated with fewer computational resources without significantly affecting the simulation's accuracy. This is dispersal between regions which are not both neighbours and in the wave-front. This is justified by the observation that all realistic dispersal kernels eventually decay after sufficient distance from the origin. For each pair of demographic stages, the dispersal between non-neighbouring regions is calculated by taking the convolution of the array of mean pre-dispersal population densities for each region with the long-distance LR dispersal kernel, a low-resolution discretisation of the dispersal kernel restricted to dispersal between non-neighbouring regions. For dispersal to low-resolution regions from neighbouring regions, we take the mean pre-dispersal population densities for the region and its neighbours and take the discrete convolution of this distribution with the local LR dispersal kernel, a low-resolution discretisation of the dispersal kernel restricted to dispersal between neighbouring regions. Finally, the contributions of low-resolution dispersal are added to the post-high-resolution dispersal population distribution. For low-resolution regions, this involves the addition of two scalars (for each pair of demographic stages). For high-resolution regions, the contribution of long distance dispersal to a region is approximated as equal for all cells, so the region's global dispersal density is added to each cell's population density.
In summary, picking and choosing which level of dispersal resolution to use, based on the change in population and distance between locations, allows gains in efficiency with little reduction in accuracy (See Appendix S2 for full details). 1. Mesh refinement is where the resolution of each region is changed to ensure that while the population distribution and the wave-front (the area with population density between the lower T HR and upper T LR thresholds) change, the mesh always has high resolution where needed. Individuals in the far field may drive population expansion (Lewis 2016), so the choice of T HR will depend on the IDE model being studied and is in general set as low as required to ensure no significant change in model prediction. The mesh refinement phase enables regions which enter the wave-front to become high resolution, and regions which leave to return to low resolution. Regions in the far-field of the wavefront where the juvenile stage's minimum population density has crossed the T HR threshold become high resolution. Regions in the wave-back which have either a juvenile population that has crossed the T LR threshold, or densities which have changed less than a further parameter T LR;2 , become low resolution. As regions with only zero intrinsic population growth rates will always have zero population, they will always remain low resolution (making the algorithm very efficient for sparse landscapes).
The algorithm was implemented in C++. An executable which runs on Windows operating systems, is downloadable from GitHub (Gilbert et al. 2016); source code for non-Windows machines is available upon request.
We measured the algorithm's accuracy in simulating a non-stagestructured spreading population compared to a non-adaptive simulation algorithm (Section 'Accuracy'), its efficiency in CPU time and RAM usage (Section 'Efficiency') and tested it by simulating hypothetical biological invasions into coniferous woodland with (i) a single demographic stage across Great Britain, and (ii) two demographic stages into the New Forest (Hampshire, UK). In both cases, we used fine-scale mapping data from the CEH Land Cover Map 2007 (LCM2007) in Section 'Examples'. We also demonstrated the accuracy of our algorithm for stage-structured populations in Appendix S5.
Throughout the following section, we simulate IDEs of the form (1) with a single demographic stage and the commonly used 2D Laplace (exponential) dispersal kernel (Carrasco et al. 2010) with spatially constant dispersal parameter a, given by kðx À y; yÞ ¼ 1 2pa 2 exp Àkx À yk a ; eqn 2 where ║xÀy║ is the distance between y and x. We choose a piece-wise linear scalar population growth function due to its simplicity. Other stage-structured population projection matrices and single-stage growth rates, including those with Allee effects, can also be used. We have tested the algorithm for a growth function with an Allee effect and found that the simulations were accurate for selected examples. For the piece-wise linear growth function, the population growth rate is constant below the carrying capacity C, with the population after the growth phase at x given by where r 0 ðxÞ is the intrinsic population growth rate at x. Without loss of generality, we take the carrying capacity C = 1 throughout this section. In Sections 'Accuracy' and 'Efficiency' we take a = 1, but in Section 'Examples' we simulate over a large real landscape, and take a = 50 m to allow wider dispersal.
When measuring the algorithm's accuracy (Section 'Accuracy') and efficiency (Section 'Efficiency'), we take the initial distribution to be along one side of the square domain. For the example of spread of a specialist into a fragmented coniferous habitat in the UK as given by the LCM 2007 (Section 'Examples'), we took the initial distribution to be localised at a single cell.

A C C U R A C Y
We tested the accuracy of the adaptive algorithm presented in Section 'Materials and methods' by comparing its output with a high-resolution non-adaptive algorithm across a range of parameters and landscape types. We varied four userdefined tolerance parameters to investigate the trade-off between accuracy and speed, as well as two model parameters to determine how the demography and landscape pattern affect the algorithm's accuracy. For a spatially homogeneous landscape and the spatially periodic linear and island landscapes presented in Fig. 1b,c, we simulated both algorithms while varying six different parameters to explore the relationship among the parameter value, the landscape type and the algorithm's accuracy. We varied: the density threshold at which a region returns to low resolution once in the wave-back T LR , the density threshold at which a region in the wave-front becomes high resolution T HR , the cell size dx, the region size Dx, the growth rate in suitable habitat r and the landscape period p (the spatial distance taken for landscape features to repeat).
We measured the accuracy using the area occupied (Appendix S3) and the furthest distance reached by the population (Appendix S4), and found similar responses to the different parameters in both cases. For standard and default parameter values, the algorithm gave accurate results for the scenarios considered. For all three landscape types, simulation errors were minimised for smaller values of dx, larger values of T LR and r and intermediate values of Dx. The relationship between T HR and error differed between landscape scenarios and there was no straightforward relationship between p and simulation error. To show that our algorithm gives accurate results for stage-structured populations, we tested it against a stage-structured example and present our results in Appendix S5.

E F F I C I E N C Y
The adaptive algorithm was tested for efficiency against a nonadaptive algorithm. We varied the parameters that have the greatest effect on the algorithm's computational complexity as these demonstrated the differences in performance between the adaptive and non-adaptive algorithms. Different values to the ones used to test the accuracy were used as performance gains are only possible for larger simulations (see Table 1). We varied the landscape/domain length, L, the number of generations, T, and the ratio M of the region width Dx to the cell width dx, M = Dx/dx (Fig. 2).
For each of the three parameters varied, 24 values were investigated under the three landscape scenarios introduced in Section 'Accuracy'. For each parameter value and landscape scenario, the average of 10 runs was taken. The following performance measures were used: 1. CPU time. The total time used by the processor on the simulation. This ignores any interruptions from other tasks, so is generally less than the actual time taken for the simulation, and was measured via the UNIX command time(). 2. Maximum RAM is measured by the maximum resident set size over the programme's run that is the amount of main memory (RAM) occupied by the process.
In Fig. 2, we plot the adaptive algorithm's performance against the performance of the non-adaptive algorithm for the three different landscape scenarios and for the three different parameters. In all cases, the adaptive algorithm was generally faster than the non-adaptive algorithm, particularly as the domain length L, total number of generations, T, and the width in cells of each region M increased, with the adaptive algorithm using up to 10 times fewer resources than the non-adaptive algorithm (Fig. 2a,c,e). The adaptive algorithm also achieved a lower maximum RAM usage than the non-adaptive algorithm, especially as the domain length L and the ratio of region width to cell width M = Dx/dx increased.

E X A M P L E S
To run the algorithm, the landscape map must be saved in the same directory as the executable, and the model and algorithm parameters (which determine the speed and accuracy of the algorithm) must be set in the parameter text file. The parameter file includes the option of specifying the location of a point release for the spreading population, or of defining an initial population through another map file. The executable and user guide are available from GitHub (Gilbert et al. 2016).

UK land cover map
To assess our algorithm's performance on high-resolution maps of large landscapes, we simulated an invasion of a hypothetical specialist non-stage-structured species that can only inhabit coniferous woodland habitat using data from the 25-m resolution CEH Land Cover Map, LCM2007, covering the whole of Britain (Morton et al. 2011). LCM2007 has 1Á 456 Â 10 9 individual cells and our ability to simulate the model over this scale demonstrates the usefulness of the algorithm (the maximum RAM used by the adaptive algorithm is seven times less than the non-adaptive).
We present the results of this simulation in Fig. 3, with the parameter values: growth rate r Ã ¼ 200 in coniferous habitat and zero elsewhere (non-coniferous terrain and sea/ocean), dispersal parameter a = 50 m, region size Dx = 2Á5 km, landscape length L 1 ¼ 1300 km and landscape width L 2 ¼ 700 km. All other model parameters took the values tabulated under default accuracy values in Table 1.

Stage-structured example
To demonstrate the opportunities our algorithm offers for modelling stage-structured spread that our algorithm offers, we simulated the spread of a population with two demographic stages (juvenile and adult) across part of the New Forest (Hampshire, UK). The hypothetical species has the same demographic and dispersal parameters as the example in Appendix S5, except its population projection matrix is given by  Table 1.
We present the results of this simulation in Fig. 4, which show the areas occupied by juveniles (Fig. 4a) and adults (Fig. 4b) at each time step.

Discussion
Spatially explicit, mechanistic population models facilitate a quantitative understanding of changing population distributions that incorporate dispersal and demographic behaviour. These models enable predictions of population spread for any species and landscape fitting the assumptions of the model. Under a limiting set of assumptions, spread rates of spatial models can be calculated analytically (Shigesada, Kawasaki & Teramoto 1986;Kot, Lewis & den Driessche 1996;Neubert & Caswell 2000), but in general, numerical simulation must be used to predict the dynamics of changing population distributions in realistic scenarios, particularly when landscapes are either inhomogeneous or aperiodic. For all spatially explicit models, the accuracy of results obtained through simulation depends on the quality of the data used as input to the model and on the resolution of the discretisation of the landscape, with more grid points providing greater accuracy, but requiring more computational resources.
In this study we have developed an adaptive mesh method which facilitates the simulation of IDE models of spatial spread efficiently and without significant losses in accuracy. This method allows ecologists to simulate spatial dynamics for cases that have previously been too computationally demanding even on computer clusters, and thus permits previously unachievable research questions to be addressed, and allows researchers to get more out of limited computational resources. In particular, we show that it is possible to simulate spread straightforwardly over very large areas with appropriate resolution. This adaptive algorithm has the additional cost of a mesh refinement phase, but this is significantly smaller than the savings in computational resources in the growth and dispersal phases. We implemented this algorithm in C++ and have made an executable available on GitHub (Gilbert et al. 2016).
For a range of parameter values and landscape types, we found that the adaptive algorithm was accurate and that for IDEs on large and/or high-resolution domains, in the tests (Section 'Accuracy' and 'Efficiency'), the adaptive algorithm used up to 10 times less CPU time and maximum RAM than the non-adaptive algorithm. The algorithm was also applied to stage-structured IDE (Appendix S5), and was found to be highly accurate.
The loss of resolution in the far-field and wave-back of the invading population introduces some error. In particular, if T LR is too low, the errors in the wave-back may affect the spreading behaviour and if T HR is too high, large errors will appear in the wave-front that may significantly affect spread. These effects are further exacerbated by landscape heterogeneity, which affects the spreading population. As with most numerical solvers (Berger & Oliger 1984), there is a trade-off between model accuracy and efficiency, and the simulation parameters must be chosen to reduce both computational resources and the errors added by the numerical scheme. The acceptable level of error will depend on the application, but when calculating wave speeds an error of a few per cent will usually be acceptable. We suggest that users use the parameter values used in our examples (Table 1), and then lower the error tolerance and see if the solution significantly changes. If a significant change occurs, then a lower error tolerance is probably required.
To further demonstrate the algorithm's efficiency, we simulated a hypothetical invasive species' spread across coniferous woodland in Great Britain using the CEH LCM2007 habitat map of the UK. Using a non-adaptive algorithm to simulate over the whole of Great Britain at a resolution of 25 9 25 m 2 would require around 300 GB of RAM, which is not presently possible on a regular desktop computer. For the same problem, the adaptive algorithm used a maximum of 40 GB RAM and therefore enables the running of simulations that would have been impossible with non-adaptive algorithms.
In particular, the adaptive algorithm enables the simulation of the spread of species over large landscapes with variation in demographic parameters or complex fragmentation patterns. Spatial variation in demography is commonly observed (Oostermeijer et al. 1996), with environmental variation having a strong relationship to parameter variation (Bullock, Silvertown & Hill 1996). Fine-scale variation has been identified as a strong determinant of spread and persistence (Carrasco et al. 2010;Bennie et al. 2013). Therefore, incorporating highresolution spatial variation into large-scale models will enable more accurate predictions as well as greater insight into the effectiveness of fine-scale ecological interventions such as landscape buffers and wildlife corridors in preventing and facilitating spread. Indeed, models of spread along corridors, stepping stones, etc. in realistic landscapes have, to date, been limited with simplifying assumptions about spatial extent, structure and dynamics (Moilanen 2011;Fennell et al. 2012;Hodgson et al. 2012). Our method will facilitate more complex modelling over scales and resolutions appropriate for conservation policy and planning.
A natural extension to the method presented in this study is the inclusion of dispersal heterogeneity (Fig. 1e) and densitydependent dispersal, which is required to simulate organisms which settle based on habitat quality or population density as well as distance (e.g. Bonte, Van Dyck & Bullock 2012). This is more computationally intensive to simulate, and simulating these IDEs adaptively would provide even greater relative improvements in speed and efficiency than for homogeneous IDEs. Without dispersal homogeneity, the integral in (1) cannot be calculated using FFTs. Instead dispersal is normally computed by multiplying the vector of discretised population densities by an N 9 N matrix, where N is the number of points in the discretisation, requiring OðN 2 Þ floating point operations. An adaptive algorithm for simulating IDEs with heterogeneous dispersal would reduce the size of the matrix toN 2 , whereN is the number of cells in regions with high-resolution added to the number of regions with low-resolution information, leading to even larger reductions in computational resources than for the homogeneous cases presented in this study. Extending the algorithm to heterogeneous and densitydependent dispersal would vastly extend the types of models that could be simulated, and is for future development.