Habitat fragmentation: Simple models for local persistence and the spread of invasive species

Journal of Theoretical Biology 310 (2012) 231–238

Contents lists available at SciVerse ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Habitat fragmentation: Simple models for local persistence and the spread of invasive species R.G. Brown a,n, A.F. James a, J.W. Pitchford b, M.J. Plank a a b

Department of Mathematics and Statistics, University of Canterbury, Private Bag 4800, Christchurch 8140, New Zealand Department of Mathematics, University of York, York Y010 5DD, United Kingdom

H I G H L I G H T S c c c

Population spread speed and persistence in a randomly fragmented environment. Continuum limit of stochastic process allows analytical insight into spatial phenomena. Implications for control of invasive species.

a r t i c l e i n f o

abstract

Article history: Received 11 January 2012 Received in revised form 25 June 2012 Accepted 26 June 2012 Available online 10 July 2012

Understanding the persistence and growth of natural populations in environments subject to random localised change is relevant both to the conservation of threatened species and to the control of invasive species. By developing and analysing simple strategic growth models in environments subject to random fragmentation events, we show that simple approximations can be used to predict invasion speeds and extinction probabilities. The rate and size of fragmentation events interact in a nonlinear way, a finding with important consequences for the efficient control of invasive species. Infrequent, large-scale fragmentation events provide more effective means of control than more frequent, smaller scale efforts. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Stochastic model Invasion Control Spatial model

1. Introduction Biodiversity conservation is apparently simple: threatened species require environments where their extinction is unlikely, and invasive species must be prevented from colonising such habitats. This oversimplification masks important environmental, intra- and interspecific factors. Nevertheless, the fact remains that biodiversity depends on local persistence and on resisting invasion. We develop simple mathematical models to quantify these key criteria in the context of natural or anthropogenic habitat fragmentation. The work is motivated by challenges surrounding the conservation of native flora, but the messages apply more widely to species whose growth can be halted by random habitat disturbance. Natural populations observed in the wild are seldom homogeneous, more commonly occurring in patches both spatially and temporally. There is an extensive body of literature describing the mechanisms responsible for patchiness, including

n

Corresponding author. Tel.: þ64 3 3642987. E-mail address: [email protected] (R.G. Brown).

0022-5193/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jtbi.2012.06.033

dispersal (Nathan and Muller-Landau, 2000; Lewis, 1997), intraand interspecific competition (Lehman and Tilman, 1997; Bolker and Pacala, 1999; Pacala and Levin, 1997), heterogeneous environments (Amarasekare and Possingham, 2001; Bergelson et al., 1993; Cornell and Ovaskainen, 2008; Dewhirst and Lutscher, 2009; Johst et al., 2002; Keymer et al., 2000; Kinezaki et al., 2003; North and Ovaskainen, 2007), and habitat disturbance (Shigesada and Kawasaki, 1997; Tilman et al., 1994). This study focuses on the latter mechanism. Environmental disturbance, or fragmentation, can arise in different ways. Natural phenomena, such as landslips and earthquakes, can destroy populations locally, as can human and animal behaviour. For example, deforestation or local application of herbicides can be considered as fragmentation events for plant species. The aim of this work is to develop a strategic understanding of population dynamics in an environment that exhibits spatially and temporally random fragmentation. We seek to evaluate the effects of fragmentation on population persistence, density and patchiness, and on invasion speed for an invading population. A further aim is to understand the extent to which space needs to be explicitly modelled in order to capture the essential population

232

R.G. Brown et al. / Journal of Theoretical Biology 310 (2012) 231–238

tracks the presence or absence of the species at a given location, rather than the density. This has the advantage that, other than the simple linear spread assumption, the model is independent of the details of dispersal and reproduction. For a detailed examination of the general conditions under which the linear assumption can be expected to hold, see Weinberger et al. (2002). Section 2 presents the stochastic model used in this paper. Section 3 develops the theory for an idealised model in one spatial dimension. Section 4 explores practical implications by testing the results against two-dimensional simulations and assessing the implications for efficient control of invasive species.

2. Stochastic model

Fig. 1. Snapshot of a two-dimensional system. Dark areas are inhabited, white areas are uninhabited.

dynamics caused by the fragmentation process. Fig. 1 depicts a snapshot from a simulated habitat, illustrating some of the challenges faced. In the simulation, the plant species grows vegetatively at a constant rate and the environment is subject to a series of random fragmentation events (or biological control actions), each of which destroy all plants within a fixed radius (see Sections 2 and 3 for details). The full simulation is available in video form as an online supplement to this article. Fig. 1 shows that the emergent spatial structures, even for highly simplified model scenarios, may be complex and difficult to quantify. However, from a management perspective, the intricacies of these spatial structures are less important than the answers to simple questions, such as 1. Under what circumstances can a plant species persist in an environment subject to fragmentation? 2. How rapidly will an invading plant species spread? 3. Can this spread be halted? 4. Faced with finite resources, what is the most efficient strategy for controlling an invasive species? To provide strategic answers to these questions, we develop the simplest possible model that can capture essential growth and fragmentation characteristics. We thereby derive a mean-field approximation that allows us to quantitatively estimate equilibria and, perhaps more importantly, to identify changes in large-scale behaviour in response to changes of the environmental parameters. Early mathematical models of invasion, for example the density-dependent Fisher–Kolmogorov model (Fisher, 1937; Kolmogorov et al., 1937) and density-independent Skellam model (Skellam, 1951), predicted a linear spread rate. The spread rate can be related to the rate of dispersal and the intrinsic rate of increase of the species (Murray, 2002). It is now well known that some invasions spread linearly, whereas others do not (Hastings et al., 2005). For example, in a review of 14 historical plant invasions (Weber, 1998), half of the invasions exhibited linear spread and half did not. Mathematical models have shown that nonlinear spread can be caused by a variety of mechanisms, for example heterogeneous environments (Kinezaki et al., 2003; Dewhirst and Lutscher, 2009), Allee effects (Lewis and Kareiva, 1993), interspecific competition (Hosono, 1998), and long-range dispersal (Kot et al., 1996). The invasion model presented here assumes that the spread is linear. The model is based on occupancy, meaning that it only

We develop a two-state model based on occupancy. The population exists in a continuous one- or two-dimensional environment R with periodic boundary conditions. At any instant, t, the environment can be partitioned into points that are occupied, OðtÞ, and points that are not occupied, R\OðtÞ. It is assumed that OðtÞ is the disjoint union of a finite number of simply connected sets, i.e. the population exists as a finite number of locally homogeneous patches. (In a two-dimensional environment it is also assumed that the boundary of each patch is piecewise smooth.) The total occupied area at time t is simply the area of OðtÞ: Z UðtÞ ¼ dA: OðtÞ

The dynamics are governed by two simple mechanisms. Firstly, the population grows by spreading at a constant speed v normally outward from the boundary of each patch. Growing patches may therefore collide and merge with other patches. Secondly, the population can be reduced by random fragmentation events which occur as a Poisson process in time and space with a constant rate of l events per unit time per unit area. When a fragmentation event occurs, all points within a distance r of the event location become instantaneously unoccupied. The parameter r is referred to as the fragmentation radius. The set of occupied points at time t þ dt, where dt is a sufficiently small time step, can therefore be written to first order as     Oðt þ dtÞ ¼ OðtÞ [ x A R inf 9xy9 r vdt \BðldtÞDðX,rÞ, ð1Þ y A OðtÞ

where B(P) is a Bernoulli random variable that is equal to 1 with probability P, 9  9 denotes the Euclidean norm, X A R is a uniform random variable corresponding to the centre of the fragmentation event and DðX,rÞ is the disc of radius r and centre X.

3. One-dimensional model In one dimension, the stochastic process can be simulated exactly using an event-based algorithm (see Appendix). A time history of a one-dimensional simulation is shown in Fig. 2. 3.1. Spatial spread rate of the population From a management perspective, the spatial extent of the population is often more important than the total occupied area. Suppose that at time t the population is contained within the interval ½a,b, a subset of the whole environment. Note that the population is unlikely to fill the interval ½a,b contiguously; typically there will be fragmentation-induced gaps in this interval. Clearly, in the absence of fragmentation, the spatial extent of the population will increase from each end of the interval ½a,b with speed v. We ask how, in this simplest case, fragmentation acts to slow or even halt this growth.

R.G. Brown et al. / Journal of Theoretical Biology 310 (2012) 231–238

The predicted spread rate v was tested empirically by numerical experiments in which the population was initialised as a small contiguous interval in the centre of the environment and the spatial extent of the population (largest distance between any two occupied points) measured after a fixed period of time T0 (Fig. 3). The assumption used to derive v , namely that a fragmentation event only removes a portion of the endmost patch and not the entire patch, breaks down near v ¼ 0. This is the reason for the discrepancy between the stochastic process and the theoretical prediction as l approaches 100.

1

%

0.8 0.6 x

%

0.4

%

0.2

0

50 t

100

50 t

100

3.2. Mean-field approximation to the one-dimensional model To make further progress, we approximate the one-dimensional version of the stochastic model by a system of two ordinary differential equations, one for the total number n(t) of distinct patches (i.e. the number of disjoint sets in OðtÞ) and one for the total occupied area u(t). In this mean-field approximation it is assumed that the environment is sufficiently large that n(t) may be treated as a continuous variable. In a small time increment dt, each patch spreads from each of its two boundaries a distance vdt. The increase in the occupied area is therefore dugrowth ¼ 2vdt. Because fragmentation events are uniformly random in space, all occupied points are equally likely to become fragmented in short time dt. The expected decrease in the occupied area due to fragmentation is therefore dufrag ¼ 2r lu. Taking the limit as dt-0 gives the following equation for the rate of change of occupied area:

Population density

1 0.8 0.6 0.4 0.2 0

0

u_ ¼ 2vn2r lu:

Fig. 2. Time-evolution of a one-dimensional environment. (a) The simulated environment against time: grey represents occupied areas; black represents unoccupied areas. Gaps caused by fragmentation events close at a linear rate and therefore appear as triangles in the diagram. (b) The corresponding total area occupied against time, indicating approximately linear initial growth toward an asymptote as quantified in the subsequent analysis. Parameter values: L ¼1, r ¼0.05, v ¼ 8:33  103 , l ¼ 1.

Consider the right-hand boundary of the population xb(t) during the time interval ½t,t þ dt. If no fragmentation event occurs within a distance r of xb(t) during this time interval, then the boundary simply extends by a distance vdt. The probability of no fragmentation event occurring within this (moving) spatial interval ½xb ðtÞr, xb ðtÞ þr is 12r ldt þ Oðdt 2 Þ. Likewise, the probability of a single such fragmentation event occurring is 2r ldt þOðdt 2 Þ and the probability of more than one such event occurring is Oðdt 2 Þ. The expected amount of population removed by such an event is r (there is a simple linear relationship between fragmentation location and the amount of population removed: fragmentation events centred near xb r remove approximately 2r, those near xb þ r remove a negligible amount, and the fragmentation location is uniformly random.) Hence, the expected change in the location of the right-hand boundary in the time interval ½t,t þ dt is 2

2

2

Eðxb ðt þ dtÞxb ðtÞÞ ¼ vdtð12r ldt þ Oðdt ÞÞrð2rldt þOðdt ÞÞ þ Oðdt Þ:

Taking the limit as dt-0 gives the expected spread rate v at each boundary as

ð3Þ

The dynamics of the patch density n are more complicated. In one dimension, there are three ways in which the number of patches can change: patches can be annihilated by a fragmentation event; patches can split into two separate patches by a fragmentation event; and patches can merge with other patches through spread. Therefore, the equation for n_ takes the form n_ ¼ splittingmergingannihilation: The rates of these three processes depend on the distribution of the sizes of occupied patches and of unoccupied gaps at time t. We denote these distributions by probability density functions pðs,tÞ and qðs,tÞ, respectively. Annihilation: Patch annihilation can only happen to patches of size s r 2r. For a patch of size s to be annihilated, the

x 10−3 5

mean stochastic v*

4

speed

0

233

th

th

2.5 , 97.5 percentiles

3 2

%

v ¼ lim %

dt-0

Eðxb ðt þ dtÞxb ðtÞÞ ¼ v2r 2 l: dt

ð2Þ

This very simple result already suggests potentially important practical consequences. An invasive spread of vegetatively propagating plants is more readily halted by a few large fragmentation events (small l, large r) than by more frequent, smaller events (large l, small r). This is, however, a highly simplified argument which ignores the structure evident in Fig. 1.

1 0

0

50

100

150

λ Fig. 3. Spatial spread of the population over a fixed period of time T 0 ¼ 90. The dashed line shows the average of 50 simulations (for each value of l); the grey lines show the 2.5th and 97.5th percentiles. The solid line shows the theoretical prediction (Eq. (2)). Other parameter values: v ¼ 5  103 , r ¼ 5  103 , L¼ 1. The initial population was the interval ½0:45,0:55.

234

R.G. Brown et al. / Journal of Theoretical Biology 310 (2012) 231–238

fragmentation event must happen in a region of size 2rs. In time dt the probability of such an event happening to the patch is ð2rsÞldt þOðdt 2 Þ. The overall annihilation rate at time t is therefore given by Z 2r annihilation ¼ nl pðs,tÞð2rsÞ ds: 0

Splitting: Patch splitting can only happen to patches of size s 42r. For a patch of size s to be split, a fragmentation event must occur in a region of size s2r. In time dt the probability of such an event happening to the patch is ðs2rÞldt þOðdt 2 Þ. The overall splitting rate is therefore given by Z 1 pðs,tÞðs2rÞ ds: splitting ¼ nl 2r

Merging: The merging rate may be obtained by considering the expected number of patches merging in time dt. Two occupied patches merge when an unoccupied gap is closed up by spread of the two patches. Consider a gap of size s r 2r at time t ¼ t 0 . This gap will close at time t ¼ t 0 þ s=ð2vÞ, unless a fragmentation event occurs within a distance r of the boundary of the gap before this time. Fragmentation events occurring within a distance r of this narrowing interval occur as a nonhomogeneous Poisson process with rate m^ ðtÞ ¼ lð2r þ s2vtÞ, where t ¼ tt 0 A ½0,s=ð2vÞ. Hence, the number of such events occurring in the vicinity of the gap during the time interval ½t 0 ,t 0 þs=ð2vÞ follows a Poisson distribution with mean Z s=ð2vÞ lsð4r þsÞ : ð4Þ mðsÞ ¼ l ð2r þ s2vtÞ dt ¼ 4v 0 The probability that a gap of size s closes without such a fragmentation event occurring is emðsÞ . Now consider the entire environment. The merging rate is defined to be Eðnumber of gaps closing in ½t,t þ dtÞ : dt dt-0

During ½t,t þ dt, only gaps of size s r 2vdt can close, and provided 2vdt o2r no new gaps in this size range can be created. Hence, for sufficiently small dt, we have Eðnumber of gaps closing in ½t,t þ dtÞ Z 2vdt qðs,tÞemðsÞ ds ¼n ¼ 2vndtqðz,tÞemðzÞ for some 0 rz r 2vdt by the first mean value theorem for integration, assuming that qðs,tÞ is continuous with respect to s for 0 r s r2vdt. As dt-0, then necessarily z-0. By equation (4), mð0Þ ¼ 0, so the merging rate is given by merging ¼ 2vnqð0,tÞ: This is intuitive as the merging rate is proportional to the density of gaps of infinitesimally small size (nqð0,tÞ) and to the spread rate (v). Combining the annihilation, splitting and merging terms gives _ the equation for n: Z 1 n_ ¼ 2vnqð0,tÞ þ ln pðs,tÞðs2rÞ ds: R1 0

vn2 : 1u

This expression for the merging rate implies that qð0,tÞ ¼

n , 2ð1uÞ

ð7Þ

i.e. the reciprocal of twice the mean gap length. In the limiting case where fragmentation events are vanishingly unlikely to overlap spatially, the distribution of gap sizes is uniform on ½0; 2r, which implies that qð0,tÞ ¼ 1=ð2rÞ. This is consistent with (7) as the unoccupied area (1u) in this special case is simply equal to the number of gaps (n) multiplied by the mean gap size (r). The dynamical system approximating the evolution of patches and gaps over time is therefore 8 < u_ ¼ 2vn2r lu, ð8Þ vn2 : n_ ¼  þ lu2r ln: 1u This system has two equilibria, the trivial equilibrium ðu ,n Þ ¼ ð0; 0Þ and   v2r 2 l rl : ð9Þ ðu ,n Þ ¼ 1, v vr 2 l %

spðsÞ ds must equal u=n gives the

u_ ¼ 2vn2r lu,

ð5Þ

n_ ¼ 2vnqð0,tÞ þ lu2r ln:

ð6Þ

Because a closed-form expression for qð0,tÞ is, in general, intractable, we approximate the merging rate by considering

%

%

If v 4 2r 2 l, this equilibrium is stable and the origin is unstable, so the population is expected to converge to persist in the long-term. At v ¼ 2r 2 l the system undergoes a transcritical bifurcation, and the two equilibria interchange stability. If v o2r 2 l, then the origin is stable and the equilibrium in Eq. (9) is not in the biologically relevant region (u ,n Z 0), so the population will go extinct. Note that this bifurcation point relates directly to the criterion for population spread derived in Section 3.1: the condition for a stable positive equilibrium is equivalent to that for the local spread of the population through unoccupied habitat. Fig. 4 shows the comparison between simulations of the stochastic process and the solution of the mean-field approximation derived above. The results are shown for stable persistence case (Fig. 4a, v 4 2r 2 l), the borderline case (Fig. 4b, v ¼ 2r 2 l) and the extinction case (Fig. 4c, v o 2r 2 l). Close to the extinction– persistence threshold, the mean-field approximation is less accurate, but exhibits the same qualitative behaviour over short to medium time scales (Fig. 4c). Over longer time scales, extinction necessarily becomes more likely in the stochastic model, magnifying discrepancies between stochastic realisations and the meanfield predictions. Fig. 5 shows the predictions of the stochastic model and meanfield approximation for a range of fragmentation rates l. The discrepancy between the solid and dashed curves is due to the frequent extinction of the stochastic model for values of l significantly lower than the mean-field extinction threshold of l ¼ 250. This type of discrepancy is a well-known feature of meanfield approximations. For example, the Kermack–McKendrick model of an infectious disease predicts an epidemic whenever the basic reproduction number (R0) is greater than 1, but the equivalent stochastic model will frequently fail to undergo an %

0

0

merging ¼

%

merging ¼ lim

Noting that mean patch size system

each patch as a ‘‘predator’’ encountering other patches at some rate to be determined. The effective search rate, or speed, of an individual patch is 2v, as the patches it encounters are moving directly toward it at speed v. The encounter rate for a particular patch is proportional to the ‘‘prey’’ density n=ð1uÞ, the reciprocal of the mean distance between patches. The individual encounter rate is thus given by 2vn=ð1uÞ. The population encounter rate is attained by multiplying by n, and then dividing by 2 to avoid counting each patch–patch encounter twice. This gives the population-level merging rate as

%

R.G. Brown et al. / Journal of Theoretical Biology 310 (2012) 231–238

u

0.8 0.6

300

1

200

240

0.8

160

180

0.6

u

n

0.2 0

0

5

10 t

n

n

0.4

120

0.4

60

0.2

80 40

u

0

0 20

15

120

n

u

1

235

0

20

40 t

1

200

0.8

160

0.6

60

0 80

120

u

n

n 0.4

80

0.2

40 u

0

0

1

2

3

4

5

6

0

t Fig. 4. Simulations of the stochastic model (solid) and solutions of the mean-field approximation (dashed): (a) in the stable persistence case; (b) on the borderline between persistence and extinction; (c) in the extinction case. Parameter values: r ¼ 103 , v ¼ 5  104 , L¼1.

1

200

0.8

160

u

0.6

120

u

n

n 0.4

80

0.2 0

40

0

50

100

150

200

0 250

λ Fig. 5. Equilibrium values of occupied area (u) and number of patches (n). Equilibrium values for the stochastic process (solid) were obtained by running a single realisation for sufficiently long (T¼300) that any initial transients had died out, and averaging 50 evenly spaced samples between T ¼300 and T¼ 400. The dashed line shows the equilibrium values predicted by the mean-field approximation. Parameter values: r ¼ 103 , v ¼ 5  104 , L ¼1.

epidemic when R0 is close to 1 (Diekmann and Heesterbeek, 2000). As seen in Fig. 5, predictions of the mean-field approximation are accurate away from the extinction threshold.

first formulate suitable nondimensional parameters that allow direct comparison between the two domains. The parameters of the one-dimensional model are the fragmentation event rate l, the fragmentation radius r, the local spreading speed v and the environment size L. We define the fragmentation size r to be the proportion of the environment affected by a fragmentation event:

r¼

2r ð1DÞ, L

r¼

pr 2 L2

ð2DÞ:

ð10Þ

To relate the spreading speeds of the two models, a meaningful parameter that is a function of v needs to be determined. We consider the limiting scenario where fragmentation events are rare and hence non-overlapping, and prescribe the spreading speeds v such that the expected equilibrium population sizes are equal. Under this assumption, the equilibrium population size is given by  mean gap area   ðmean number of gapsÞ u ¼ 1 area environment mean gap area   ðgap lifetimeÞ  ðevent rate per unit timeÞ: ¼ 1 environment area %

In both one and two dimensions, the gap lifetime (i.e. the time taken for a gap to fill) is given by r=v. So ! Z r  1 v r=v r2 l ðlLÞ ¼ 1 ð1DÞ 2ðrvtÞ dt u ¼ 1 L r t¼0 v v %

and 4. Two-dimensional model u ¼ 1 %

Extending the model into two spatial dimensions, while reducing its analytical tractability, allows more practical insight into the spread and control of invasive species. In order to ascertain whether the conclusions drawn from the one-dimensional results carry over to two dimensions, it is necessary to compare one- and two-dimensional simulations using similar parameter ranges. The parameters of the one-dimensional model do not translate directly into a two-dimensional setting, so we

1 L2

v r

Z

r=v

! r  pr3 l ðlL2 Þ ¼ 1 ð2DÞ: v 3v

pðrvtÞ2 dt

t¼0

ð11Þ

Comparing these two expressions, a suitable dimensionless parameter to capture the ‘‘speed’’, Z, is defined to be

Z¼

v ð1DÞ, r2 l

Z¼

3v

pr3 l

ð2DÞ:

ð12Þ

Ensuring that the values of r, Z, and l are the same in the two forms of the model enables the one- and two-dimensional results

236

R.G. Brown et al. / Journal of Theoretical Biology 310 (2012) 231–238

to be directly compared. Note that the equilibrium of the system (8) can also be recast in these parameters as   Z2 2 1, : ð13Þ ðu ,n Þ ¼ Z1 Zr %

where A(t) is the area of the convex hull at time t, gives an estimate of the linear extent of the population. The effective spread speed is thus defined to be

%

r est r 0 1 ¼ ð T0 T0

vest ¼ %

ð14Þ

Fig. 8 shows a contour plot of the effective spread speed vest against l and r. Unsurprisingly, an increase in either l or r in isolation reduces vest , and hence increases the probability of extinction. However, these factors do not contribute equally. The superimposed curves in Fig. 8 represent curves of equal control effort (constant lr 2 ). The curves of constant effort cross the spread speed contours, showing that control strategies with equal overall effort can lead to different outcomes. Strategies with large r and small l (top-left of Fig. 8) are more likely to result in %

%

1 0.8 0.6 u

A level-set method (Sethian, 1999) was implemented to simulate the stochastic process in two dimensions (see Appendix). A snapshot of a two-dimensional simulation (with parameter values r ¼ 0:1, Z ¼ 3:333, l ¼ 1 corresponding to the one-dimensional simulation shown in Fig. 2) is depicted in Fig. 1. Different sized round gaps can easily be made out, together with regions where gaps have overlapped. A movie of the full two-dimensional simulation is available as an online supplement to this article. To compare the behaviour of the two-dimensional stochastic model with the deterministic model, r, Z, and l were chosen to correspond to the one-dimensional simulations depicted in Figs. 4 and 5. The results are depicted in Figs. 6 and 7. Note that both the time dynamics and equilibrium prediction are again qualitatively similar to the deterministic model. As with the onedimensional simulations (Figs. 4 and 5), there is a discrepancy in the long-term behaviour near the persistence–extinction threshold as the stochastic system is more likely to go extinct. In the context of control, the amount of control effort is E ¼ pr 2 l. For fixed E, we have Z ¼ 3v=ðrEÞ, so Z can be increased (and hence equilibrium population density decreased) by decreasing r. This leads to the same conclusion as was drawn in Section 3.1 for the one-dimensional model: for a fixed amount of control effort, larger, less frequent fragmentation events (large r, small l) are more effective in reducing population density than smaller, more frequent ones (small r, large l). To test this hypothesis, we calculated the effective spread speed of a population for a range of values of the fragmentation rate l and radius r. This speed was measured by running a twodimensional simulation with an initial population consisting of a circle of radius r0. If we assume that the population spread is approximately isometric, then the convex hull of the resulting pffiffiffiffiffiffiffiffiffiffiffiffiffiffi population will be approximately circular, and r est ¼ AðtÞ=p,

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AðT 0 Þ=pir 0 Þ:

0.4 0.2 0

0

50

100

150 λ

200

250

300

Fig. 7. Equilibrium values of occupied area (u) for the 2D stochastic model and the deterministic 1D model (13) with corresponding parameters. Equilibrium values for the stochastic process (solid) were obtained by running a single realisation for sufficiently long (T ¼300) that any initial transients had died out, and computing the mean population density between T¼ 300 and T¼ 400. The dashed line shows the equilibrium values predicted by the mean-field approximation. Parameter values were chosen to give the same values of the dimensionless parameters, r and l, as Fig. 5 (r ¼ 0:002, Z ¼ 500=l): r ¼ 2:523  102 , v ¼ 8:410  104 , L¼ 1.

1

0.8

0.8

0.6

0.6

u

u

1

0.4

0.4

0.2

0.2

0

0

5

10 t

15

0

20

0

20

40 t

60

80

1 0.8 u

0.6 0.4 0.2 0

0

1

2

3 t

4

5

6

Fig. 6. Simulations of population density in the two-dimensional stochastic model (solid) and solutions of the mean-field approximation (dashed): (a) in the stable persistence case; (b) on the borderline between persistence and extinction; (c) in the extinction case. Parameter values were chosen to match those in Fig. 4: r ¼ 2:523  102 , v ¼ 8:410  103 , L¼ 1.

R.G. Brown et al. / Journal of Theoretical Biology 310 (2012) 231–238

x 10−3

0.1

4

0.09 r 2λ = 0.06

0.08

3.5 3

0.06

2.5

r

v*

0.07

r 2λ = 0.04

0.05

2

0.04

1.5 2

r λ = 0.02

0.03

1

0.02 0.01

0.5 0

10

20

30

40

50

0

λ Fig. 8. Effective population spread speed. Lighter colours correspond to a faster spreading population. In the presence of no fragmentation, v ¼ 5  103 . The white dashed lines correspond to constant effort r 2 l ¼ const. Parameter values: T 0 ¼ 80, r 0 ¼ 0:1, v ¼ 5  103 . %

extinction than strategies with small r and high l (bottom-right of Fig. 8). This confirms the hypothesis from the deterministic model that large but infrequent fragmentation events have a more severe impact on population density than small, frequent ones.

5. Discussion There is an array of mathematical models of population dynamics that incorporate space to different degrees. The extent to which space is incorporated is a trade off between capturing spatially dependent phenomenon while maintaining some analytical tractability. Even apparently non-spatial models such as the Ricker map or nonlinear predator–prey systems can be regarded as accounting for spatial effects implicitly, for example through the processes defining the reproductive map (Sumpter and Broomhead, 2001) or the functional form of the predator–prey encounter rates (Pitchford and Brindley, 2001). Spatially implicit models such as metapopulations (Hanski, 1999) or pair approximations and spatial moment models (Dieckmann et al., 2000) represent a step towards modelling spatial effects. In these models there are no explicit space variables, but neighbourhood-based interactions are captured by either grouping populations into local metapopulations or by modelling the distribution of distances between individuals. Finally, space can be explicitly included using, for example, stochastic individual-based models or cellular automata to mimic real landscapes (Balzter et al., 2009; Dytham, 1995) or mathematically idealised PDE representations (Murray, 2002; Kot, 2001). The approach developed here is most closely related to pair approximation and moment methods, although no assumptions on high moments or closure terms are involved, and space is incorporated without granularity. This allows us to be confident that our conclusions arise from the spatiotemporal structures induced by fragmentation and vegetative growth, rather than any hidden complexity or numerical artefact caused by the modelling strategy. We have shown that the spread of invasive species can be slowed down or halted by random environmental fragmentation events. Similarly, a sufficient level of habitat fragmentation can eradicate local populations. However, what is less obvious is the way in which the frequency and size of fragmentation events interact. What emerges from this study is a rule of thumb for

237

conservationists: species are likely to be more seriously affected by infrequent but large fragmentation events than by more frequent, smaller events. This rule of thumb applies both to the conservation of threatened species, where one seeks to mitigate the effects of natural and anthropogenic disturbance, and to the eradication of invasive species where the aim is to deploy limited resources as efficiently as possible. The notion of spatially random interventions to control the spread of invasive species may appear unrealistic. However, where control efforts involve the participation of many independent agents, the analogy is appropriate. For example, efforts to prevent the spread of alien plant species in New Zealand typically involve public engagement (Kriticos et al., 2004; Bourdˆot et al., 1992), where members of the public are encouraged to clear particular species from their locations. The result is effectively random (because not all members of the public respond) and has a characteristic spatial scale based on the amount of land a particular person can clear. Our results indicate that public engagement policies based on the coordinated clearances of larger areas are likely to be more effective, even if budgetary or other constraints dictate that such large scale efforts are carried out less frequently. Even for more realistic and complex management scenarios (Buckley et al., 2005, 2004; Shea, 1998), the simple findings of this model can inform the manner in which limited control resources are deployed. The model we have used assumes that the boundary of the region occupied by the plant species advances at a constant speed. Although many invasions do spread in this way (Skellam, 1951; Weber, 1998), this assumption is likely to be unrealistic in cases where there is a heterogeneous environment, long-range dispersal, or interactions between different species (Kot et al., 1996; Dewhirst and Lutscher, 2009; Weinberger et al., 2002). Nevertheless, such factors could be readily incorporated into the twodimensional stochastic model. A heterogeneous environment could be represented by a spatially varying function for the spread rate. Long-range dispersal could be modelled, via a suitable dispersal kernel, by the formation of new colonies, spatially separated from the existing population. A simple two-species model could assume that the species are spatially mutually exclusive and that each species advances into unoccupied space at its own spread rate. All of these model generalisations could be simulated by a level-set method (Sethian, 1999). Coexisting species and Allee effects, which can also be responsible for nonlinear spreading (Lewis and Kareiva, 1993), would require a different approach that includes population density, rather than simply occupancy.

Acknowledgments We are grateful to NZIMA for postdoctoral support (RGB), University of York pump-priming funds (RGB) and a University of Canterbury Erskine Fellowship (JWP).

Appendix A A.1. One-dimensional simulations Every patch is stored individually, and the algorithm jumps discretely in time from event to event. The algorithm proceeds as follows: 1. Set t ¼0, initialise patches, compute t next ¼ t þT: 2. Compute tmerge ¼(distance between two closest patches)/2v. 3. if t þ t merge ot next (a) extend each boundary of each patch by vtmerge , (b) merge the two closest patches,

238

R.G. Brown et al. / Journal of Theoretical Biology 310 (2012) 231–238

(c) set t ¼ t þ t merge . else (a) extend each boundary of each patch by vðt next tÞ, (b) process a fragmentation event at a location chosen randomly from Uð0,LÞ, (c) set t ¼ t next , (d) set t next ¼ t þ T. 4. Goto 2. In the above algorithm T  ExpðlÞ represents an exponentially distributed random variate with rate parameter l. A.2. Two-dimensional simulations In two dimensions, a level-set method can be used to track the population’s moving boundaries. The boundaries are described by the zero level set of a function f: a point (x,y) is occupied at time t if and only if fðx,y,tÞ Z 0. Since the growth of a population patch is assumed to be at speed v normal to its boundary, the function f satisfies the following PDE (Sethian, 1999):

ft ¼ vJrfJ:

ð15Þ

A fragmentation event of radius r at x0 may be represented by the following mapping of f:

f/minðf,gðx,x0 ,rÞÞ,

ð16Þ

where gðx,x0 ,rÞ is a function whose zero level set is a circle of radius r centred at x0. We use a cone: gðx,x0 ,rÞ ¼ Jxx0 Jr. We obtain numerical solutions of Eq. (15) using an explicit, upwind finite difference scheme (see Sethian, 1999, Chapter 6) on a rectangular domain with periodic boundary conditions. We choose the time step dt to be dt ¼ minðr=ð10vÞ,1=ð5lÞÞ so that it is small enough to capture individual fragmentation events and the motion of the patches. Fragmentation events are moved slightly to be aligned with time steps; if more than one event occurs during one time step they effectively take place simultaneously.

Appendix B. Supplementary material Supplementary data associated with this article can be found in the online version of http://dx.doi.org/10.1016/j.jtbi.2012.06.033.

References Amarasekare, P., Possingham, H., 2001. Patch dynamics and metapopulation theory: the case of successional species. J. Theor. Bio. 209 (3), 333–344. ¨ Balzter, H., Braun, P.W., Kohler, W., 2009. Cellular automata models for vegetation dynamics. Ecol. Modell. 107 (2–3), 113–125. Bergelson, J., Newman, J.A., Floresroux, E.M., 1993. Rate of weed spread in spatially heterogeneous environments. Ecology 74, 999–1011. Bolker, B.M., Pacala, S.W., 1999. Spatial moment equations for plant competition: understanding spatial strategies and the advantages of short dispersal. Am. Nat. 153 (6), 575–602. Bourdˆot, G., Hurrell, G., Saville, D., 1992. Eradication of nassella tussock (Nassella trichotoma), an unlikely outcome of grubbing. N.Z. J. Agric. Res. 35, 245–252. Buckley, Y.M., Rees, M., Paynter, Q., Lonsdale, M., 2004. Modelling integrated weed management of an invasive shrub in tropical Australia. J. Appl. Ecol. 41 (3), 547–560.

Buckley, Y.M., Brockerhoff, E., Langer, L., Ledgard, N., North, H., Rees, M., 2005. Slowing down a pine invasion despite uncertainty in demography and dispersal. J. Appl. Ecol. 42 (6), 1020–1030. Cornell, S.J., Ovaskainen, O., 2008. Exact asymptotic analysis for metapopulation dynamics on correlated dynamic landscapes. Theor. Popul. Biol. 74 (3), 209–225. Dewhirst, S., Lutscher, F., 2009. Dispersal in heterogeneous habitats: thresholds, spatial scales, and approximate rates of spread. Ecology 90 (5), 1338–1345. Dieckmann, U., Law, R., Metz, J.A.J., 2000. The Geometry of Ecological Interactions, vol. 1. Cambridge University Press. Diekmann, O., Heesterbeek, J.A.P., 2000. Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation. John Wiley & Sons. Dytham, C., 1995. The effect of habitat destruction pattern on species persistence: a cellular model. Oikos 74 (2), 340–344. Fisher, R.A., 1937. The wave of advance of advantageous genes. Ann. Eugen. 7, 355–369. Hanski, I., 1999. Metapopulation Ecology, vol. 399. Oxford University Press. Hastings, A., Cuddington, K., Davies, K.F., Dugaw, C.J., Elmendorf, S., Freestone, A., Harrison, S., Holland, M., Lambrinos, J., Malvadkar, U., Melbourne, B.a., Moore, K., Taylor, C., Thomson, D., 2005. The spatial spread of invasions: new developments in theory and evidence. Ecol. Lett. 8 (1), 91–101. Hosono, Y., 1998. The minimal speed of travelling fronts for a diffusive Lotka Volterra competition model. Bull. Math. Biol. 60, 435–448. Johst, K., Brandl, R., Eber, S., 2002. Metapopulation persistence in dynamic landscapes: the role of dispersal distance. Oikos 98 (2), 263–270. Keymer, J.E., Marquet, P.A., Velasco-Hernandez, J.X., Levin, S.A., 2000. Extinction thresholds and metapopulation persistence in dynamic landscapes. Am. Nat. 156 (5), 478–494. Kinezaki, N., Kawasaki, K., Takasu, F., Shigesada, N., 2003. Modeling biological invasions into periodically fragmented environments. Theor. Popul. Biol. 64 (3), 291–302. Kolmogorov, A., Petrovsky, I., Piscounov, N., 1937. e´tude de l’e´quation de la diffusion avec croissance de la quantite´ de matie re et son application a´ un proble me biologique. Moscow Univ. Bull. Math. 1, 1–25. Kot, M., 2001. Elements of Mathematical Ecology. Cambridge University Press. Kot, M., Lewis, M.A., van den Driessche, P., 1996. Dispersal data and the spread of invading organisms. Ecology 77, 2027–2042. Kriticos, D., Lamoureaux, S., Bourdˆot, G., Pettit, W., 2004. Nassella tussock: current and potential distributions in New Zealand. N.Z. Plant Prot. 57, 81. Lehman, C.L., Tilman, D., 1997. Competition in spatial habitats. In: Tilman, D., Kareiva, P. (Eds.), Spatial Ecology. Princeton University Press, pp. 185–203. Lewis, M.A., 1997. Variability, patchiness, and jump dispersal in the spread of an invading population. Biol. Invas. 30, 46–69. Lewis, M.A., Kareiva, P., 1993. Allee dynamics and the spread of invading organisms. Theor. Popul. Biol. 43, 141–158. Murray, J.D., 2002. Mathematical biology: I. An introduction. In: Antman, S.S., Marsden, J.E., Sirovich, L., Wiggins, S. (Eds.), Interdisciplinary Applied Mathematics, third ed., vol. 1. Springer. Nathan, R., Muller-Landau, H., 2000. Spatial patterns of seed dispersal, their determinants and consequences for recruitment. Trends Ecol. Evol. 15 (7), 278–285. North, A., Ovaskainen, O., 2007. Interactions between dispersal, competition, and landscape heterogeneity. Oikos 116 (7), 1106–1119. Pacala, S.W., Levin, S.A., 1997. Biologically generated spatial pattern and the coexistence of competing species. In: Tilman, D., Kareiva, P. (Eds.), Spatial Ecology the Role of Space in Population Dynamics and Interspecific Interactions, vol. 30. Princeton University Press, pp. 204–232. (Chapter 9). Pitchford, J.W., Brindley, J., 2001. Prey patchiness, predator survival and fish recruitment. Bull. Math. Biol. 63 (3), 527–546. Sethian, J.A., 1999. Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision and Materials Science. Cambridge University Press. Shea, K., 1998. Management of populations in conservation, harvesting and control. Trends Ecol. Evol. 13 (9), 371–375. Shigesada, N., Kawasaki, K., 1997. Biological Invasions: Theory and Practice. Oxford Series in Ecology and Evolution. Oxford University Press. Skellam, J.G., 1951. Random dispersal in theoretical populations. Biometrika 38, 196–218. Sumpter, D.J., Broomhead, D.S., 2001. Relating individual behaviour to population dynamics. Proc. R. Soc. B: Biol. Sci. 268 (1470), 925–932. Tilman, D., May, R.M., Lehman, C.L., Nowak, M.A., 1994. Habitat destruction and the extinction debt. Nature 371 (6492), pp. 65–66, Nature Publishing Group. Weber, E., 1998. The dynamics of plant invasions: a case study of three exotic goldenrod species (Solidago L.) in Europe. J. Biogeogr. 25 (1), 147–154. Weinberger, H.F., Lewis, M.A., Li, B., 2002. Analysis of linear determinancy for spread in cooperative models. J. Math. Biol. 45, 183–218.