Spatiotemporal complexity of patchy invasion in a predator-prey system with the Allee effect.

ARTICLE IN PRESS

Journal of Theoretical Biology 238 (2006) 18–35 www.elsevier.com/locate/yjtbi

Spatiotemporal complexity of patchy invasion in a predator-prey system with the Allee effect. Andrew Morozova, Sergei Petrovskiia,, Bai-Lian Lib a

Shirshov Institute of Oceanology, Russian Academy of Science, Nakhimovsky Prosp. 36, Moscow 117218, Russian Federation Ecological Complexity and Modeling Laboratory, Department of Botany and Plant Sciences, University of California, Riverside, CA 92521-0124, USA

b

Received 5 January 2005; received in revised form 29 April 2005; accepted 2 May 2005 Available online 6 July 2005

Abstract Invasion of an exotic species initiated by its local introduction is considered subject to predator–prey interactions and the Allee effect when the prey growth becomes negative for small values of the prey density. Mathematically, the system dynamics is described by two nonlinear diffusion–reaction equations in two spatial dimensions. Regimes of invasion are studied by means of extensive numerical simulations. We show that, in this system, along with well-known scenarios of species spread via propagation of continuous population fronts, there exists an essentially different invasion regime which we call a patchy invasion. In this regime, the species spreads over space via irregular motion and interaction of separate population patches without formation of any continuous front, the population density between the patches being nearly zero. We show that this type of the system dynamics corresponds to spatiotemporal chaos and calculate the dominant Lyapunov exponent. We then show that, surprisingly, in the regime of patchy invasion the spatially average prey density appears to be below the survival threshold. We also show that a variation of parameters can destroy this regime and either restore the usual invasion scenario via propagation of continuous fronts or brings the species to extinction; thus, the patchy spread can be qualified as the invasion at the edge of extinction. Finally, we discuss the implications of this phenomenon for invasive species management and control. r 2005 Elsevier Ltd. All rights reserved. Keywords: Biological invasion; Predator–prey system; Patchy spread; Allee effect; Spatiotemporal chaos

1. Introduction Understanding of patterns and mechanisms of species spatial dispersal is an issue of significant current interest in conservation biology and ecology. It arises from many ecological applications; in particular, it plays a major role in connection to biological invasion and epidemic spread (Drake et al., 1989; Hengeveld, 1989; Murray, 1989; Shigesada and Kawasaki, 1997). A variety of theoretical approaches has been developed Corresponding author. Current address: 114 Winchester Gardens, Birmingham B31 2QB, UK. Tel.:+44 121 258 0899. E-mail addresses: [email protected] (A. Morozov), [email protected] (S. Petrovskii), [email protected] (B.-L. Li).

0022-5193/$ - see front matter r 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jtbi.2005.05.021

and a considerable progress has been made during the last decade. However, many aspects related to species dispersal have never been properly addressed yet. Regarding the spread of exotic species, a problem of high practical importance is how to create an effective program of management and control of invasive species (Andow et al., 1990; Sakai et al., 2001; Fagan et al., 2002). Such a program must include both good understanding of the mechanisms underlying species spread and an optimal monitoring strategy. In its turn, a relevant strategy is likely to be different for different species (e.g. depending on whether a given species is detectable with satellite-based remote sensing imagery data) and should be also based on the knowledge of the pattern of spread. For instance, in many cases invasion

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

of exotic species takes place via propagation of a population front separating the areas where given species is absent, i.e. in front of the front, from the areas where it is present in considerable densities, i.e. in the wake of the front. In this case, under somewhat idealized assumption that the invasion is going isotropically so that its rate does not depend on the direction of spread, the advance of invading species can be monitored using a relatively small number of on-site observers or sampling stations situated along a certain line coming out of the place of original species introduction. However, reality is often much more complicated. As a result of environmental heterogeneity the rates of invasion can be significantly different in different directions (Shigesada and Kawasaki, 1997). Also, due to the impact of both environmental and biological factors, the pattern of spread can be more complicated than simple population front. There is growing evidence (Davis et al., 1998; Kolb et al., 2004; Swope et al., 2004) that, in some cases, invasion of exotic species takes place through dynamics of separate population patches not preceded by propagation of a continuous population front, see also, Shigesada and Kawasaki (1997), Lewis (2000); and Lewis and Pacala (2000) and the references therein. Below we will call that pattern of spread a patchy invasion. Obviously, in the case of patchy invasion an adequate monitoring strategy should be more complicated and likely include many more observers. Thus, distinguishing between the situations when species spread takes place via propagation of a continuous population front and when it happens via patchy invasion, as well as identification of factors enhancing or hampering patchy invasion are problems of significant practical and theoretical importance. The origin of patchy invasion is often seen either in environmental heterogeneity (cf. Murray, 1989) or in the impact of stochastic factors (Lewis, 2000; Lewis and Pacala, 2000). Indeed, the whole dynamics of ecological communities appears as a result of interplay between numerous deterministic and stochastic factors. However, the importance of stochasticity should not be overestimated. Sometimes stochastic models and deterministic models lead to qualitatively similar spatiotemporal patterns, although reached through different mechanisms (cf. Kawasaki et al., 1997; Mimura et al., 2000). For a rather general model of marine ecosystem, Malchow et al. (2002) showed that there exists a critical level of noise so that for undercritical noise the system is more driven by deterministic factors. Even in the supercritical case, when noise can change the system dynamics considerably, it was shown that intrinsic spatial scales of the system are still controlled by deterministic mechanisms (Malchow et al., 2004).

19

Recently, it was shown by Petrovskii et al. (2002a,b) that patchy invasion can arise in a fully deterministic predator–prey system as a result of the Allee effect, i.e. of a threshold phenomenon when the population growth rate becomes negative for low population density (Allee, 1938; Dennis, 1989; Courchamp et al., 1999). Deterministic patchy invasion was shown to correspond to the invasion at the edge of extinction (Petrovskii and Venturino, 2004; Petrovskii et al., 2005c) so that a small finite variation of the system parameters either restores usual population front propagation scenario or brings the species to extinction. Moreover, it was shown that the system dimensionality is a crucial point and the patchy invasion in two spatial dimensions corresponds to species extinction in the corresponding 1D system (Petrovskii and Venturino, 2004; Petrovskii et al., 2005c). The above papers, however, left many questions open. In this paper, we make a detailed study of the deterministic patchy invasion in a predator–prey system where the prey growth is damped by the Allee effect. In particular, the following issues are addressed: (i) what is the succession of invasion regimes in response to variation of an ecologically meaningful controlling parameter and how the invasion speed depends on the type of spread, (ii) what is the degree of spatiotemporal complexity corresponding to the regime of patchy spread and (iii) what are the ecological implications of the patchy invasion. We show that, although the invasion speed is much lower in the regime of patchy spread than it is in the usual regime(s) of continuous front propagation, the patchy spread provides a scenario of species invasion below the survival threshold. Also, we show that deterministic patchy invasion corresponds to spatiotemporal chaos and estimate the value of the dominant Lyapunov exponent.

2. Mathematical model We consider 2D dynamics of a predator–prey system described by two partial differential equations of diffusion–reaction type (Nisbet and Gurney, 1982; Murray, 1989; Holmes et al., 1994; Shigesada and Kawasaki, 1997; Medvinsky et al., 2002):
2  qP q P q2 P ¼ D1 þ þ F ðPÞ
f ðPÞZ, (1) qT qX 2 qY 2
2  qZ q Z q2 Z ¼ D2 þ þ kf ðPÞZ
MZ, qT qX 2 qY 2

(2)

Here P ¼ PðX ; Y ; TÞ and Z ¼ ZðX ; Y ; TÞ are densities of prey and predator, respectively, at moment T and position ðX ; Y Þ. The function F ðPÞ stands for the intrinsic prey growth, f ðPÞ is the predator trophic

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

20

response, k is the food utilization coefficient, D1 and D2 are diffusion coefficients and MZ describes predator mortality. We assume that the dynamics takes place in a homogeneous environment, which means that none of the parameters exhibit any dependence on space. In the most part of this paper, we consider D1 ¼ D2 ; the case of unequal diffusivities will be briefly addressed in the last section. We assume that the predator trophic response is of Holling type II and choose the following parameterization (Holling, 1965; Murray, 1989): P f ðPÞ ¼ A , (3) PþB where B is the half-saturation density and A is the maximum predation rate. Considering Eqs. (2) and (3) together, it is readily seen that the product Ak has the meaning of predator maximum growth rate. In the case of strong Allee effect, a convenient parameterization for the prey growth is as follows (cf. Lewis and Kareiva, 1993; Owen and Lewis, 2001): F ðPÞ ¼ mPðP
P0 ÞðK
PÞ,

(4)

where m is a coefficient, K is the prey carrying capacity, and P0 is the prey survival threshold (we assume 0oP0 oK), so that at low densities 0oPoP0 its growth becomes negative. The value of P0 is a characteristic of the strength of the Allee effect: the less P0 is, the less prominent is the Allee effect. Introducing, for convenience, dimensionless variables u ¼ P=K;

v ¼ Z=ðkKÞ;

t ¼ aT;

x ¼ X ða=DÞ1=2 ,

y ¼ Y ða=DÞ1=2 , where D ¼ D1 ¼ D2 and a ¼ AkK=B, from Eqs. (1)–(4) we arrive at the following equations: @u @2 u @2 u uv ¼ , þ þ guðu
bÞð1
uÞ
@t @x2 @y2 1 þ au

(5)

@v @2 v @2 v uv ¼
dv, þ þ @t @x2 @y2 1 þ au

(6)

where a ¼ K=B, b ¼ P0 =K, g ¼ K 2 mB=a and d ¼ M=a are dimensionless parameters. Note that, since the model (5)–(6) does not take into account possible existence of alternative food sources, predator cannot survive without prey and extinction of prey leads to extinction of both species. Biological invasion usually starts with a local introduction of exotic species; thus, relevant initial conditions for system (5)–(6) should be described by functions of compact support when the density of one or both species at the initial moment of time is non-zero only inside a certain domain. The shape of the domain and the profiles of the population densities can be different in different cases. In this paper, however, we are primarily concerned with the large-time dynamics of

the system when small details of the initial species distribution are not likely to play an important role (provided that the initial population size is large enough to ensure invasion success, see the beginning of the next section). For that reason, we consider the initial conditions in the form of elliptic patches:

 x
x1 2 y
y1 2 þ uðx; y; 0Þ ¼ u0 for D11 D12 p1; otherwise uðx; y; 0Þ ¼ 0,

ð7Þ



 x
x2 2 y
y2 2 þ vðx; y; 0Þ ¼ v0 for D21 D22 p1 otherwise vðx; y; 0Þ ¼ 0,

ð8Þ

where D11, D12, D21, D22, x1, y1, x2, y2, u0, v0 are parameters with obvious meaning. In the corresponding 1D case, the system (1)–(2) without Allee effect is known to exhibit a variety of travelling population waves (Kolmogorov et al., 1937; Fisher, 1937; Aronson and Weinberger, 1978; Volpert et al., 1994; Berezovskaya and Karev, 1999), some of them are linked to spatiotemporal pattern formation, cf. Sherratt et al. (1995), Petrovskii and Malchow (2000), Sherratt (2001) and Malchow and Petrovskii (2002). The impact of the Allee effect can modify the system dynamics significantly resulting in new exotic regimes such as standing chaotic patches and travelling population pulses (Morozov, 2003; Morozov et al., 2004; Petrovskii et al., 2005a). In two spatial dimensions, the system dynamics appears to be even richer. Along with regimes giving an immediate generalization of the regimes observed in 1D case, there appear more complicated patterns that cannot exist in 1D case. In the next section, the 2D dynamics of the predator–prey system with the Allee effect will be considered in much detail with a special attention to the patchy invasion and its relation to other regimes of the system dynamics.

3. Results of computer simulations Eqs. (5)–(6) with initial conditions (7)–(8) were solved numerically by finite-difference method. Computer experiments were run in a square numerical domain L  L where L ¼ 300, cf. Figs. 1–4. In most cases, we used the explicit scheme. In order to avoid numerical artifacts we checked the sensitivity of the results to the choice of the time and space steps and their values have been chosen sufficiently small. Also, some of the results were reproduced by means of using more advanced alternate directions scheme. Both numerical schemes are standard hence we do not describe them here; details and particulars can be found in Thomas (1995). Before proceeding to the simulation results, a certain remark is necessary regarding the choice of the initial

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

Fig. 1. Snapshots of the prey density (in dimensionless units) obtained at different times for parameters a ¼ 0:1, b ¼ 0:22, g ¼ 3, d ¼ 0:63. White and black/gray colors correspond to species absence and species presence at high density, respectively. Behind the expanding circular front, the species is distributed homogeneously. Predator density exhibits similar properties.

21

Fig. 2. Snapshots of the prey density obtained at different times for d ¼ 0:51, other parameters as in Fig. 1. White and black/gray colors correspond to species absence and species presence at high density, respectively. A prominent patchy species spatial distribution appears in the wake of the expanding front. Predator density exhibits similar properties.

ARTICLE IN PRESS 22

A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

Fig. 3. Snapshots of the prey density obtained at different times for d ¼ 0:43, other parameters as in Fig. 1. White and black/gray colors correspond to species absence and species presence in high density, respectively. Note that the species is absent in the wake of the expanding front. Predator density exhibits similar properties.

conditions. It is well known that, in case population growth is affected by the strong Allee effect, not every species introduction leads to successful invasion (Lewis and Kareiva, 1993; Petrovskii et al., 2005a). There exists a minimum viable population size so that invasion success can only be guaranteed when the initially invaded area and the initial population density of the alien species are not too small. Correspondingly, since in this paper we are primarily concerned with the dynamics of species spread, in our numerical simulations parameters u0, v0 and Dij in Eqs. (7)–(8) have always been chosen sufficiently large. It must be mentioned that, since the system under study depends on a relatively large number of parameters (eight for the original system (1)–(4) and four for its dimensionless version (5)–(6), its detailed numerical investigation in the whole parameter space is virtually impossible. Instead, we choose one controlling parameter and consider the changes in the invasion scenario subject to its variation. In this paper, we choose predator mortality d as the controlling factor, and keep all other parameters fixed. The choice of d as the controlling factor is justified by the results of field observations showing that intensity of predation affects the rate of prey invasion (Fagan and Bishop, 2000). Also, this choice is reasonable from the point of prospective invasive species management because, in a real ecosystem, the rate of predator mortality can be relatively easy controlled by means of additional harvesting. Results of our computer simulations show that, in the predator–prey system with the Allee effect for prey, there are several qualitatively different patterns of invasion, see also, Petrovskii et al. (2002a,b), Morozov (2003), Morozov et al. (2004) and Petrovskii et al. (2005a). Typically, for large values of d the system exhibits propagation of population fronts with stationary spatially homogeneous species distribution behind, see Fig. 1 showing snapshots of the prey density (the predator density has similar features) obtained for parameters a ¼ 0:1, b ¼ 0:22, g ¼ 3, d ¼ 0:63 and the initial conditions (7)–(8) with D211 ¼ 12:5, D212 ¼ 12:5, D221 ¼ 5, D222 ¼ 10, x1 ¼ 153:5, y1 ¼ 145, x2 ¼ 150, y2 ¼ 150, u0 ¼ 1 and v0 ¼ 0:2. For somewhat smaller d, this regime changes to travelling fronts with spatiotemporal oscillations in the wake, cf. Fig. 2 obtained for a ¼ 0:1, b ¼ 0:22, g ¼ 3, d ¼ 0:51 and the same initial conditions. (Note that these two patterns are also typical for a predator–prey system without Allee effect, see Sherratt et al. (1995) and Petrovskii and Malchow (2000).) In both of these cases, the populations are absent in front of the front and they are present in considerable densities behind the front. Apparently, both patterns correspond to successful invasion. A further decrease in d changes it to travelling population pulses (in 1D case), or rings (in 2D case) travelling/expanding from the place of the species

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

23

Fig. 4. Snapshots of the prey densities obtained at different times for the regime of patchy invasion, d ¼ 0:42, other parameters as in Fig. 1. White and black/gray colors correspond to species absence and species presence in high density, respectively. The species spreads without forming a continuous travelling front. Predator density exhibits similar properties.

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

Z

0 L

Z

0

0

4

10

1

5

2 3

0

(a)

0 x 10

100

200

300

400

500

600

t

3

16

12

2 1 8

3 4

0

(b)

0

200

400

600

800

1000

t

Fig. 5. Temporal variation of the prey population size for different invasion scenarios: (a) regimes of species invasion through propagation of continuous travelling fronts, curve 1 for d ¼ 0:51, curve 2 for d ¼ 0:46, curve 3 for d ¼ 0:43; (b) patchy invasion regime, curve 1 for d ¼ 0:42, curve 2 for d ¼ 0:417, curve 3 for d ¼ 0:414, respectively. Other parameters are same as in Fig. 1. Note different order of magnitude in (a) and (b).

cf. Fig. 3. Fig. 5b shows U(t) and V(t) obtained for patchy invasion, curve 1 for d ¼ 0:42, curve 2 for d ¼ 0:417, curve 3 for d ¼ 0:414. Thus, not only the spatial distribution of the population density exhibits prominent spatial irregularity (cf. Fig. 4) but also the population sizes experience irregular temporal fluctuations. This irregularity will come into the focus of our study in Section 3.2. 3.1. Calculation of the invasion speed

L

vðx; y; tÞ dy dx.

V ðtÞ ¼

x 10

15

Population Size of Prey

introduction, see Fig. 3 obtained for a ¼ 0:1, b ¼ 0:22, g ¼ 3, d ¼ 0:43 and the same initial conditions. In this case, the species are absent both in front of the front and behind the front. This type of the system dynamics is somewhat paradoxical from ecological point of view because successful species spread ( ¼ propagation of population front of invasive species) yet leads to invasion failure ( ¼ no species behind the front). In all above cases, the species spread takes place via propagation of a continuous travelling front. However, the impact of the strong Allee effect in the system (1)–(2) can yield a completely different pattern of population spread, in particular, the regime of patchy invasion (Petrovskii et al., 2002a,b). An example of patchy invasion is shown in Fig. 4 (obtained for parameters a ¼ 0:1, b ¼ 0:22, g ¼ 3, d ¼ 0:42 and the initial conditions the same as in Figs. 1–3). At the early stage of the system dynamics, the species spread is characterized by formation of a continuous round front (see Fig. 4a) similar to what is observed for other regimes, cf. Fig. 2a. At later time, however, the front breaks to pieces, cf. Figs. 4b and 4c, and at further stages of invasion the population spreads over space via irregular motion of separate patches. No continuous front arises again, cf. Figs. 4d–4f. The patches move, merge, disappear or produce new patches, etc. After the population patches invade over the whole domain, the spatiotemporal dynamics of the system does not change and the spatial distribution of species at any time is qualitatively similar to the one shown in Fig. 4f. We want to emphasize that, although the parameter range where the patchy invasion can be observed is relatively narrow compare to the range where invasion takes place through propagation of continuous travelling fronts, the parameter set chosen for Fig. 4 is not at all unique and the patchy invasion can be observed for other parameter values as well (cf. Petrovskii et al., 2002a; Petrovskii and Venturino, 2004; Petrovskii et al., 2005c). The invasion of introduced species is naturally followed by a gradual increase in the population sizes U and V: Z LZ L UðtÞ ¼ uðx; y; tÞ dy dx,

Population Size of Prey

24

ð9Þ

0

Fig. 5 shows U and V calculated at different time for parameter values corresponding to different invasion regimes. Fig. 5a performs U(t) and V(t) obtained for invasion through propagation of continuous travelling fronts, curve 1 for d ¼ 0:51, curve 2 for d ¼ 0:46, curve 3 for d ¼ 0:43. Here curve 1 corresponds to travelling fronts with oscillations in the wake, cf. Fig. 2, and curves 2 and 3 correspond to pulse/rings propagation,

When managing invasion of exotic species, one of the most important issues is to estimate the rate of their spread and to identify the main factors the rate can be affected by. For single-species models of population dynamics, this problem is well studied, cf. Mikhailov (1990) and Volpert et al. (1994). However, the problem becomes much more difficult when the spread of invading species is complicated by interspecies interactions, e.g. through predator–prey relations. An analytical solution of this problem is hardly possible (but see

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

where U(t) is the population size of prey given by Eq. (9). Definition (10) is thus similar to the definition of the center of mass in mechanics. The radius R of the domain can be now defined as the maximum distance between the center (xC,yC) and all the points in space where the prey density is higher then a certain prescribed value u*: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 (11a) RðtÞ ¼ maxOðtÞ ðx
xC ðtÞÞ2 þ y
yC ðtÞ , where OðtÞ ¼ fðx; yÞjuðx; y; tÞ4u g and the value of u* should be reasonably small. According to Eqs. (9)–(10), in general, the position of the center can change in the course of the system dynamics. On the other hand, it is intuitively clear that the center of the invaded area should be related to the place of original species introduction, i.e. the initial

conditions. Applying definition (9)–(10) to the initial condition (7), we immediately obtain that xC ð0Þ ¼ x1 , yC ð0Þ ¼ y1 . Thus, an alternative definition of the area radius can be as follows: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2ffi 2 RðtÞ ¼ maxOðtÞ ðx
x1 Þ þ y
y1 . (11b) We denote by R1 ðtÞ the radius of the domain given by Eq. (11b) and by R2 ðtÞ the radius given by Eq. (11a). Fig. 6a shows R1 ðtÞ (solid) and R2 ðtÞ (dashed) calculated for u ¼ 0:025 for the three different regimes shown in Figs. 2–4, i.e. for travelling fronts with irregular oscillations in the wake (I), travelling pulses (II) and patchy spread (III). Note that, contrary to intuitive expectations, definition (11b) of the domain radius (cf. the solid curve) appears to be more appropriate for speed estimation because it produces a curve very close to a straight line, at least, for sufficiently

150

II

I

III

R1(t), R2(t)

100

50

0

(a)

0

50

100

150

200

250

t

300

350

400

450

500

550

0.7 0.6

II

III

I

0.5

C1 , C2

Petrovskii et al., 2005b) and a thorough numerical investigation is very difficult as well because of high dimensionality of the corresponding parameter space. Let us note, however, that in a real ecosystem the rate of spread depends on numerous factors that unlikely can be all taken into account by a single mathematical model. From this perspective, it seems to be more important to understand the structural dependence of the invasion speed on the type of spread (e.g. patchy or not patchy) rather than to find out its dependence on all model parameters. For the invasion regimes showing continuous travelling fronts, estimation of the invasion speed does not bring much difficulty. In an idealized case of purely homogeneous environment when the shape of invaded area is very close to a circle, it can be readily obtained from the rate of increase of the radius R of the area. For a more realistic case of inhomogeneous environment when the shape of the invaded area can be rather distorted, a special averaging procedure was developed which is still essentially based on the existence of the continuous travelling front separating invaded and noninvaded areas, see Shigesada and Kawasaki (1997) for details and further references. In case of patchy invasion, however, these approaches do not immediately apply because it is hardly possible to trace an exact border between invaded and non-invaded areas. In order to overcome this difficulty, we have developed two different methods to estimate the invasion speed that are not based on existence of continuous travelling front. The idea of our approach is as follows. First, we define the center (xC,yC) of the invaded area by means of the following relations: Z LZ L 1 xC ðtÞ ¼ xuðx; y; tÞ dy dx, UðtÞ 0 0 Z LZ L 1 yuðx; y; tÞ dy dx, ð10Þ yC ðtÞ ¼ UðtÞ 0 0

25

0.4 0.3 0.2 0.1 0

(b)

0.42

0.44

0.46

0.48

0.50

δ

Fig. 6. (a) The radii R1 (solid) and R2 (dashed) of invaded area vs time obtained for three different types of invasion regimes (see Figs. 2–4) by using two different methods, cf. Eqs. (11b) and (11a), respectively. Curves I are obtained for parameters of Fig. 2 (patterns on the wake), curve II for parameters of Fig. 3 (travelling pulses), curve III for parameters of Fig. 4 (patchy invasion). Other parameters are the same as in Fig. 1; (b) speed of invasion calculated for different d based on Eq. (12) with R ¼ R1 (dots) and R ¼ R2 (asterisks). Parameter domains I and II correspond to propagation of continuous travelling fronts with irregular spatiotemporal patterns and species extinction in the wake, respectively; domain III corresponds to patchy invasion.

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

26

large t. Oscillations in the dashed curve indicate that the position of the dynamical center of the invaded domain given by Eqs. (9)–(10) fluctuates with time, likely due to local oscillatory population dynamics. The average slopes of the curves give an estimate of the invasion speed. Let us note, however, that, generally speaking, in a 2D case with axial symmetry the rate of invasion depends on the radius of the growing domain, i.e. C ¼ CðRÞ, and the stationary value can only be reached in the large-time limit when R tends to infinity (Mikhailov, 1990). Indeed, it is easily seen that the process has different stages, especially, in the regimes II and III. Thus, to estimate the invasion speed we use the following equation: 00

0

¯ ¼R
R , C t00
t0

(12)

where R0 should be chosen sufficiently large. Accounting for the fact that the radius of the domain grow non 0 can  0 0 0 0 monotonously, we define t ¼ min t jt : Rðt i i i ¼ÞR ;   t00 ¼ min t00i jt00i : Rðt00i Þ ¼ R00 ; i.e., t0 and t00 are the moments when the radius of the growing domain for the first time reaches prescribed values R0 and R00 , respectively. To estimate the speed of invasion by means of Eq. (12), we chose R0 ¼ 100 and R00 ¼ 150; for R(t)o100, the impact of transients and dependence of the speed on R is still strong and for R00 4150 the species nearly reaches the domain border. In order to reveal the ecological implications of the patchy invasion, we compare the invasion speed for different regimes using predator mortality d as a controlling parameter, cf. the beginning of Section 3. Fig. 6b shows the speed C¯ 1 ðdÞ (dots) and C¯ 2 ðdÞ (asterisks) calculated for R1 ðtÞ, R2 ðtÞ, respectively (other parameters are the same as in Figs. 1–4). The symbols I, II and III, show which regime of invasion takes place for given range of d. For do0:412, extinction of both species takes place for any initial conditions; for d40:61, regime I (travelling fronts with oscillations in the wake) changes to propagation of fronts with homogeneous species distribution in the wake. Thus, the two methods to estimate the invasion speed ¯1 C ¯ 2  C, ¯ where C ¯ is the true give very close values, C ¯ value of the speed. Our numerical data show that C exhibits a clear tendency to decrease with a decrease in predator mortality d. This result is in a good agreement with data of field observations (Fagan and Bishop, 2000) showing that the speed of invasion is the less the higher is the predation. Fig. 6b shows that the transition between the patterns of species spread with continuous travelling fronts, i.e. types I and II, does not lead to any significant change in ¯ On the contrary, the transition from domain II to C. domain III, i.e. from a pattern of spread via propagation of continuous population front to the patchy invasion, is

characterized by a remarkable drop in the invasion speed. This dramatic decrease in C¯ can be interpreted heuristically as a transition to a different mechanism of species spread which will be discussed in Section 4. It should be mentioned that, for d 2 ½0:412; 0:427, the choice of initial conditions is important. Our numerical simulations show that the same parameters in Eqs. (7)–(8) that lead to patchy invasion for a certain d from this range can lead to species extinction for a slightly different value d0 . However, for this new value d0 another parameter set in Eqs. (7)–(8) can be found that lead again to the patchy spread. This singular dependence of the system dynamics on the initial conditions brings a certain difficulty for a regular study of the patchy invasion, e.g. when investigating invasion speed dependence on given controlling parameter. One possible solution can be reached by tuning the initial conditions for each new d. Our simulations show that parameter changes in Eqs. (7)–(8) do not affect much the value of invasion speed as far as the type of spread is retained. Another approach to the choice of appropriate initial conditions, which we actually used in our study, is based on the observation that, in case the patchy spread has been initiated, small variation of d normally does not destroy it. Thus, it is convenient to use an embryo of the patchy distribution obtained for a particular d at a particular moment of time as the initial condition in the simulations for other values of d. Specifically, for that purpose we used the pattern shown in Fig. 4c. 3.2. Chaotic properties of the patchy dynamics Irregular temporal fluctuations of the population sizes shown in Fig. 5b seem to indicate that the system dynamics is chaotic. It should be noted that the corresponding homogeneous system, i.e. Eqs. (1)–(2) without diffusion terms, cannot exhibit behavior more complex than periodic. Thus, chaos in the system (1)–(2) is an essential result of the system dynamics in space and should be qualified as spatiotemporal (Pascual, 1993; Sherratt et al., 1995; Petrovskii and Malchow, 1999, 2001; Sherratt, 2001; Petrovskii et al., 2003). An apparent irregularity of the system behavior, however, does not necessarily correspond to chaotic dynamics and a more careful analysis is needed. That can be done in different ways. For instance, one of the basic properties of deterministic chaos is its sensitivity to the initial conditions so that the distance between the solutions corresponding to perturbed and unperturbed initial conditions grows with time (Nayfeh and Balachandran, 1995). In order to check whether sensitivity to the initial conditions takes place for the regime of patchy invasion, we varied parameters in Eqs. (7)–(8). Fig. 7 shows the prey density vs. time for the unperturbed system (solid line) and the perturbed system (dashed line) obtained at

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

The procedure is as follows. At certain moment t0 4tinv , the following small perturbation is applied to the system:


2pðx þ yÞ u1 ðx; y; t0 Þ ¼ uðx; y; t0 Þ 1 þ e cos L0

 2pðx
yÞ
e cos ð15Þ L0

1

Prey Density

0.8

0.6

0.4

0.2

0

27

0

200

400

t

600

800

1000

Fig. 7. Temporal oscillations of the prey density at a fixed point, x ¼ 155, y ¼ 145 subject to a small perturbation of the initial conditions, solid curve for the undisturbed system, dashed curve for the disturbed system.

a fixed point in space, x ¼ 155, y ¼ 145, for the same parameter values as in Fig. 4. The perturbed system was obtained from the unperturbed one by multiplying u0 by factor 1.00001. It is readily seen that up to t ¼ 550 the difference between the two solutions is very small; however, for t4550 the difference grows rapidly and the discrepancy soon becomes on the order of the solutions themselves. Qualitatively similar results were also obtained for perturbation of other parameters and for other positions in space. A more rigorous definition of deterministic chaos, however, requires the discrepancy between the undisturbed and disturbed systems not just to be increasing but to be increasing exponentially. There is a number of methods to calculate the dominant Lyapunov exponent lmax (cf. Wolf et al., 1985; Kantz and Schreiber, 1997). In this paper, we have estimated the value of lmax based on its definition. In order to take into consideration both dynamical variables, i.e. the prey and predator densities, and also to take into account the spatial aspect, we analyse the behavior of the following two values:   hu ðtÞ ¼ uðx; y; tÞ
u1 ðx; y; tÞ,   hv ðtÞ ¼ vðx; y; tÞ
v1 ðx; y; tÞ, ð13Þ

v1 ðx; y; t0 Þ ¼ vðx; y; t0 Þ.

(16)

The dynamical variables of perturbed and unperturbed systems are then used to construct the discrepancies hu and hv . Fig. 8 shows lnðhu Þ and lnðhv Þ vs. time calculated for e ¼ 10
9 , t0 ¼ 1250, L0 ¼ 100 and for the same parameters as in Fig. 4. The straight lines are obtained by the least square approximation. It is easily seen that the distances hu and hv between initially close solutions u,v and u1 ,v1 grow exponentially with time (up to small fluctuations). The slope of the straight lines (which appears to be the same for both lines) gives an estimate for the dominant Lyapunov exponent, lmax ¼ 0:035  0:001. Another evidence of chaotic nature of the patchy dynamics can be obtained from the behavior of relevant power spectra (Nayfeh and Balachandran, 1995). We apply this method to the time series of the population density obtained in a fixed point x ¼ 160, y ¼ 155, see Fig. 9a (parameters are the same as in Fig. 4). As well as the population sizes, the local dynamics exhibits apparent irregularity giving an indication of deterministic chaos. This irregularity perhaps can be seen even better in the phase plane of the local population densities, cf. Fig. 9b obtained for a somewhat longer period of time compare to Fig. 9a. Indeed, the power spectrum of the prey density time series, see Fig. 9c, has properties typical for chaotic dynamics. In particular, it shows the rate of decay slower than exponential (the

2 0

where

-2

Z

L

w2 ðx; yÞ dx dy

0

1=2

-4

.

(14)

0

ln (hv)

L

ln (hu ),


Z kwk ¼

-6 -8 -10

Strictly speaking, the analysis of the system sensitivity to phase perturbations only can be applied to stationary (in the statistical sense) time series. Meanwhile, until the species invade the whole domain, the time series are apparently transient, cf. Fig. 5. Thus, to estimate the dominant Lyapunov exponent, the system dynamics was studied for the post-invasion stage, i.e. for t4tinv where tinv is the time that takes the species to spread over the whole area. For the parameters of Fig. 4, tinv 1000.

Straight line approximation: Prey Predator

-12 -14 -16 1250

1350

1450

1550

1650

1750

t

Fig. 8. Discrepancies hu ðtÞ and hv ðtÞ between two initially closed system trajectories, cf. Eqs. (13)–(16), in the regime of patchy invasion. Parameters are the same as in Fig. 4. The slope of the straight line gives the value of the dominant Lyapunov exponent, lmax ¼ 0:035  0:001.

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

28

correlation between different patches in space, we calculate the spatial cross-correlation function. Correlation between temporal variations of prey species density at two different sites in space as a function of the intersite distance is given, in a rather general case is given by (Nisbet and Gurney, 1982; Nayfeh and Balachandran, 1995)

Prey and Predator density

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 3000

3500

4000

(a)

4500

5000

t 0.8

1 1 Kðx; y; x0 ; y0 Þ ¼ lim T!1 T su su0 Z t0 þT ðuðx0 ; y0 ; tÞ
u¯ 0 Þðuðx; y; tÞ
u¯ Þ dt, 

Predator density

t0 0.6

ð17Þ where (x,y) and ðx0 ; y0 Þ give the position of the sites, u¯ and u¯ 0 are the corresponding mean prey densities, su and su0 are the standard deviations. In Eq. (17), formally, the value of K depends on the position of the sites. Assuming, based on our numerical simulations, that for sufficiently large time the system becomes statistically homogeneous, the cross-correlation function Kðx0 ; y0 ; x; yÞ is a function only of the distance s between the two points. Then, from Eq. (17) we obtain: Z 1 1 t0 þT ðuðx0 ; y0 ; tÞ
u¯ 0 Þðuðx; y; tÞ
u¯ Þ dt, KðsÞ ¼ lim  2 T!1 T s t0

0.4

0.2

0.2

0

(b)

0.4

0.6

0.8

Prey density

2 0

Power

-2 -4 -6

(18)

-10

(c)

0

0.1

0.2

0.3

0.4

1/Time

Fig. 9. Dynamics of the population densities calculated in a fixed point, x ¼ 160, y ¼ 155 for the regime of patchy invasion (parameters are the same as in Fig. 4): (a) density of prey (solid) and predator (dashed) vs. time and (b) local phase plane of the system and (c) power spectrum of the oscillating prey density (semilogarithmic).

exponential decay would correspond to a straight line envelope). The power spectrum of predator density possesses similar properties. We want to emphasize that the results presented in Fig. 9 are not point-specific; our numerical results show that the same type of behavior is observed for any other position in space. The results of the above analysis demonstrate that the dynamics of the system (5)–(6) in the regime of patchy spread (cf. Fig. 4) is chaotic. It should be mentioned that, although solid evidence is obtained only for the post-invasion stage when the time series become stationary, the sensitivity to the initial conditions that clearly manifest itself already during the invasion stage, see Fig. 7, makes it possible to conclude that the dynamics is chaotic at earlier stages as well. The above results concern the temporal variations of the species density. Meanwhile, the spatial aspect is very important as well, cf. the beginning of this section. In order to address this issue and to reveal possible

where su0 ¼ su ¼ s, u¯ 0 ¼ u¯ and s¼ ffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðx
x0 Þ2 þ ð y
y0 Þ2 : Fig. 10 shows function K(s) calculated according (18) with t0 ¼ 2000 and T ¼ 10000 along the straight line x ¼ 150, y ¼ 25 þ s. Parameters are the same as in Fig. 4. The behavior of the crosscorrelation function is typical for chaotic dynamics

1.2 1 0.8

K (s)

-8

0.6 0.4 0.2 0 -0.2 0

40

80

120

s

160

200

240

Fig. 10. Cross-correlation function K(s), see Eq. (18), calculated for the regime of patchy invasion. Parameters are the same as in Fig. 4.

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

(Nayfeh and Balachandran, 1995). The first zero of K(s) gives the correlation length of the system: starting from that distance the local temporal fluctuations can be regarded as independent. For parameters of Fig. 4, Lcorr  18: Thus, the whole area L  L appears to be dynamically split into NðL=Lcorr Þ2 sub-domains with virtually independent temporal behavior. Note that Eqs. (17)–(18) can be written in terms of predator density as well. For the sake of brevity we do not present corresponding results here because they do not lead to any new insights: the properties of the predator-based cross-correlation function are very similar to those shown in Fig. 10. Since the dynamics of the system (5)–(6) is essentially spatiotemporal, it also makes sense to look at its properties in terms of the global dynamical variables that would explicitly take into account the spatial dimensions. The simplest variables of that kind seem to be the population sizes U and V or the spatially 0.3

u t, v

t

0.25 0.2 0.15 0.1 0.05 0 3000

(a)

3500

4000

4500

29

averaged densities, huit ¼ L
2 U; hvit ¼ L
2 V : Time series of these values are shown in Fig. 11a. As it could be expected, the dynamical splitting of the system to a number of quasi-independent oscillators results in a much smaller amplitude of the oscillations of the spatially averaged values than it was for the local densities, cf. Figs. 9a and 11a. (A similar type of spacedependence was observed earlier for another system, see Bascompte and Sole` (1994).) As well as the local variables, the oscillations of the spatially averaged densities exhibit prominent irregularity, see Fig. 11b. The properties of the corresponding power spectra, see Fig. 11c, indicates that this irregularity arises as a result of chaotic dynamics. In conclusion of this section, we want to mention that, although our attention has been mostly focused on the patchy invasion, it is not the only regime in the system (5)–(6) that exhibit chaotic dynamics. In particular, analysis of our numerical data indicates that irregular oscillations arising behind the continuous travelling front, cf. Fig. 2, also correspond to spatiotemporal chaos. That regime of the system dynamics, however, has already been observed and studied in much detail for a predator–prey system without Allee effect (Pascual, 1993; Sherratt et al., 1995; Petrovskii and Malchow, 1999, 2000, 2001; Petrovskii et al., 2003), and thus will not be considered in this paper.

5000

t 0.16

4. Discussion

0.15 0.14

0.12

v

t

0.13

0.11 0.10 0.09 0.08 0.07

0.10

0.12

0.14

0.16

u

(b)

0.18

0.20

0.25

0.3

0.22

t

2 0

Power

-2 -4 -6 -8 -10 -12

(c)

0

0.05

0.1

0.15

0.2

0.35

0.4

1/Time

Fig. 11. Dynamics of the spatially average species densities ou4t and ov4t for the regime of patchy invasion (parameters are the same as in Fig. 4): (a) spatially average density of prey ou4t (solid) and predator ov4t (dashed) vs. time and (b) phase plane of the average species densities and (c) The power spectrum of the oscillating spatially average prey density (semilogarithmic).

In this paper, we have considered regimes of biological invasion in a 2D predator–prey system where the prey growth is damped by the Allee effect. It was shown in our earlier work (Petrovskii et al., 2002a,b) that, as a result of the strong Allee effect, along with already well-known regimes of spread such as invasion of species via propagation of travelling population front with either spatially homogeneous or patchy species distribution in the wake, this system exhibits a completely new invasion scenario when the species spread is not preceded by propagation of any continuous front. Instead, the populations spread over space via irregular motion of separate patches of high population density, the population density between patches being nearly zero, cf. Fig. 4. We want to emphasize that this patchiness is not induced by any environmental heterogeneity and is self-organized. Note that the patterns shown in Fig. 4 were observed in the case of equal species diffusivity. Thus, the pattern formation is not the result of classical Turing mechanism that requires sufficiently different diffusion coefficients (Segel and Jackson, 1972) and has a different origin. We also want to mention that our model is fully deterministic and thus the patchy spread cannot be

ARTICLE IN PRESS 30

A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

attributed to the impact of stochastic factors (cf. Lewis, 2000; Lewis and Pacala, 2000). The above-mentioned papers by Petrovskii et al. (2002a,b) were more concerned with reporting the new phenomenon rather than with its comprehensive investigation. In this paper, the phenomenon of deterministic patchy invasion was considered in much detail. In particular, two important issues were addressed: (i) how the type of spread, i.e. with or without continuous travelling front, affects the speed of invasion and (ii) what is the degree of corresponding spatiotemporal complexity of the system dynamics. Concerning the first issue, we showed by means of extensive numerical experiments that, choosing the predator mortality d as a controlling parameter, the regime of patchy invasion corresponds to a remarkable drop in the invasion speed, cf. Fig. 6b. A sufficiently large decrease in d brings the invasion speed to zero. This fact is in a very good agreement with other recent results (Petrovskii and Venturino, 2004; Petrovskii et al., 2005c) showing that the patchy spread gives the scenario of species invasion at the edge of extinction: a reasonably small variation of the system parameters either restore the usual invasion regime via propagation of a continuous travelling front or bring the species to extinction. Moreover, results of our computer experiments, see Petrovskii and Venturino (2004) and Petrovskii et al. (2005c) for details, show that the dimensionality of space plays a crucial role: for the parameters when the system exhibits the patchy spread in two spatial dimensions, the species go extinct in the corresponding 1D system regardless the choice of the initial conditions. A heuristic description of invasion failure in 1D case is as follows. At early stages of species spread, the travelling front of the prey density propagates into empty space being followed by a travelling front of predator density, the dynamics being similar to the socalled predator–prey pursuit scenario (Murray, 1989). The speed of prey wave appears to be lower than that of the predator wave; as a result, prey is caught up by predator. After some oscillations, predator decreases the prey density below the survival threshold b at every location in space, and extinction of both species takes place. In the 2D system, however, the patch border is curvilinear; thus, prey can escape through the lateral sides and create separate patches. Each new patch then starts growing through formation and expansion of circular population front (cf. Fig. 4a) until prey is caught by predator, and this scenario occurs again and again resulting in a prominent patchy structure. The above description also helps to understand why the rate of species invasion in the regime of patchy spread appears to be smaller than it is in the regimes of continuous front propagation. First, the travelling front

of prey is periodically stopped by predator. Second, in every case formation of new patches takes a certain time. Finally, the growth of the newly formed patches normally leads to species spread not in the radial direction but at different angles. Note that the effect of patch border curvature is positive for prey and negative for predator. Indeed, the larger is curvature (i.e. the smaller is the patch radius) the higher is the rate of population density decrease due to its out-flux through the patch border. A decrease in the population densities inside given patch diminishes the impact of predation on the prey growth and also diminishes the growth rate of predator, cf. Eqs. (5)–(6) where the term describing interspecies interaction has different sign in the equation for prey and in the equation for predator. Thus, prey has more chances to survive in small patches where the border curvature is large than in large patches where the border shape is close to a straight line. Regarding the complexity of the system dynamics in the regime of patchy invasion, we showed that it is chaotic both for local densities and for spatially averaged values, cf. Figs. 9 and 11. Results of our numerical experiments (omitting details for the sake of brevity) indicate that the correlation dimension of the strange attractor is relatively high, Dcorr 46. For comparison, in a similar predator-prey system but with the logistic prey growth the correlation dimension Dcorr was found to be between 3 and 4 (Pascual, 1993). Thus, it can be inferred that the impact of the Allee effect significantly increases system complexity. 4.1. Statistical properties of the patchy dynamics As it was mentioned above, spatiotemporal chaos in the predator–prey system with the Allee effect has been observed for some other regimes as well (Morozov et al., 2004; Petrovskii et al., 2005a), e.g. in the case of pattern formation in the wake of a travelling population front, see Fig. 2. (Note that the patterns in the wake looks qualitatively similar to the patterns arising as a result of patchy invasion, cf. the middle of Figs. 2c and 4c.) Thus, chaoticity of the system dynamics alone cannot be used to distinguish between the two regimes. The question that remains open is whether the absence of a continuous travelling front is the only difference, albeit important, between the two regimes of irregular spatiotemporal pattern formation, cf. Figs. 2 and 4. In order to address this issue, we have to make a deeper insight into the properties of the corresponding dynamics. High complexity of the system behavior (which is quantified by the high correlation dimension of the attractor) indicates that it is, in some sense, rather close to stochastic dynamics. Indeed, it can be shown that, on the time-scale larger than the correlation time tcorr of the system, its dynamics can be reproduced by

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35 35

Distribution function

30 25 20 15 10 5 0 0.1

0.12

0.14

0.18

0.16

0.2

0.22

Average Prey Density

(a) 16 14

Distribution function

means of purely stochastic process (A.Y.Morozov, unpublished manuscript). These observations seem to make possible to describe the system dynamics in terms of probability of its different states. In order to apply this approach, we first estimate the correlation time of the system by means of calculating the autocorrelation functions and finding their first zero. Then we generate the time series of local population densities (calculated in a fixed point x ¼ 160, y ¼ 155) and of their spatially averaged values, the time lag between any two consequent terms in these series being chosen equal to tcorr . Assuming that any two measurements in these series are independent, we then restore the probability distribution functions (PDF) of system states. Since the time series obtained for prey and predator exhibit similar properties, below we show only the results obtained for prey dynamics. Fig. 12a shows the PDF of the average prey density obtained for the patchy invasion regime for the same parameter set as in Figs. 4 and 9 (for these parameters, tcorr ¼ 35). One can easily see that the shape of the PDF is very close to the Gaussian law. Indeed, the envelope shown by the thick solid line gives the normal distribution where the values of the mean mu and the standard deviation su are obtained from the same time series; here mu ¼ 0:168 and su ¼ 0:012. The result that the dynamics of the spatially averaged values is very well described by the normal distribution is easily understood if we recall that, in the patchy regime, the system is virtually reduced to an ensemble of quasi-independent oscillators, cf. the lines below Eq. (18). Moreover, it can be expected that the PDFs of spatially averaged values constructed for the irregular patterns preceded by continuous front propagation, cf. Fig. 2, will be very close to the normal distribution as well because, as it was mentioned above, the system in this regime also exhibits spatiotemporal chaos. Fig. 13a shows the average prey density PDF obtained for the parameters of Fig. 2, thick solid curve giving the normal distribution. Other useful information can be obtained from the PDFs of the local densities. Figs. 12b and 13b shows the PDFs of the prey density calculated for the patchy spread (Fig. 4) and for patterns in the wake of the population front (Fig. 2), respectively. Note that, in both cases, the value of the survival threshold is the same, b ¼ 0:22. It is readily seen that, in the case of pattern formation preceded by front propagation, high values of the prey density are more probable than low values. On the contrary, in the case of the patchy spread, low values of prey density are more probable than high ones. A brief inspection of Fig. 12b immediately leads to a somewhat counter-intuitive result that, in the regime of patchy spread, the probability to detect the prey density below the survival threshold is higher than the

31

12 10 8 6 4 2 0 0

(b)

0.2

0.4

0.6

0.8

1

Local Prey Density

Fig. 12. Probability distribution functions of the system states calculated for the patchy invasion, parameters are the same as in Fig. 4: (a) probability distribution function of the spatially average prey density ou4t. Solid line shows normal distribution with the same mean and standard deviation and (b) Probability distribution function of the local prey density at a fixed point ðx ¼ 160; y ¼ 155Þ.

probability to find it above the threshold. Moreover, as it is immediately inferred from Figs. 12a and 13a, the spatially averaged prey density appears to be well above and well below the survival threshold for the regimes of invasion with and without continuous population front, respectively. 4.2. Impact of environmental heterogeneity All the above results have been obtained under an idealized assumption of environmental homogeneity. In reality, however, species spread often takes place in heterogeneous environment or even fragmented landscape. An important question thus arises whether the self-organized patchy invasion can be observed when the model parameters become space-dependent. Apparently, it is a rather diversified issue because the actual outcome of the interplay between intrinsic population dynamics and external forcing depends on the nature, geometrical shape and the magnitude of environmental fluctuations, cf. Sherratt et al. (2003). Hence, in this paper we make only an early insight into the matter in

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

32 35

Distribution function

30 25 20 15 10 5 0 0.38

0.4

(a)

0.42

0.44

0.46

0.48

Average Prey Density 2.5

Distribution function

2.0

(b)

4.3. Concluding remarks

1.5

1.0

0.5

0

rather severe conditions for prey, cf. the white spots where do0:41 (which would correspond to species extinction in spatially homogeneous case), invasion appears to be successful and species spread takes place via essentially the same scenario of formation and motion of separated patches. Qualitatively, the dynamics looks very similar to the patchy spread in homogeneous environment, cf. Fig. 4, typical size and shape of the spatial structures being close to those obtained earlier. Note that the characteristic width of emerging patches is several times smaller than the period of spatial variations in d, which indicates that the patchy spread is still controlled by the intrinsic mechanism rather than by environmental heterogeneity.

0

0.2

0.4

0.6

0.8

1

Local Prey Density

Fig. 13. Probability distribution functions of the system states calculated for the regime of invasion via propagation of continuous population front with spatiotemporal patterns in the wake, parameters are the same as in Fig. 2: (a) probability distribution function of the spatially average prey density ou4t. Solid line shows normal distribution with the same mean and standard deviation and (b) Probability distribution function of the local prey density at a fixed point ðx ¼ 160; y ¼ 155Þ.

order to demonstrate robustness of the patchy invasion; a detailed study of this problem will be done elsewhere. In order to keep it consistent with the earlier analysis, we assume that it is the predator mortality d that depends on position in space while other parameters remain constant. Specifically, we consider heterogeneity of the following type:



2pðx þ yÞ 2pðx
yÞ dðx; yÞ ¼ d0 þ D sin
D sin , L1 L1 (19) where d0 and D are parameters with the meaning of the average mortality and the amplitude of spatial variations, respectively, and L1 quantifies the scale of the variations. Function dðx; yÞ is shown in Fig. 14a for d0 ¼ 0:42, D ¼ 0:02, L1 ¼ 125:67. We consider species invasion in heterogeneous environment (19) for the same value of parameters a, b, g and the same initial condition as in Figs. 1–4. Snapshots of the prey density at different moments of time are shown in Figs. 14b–14e. In spite of the presence of areas with

Apart from environmental homogeneity, another assumption that has been used throughout this paper is that diffusivity is the same for prey and for predator. It should be emphasized, however, that existence of patchy spread is not restricted to this case. In our numerical simulations, we checked whether the main results remain qualitatively the same when D2/D1 is not equal to unity. Computer simulations show that, at least in the range 0:7oD2 =D1 o1:5, the patchy spread with the properties described above can be observed as well, although parameters a, b, g, d should be chosen slightly different. An inspection of the spatial patterns observed during the patchy invasion shows that there is an intrinsic spatial scale, which corresponds to the typical size of the patches, cf. Figs. 4 and 10. For the given parameters, it appears to be, in dimensionless units, somewhere between 10 and 20. A question of interest is to what spatial scale and/or to what invasive species it may correspond in real population communities. According to the definition of dimensionless variables, see the lines above Eqs. (5)–(6), the relation between the dimensional R and dimensionless r spatial scales is given as R ¼ rðD=AkaÞ1=2 ; where a ¼ K=B. Note that, in our simulations, parameter a is always fixed at a hypothetical value of 0.1, which means that the effect of saturation in predator response is insignificant. As for D and Ak (here we recall that Ak has the meaning of predator maximum growth rate), they can be different for different species. As an example, we consider the vole–weasels interaction, cf. Sherratt et al. (2002), with D ¼ 0:2 km2 year
1 and Ak ¼ 2:7 year
1 as typical values. We then obtain that R lies between 8.5 and 17 km, which seems to be ecologically reasonable. Note that an increase/decrease in D or Ak as much as 10 times corresponds to only about three times increase/decrease in R; thus, an estimate for R to be between a few kms and a few dozens kms is likely to remain valid for many other terrestrial species.

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

33

Fig. 14. Patchy invasion in heterogeneous environment: (a) spatial dependence of d as given by Eq. (19) for parameters d0 ¼ 0:42, D ¼ 0:02, L1 ¼ 125:67 and (b)–(e) snapshots of the prey density obtained at different times. White and black/gray colors correspond to prey absence and prey presence at high density, respectively. Other parameters are the same as in Fig. 4.

ARTICLE IN PRESS 34

A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35

From the point of prospective application of our results to invasive species management and control, the singular dependence of invasion success on predator mortality is an important result. High predator mortality corresponds to successful invasion (cf. Figs. 1, 2 and also 6b): predator is too weak to stop the invasion of prey. Gradual decrease in d (which corresponds to an increase in predation) finally leads to invasion failure when propagation of continuous travelling fronts with either homogeneous or patchy species distribution behind changes to propagation of travelling population pulses where population density in the wake is zero, cf. Fig. 3. However, a further decrease in d restores successful invasion in the form of the patchy spread. This result seems to have important implications for managing invasive species and the problem of biological control. Indeed, current approaches to invasive species control are often based on explicit or implicit assumption that the higher is the press on invasive species, the more probable is its invasion failure. Our results, however, indicate that, instead of this monotone dependence, there can exist an optimal magnitude of controlling efforts.

Acknowledgements This research was partially supported by the US National Science Foundation under Grants DEB0409984 and DEB-0080529, by Russian Foundation for Basic Research under Grants 03-04-48018 and 0404-49649, and by the University of California Agricultural Experiment Station. References Allee, W.C., 1938. The Social Life of Animals. Norton and Co., New York. Andow, D.A., Kareiva, P.M., Levin, S.A., Okubo, A., 1990. Spread of invading organisms. Landscape Ecol. 4, 177–188. Aronson, D.G., Weinberger, H.F., 1978. Multidimensional nonlinear diffusion arising in population genetics. Adv. Math. 30, 33–76. Bascompte, J., Sole`, R.V., 1994. Spatially induced bifurcations in single-species population dynamics. J. Anim. Ecol. 63, 256–264. Berezovskaya, F.S., Karev, G.P., 1999. Bifurcations of travelling waves in population taxis models. Phys.–Usp. 169, 1011–1024. Courchamp, F., Clutton-Brock, T., Grenfell, B., 1999. Inverse density dependence and the Allee effect. Trends Ecol. Evol. 14, 405–410. Davis, M.B., Calcote, R.R., Sugita, S., Takahara, H., 1998. Patchy invasion and the origin of a Hemlock-Hardwoods forest mosaic. Ecology 79, 2641–2659. Dennis, B., 1989. Allee effect: population growth, critical density and the chance of extinction. Nat. Resources Modeling 3, 481–538. Drake, J.A., Mooney, H.A., di Castri, F., Groves, R.H., Kruger, F.J., Rejmanek, M., Williamson, M. (Eds.), 1989. Biological Invasions: Global Perspective. Wiley, New York. Fagan, W.F., Bishop, J.G., 2000. Trophic interactions during primary succession: herbivores slow a plant reinvasion at Mount St. Helens. Am. Nat. 155, 238–251.

Fagan, W.F., Lewis, M.A., Neurbert, M.G., van der Driessche, P., 2002. Invasion theory and biological control. Ecol. Lett. 5, 148–157. Fisher, R., 1937. The wave of advance of advantageous genes. Ann. Eugenics 7, 355–369. Hengeveld, R., 1989. Dynamics of Biological Invasions. Chapman & Hall, London. Holling, C.S., 1965. The functional response of predator to prey density and its role in mimicry and population regulation. Mem. Entomol. Soc. Canada 45, 1–60. Holmes, E.E., Lewis, M.A., Banks, J.E., Veit, R.R., 1994. Partial differential equations in ecology: spatial interactions and population dynamics. Ecology 75, 17–29. Lewis, M.A., 2000. Spread rate for a nonlinear stochastic invasion. J. Math. Biol. 41, 430–454. Lewis, M.A., Kareiva, P., 1993. Allee dynamics and the spread of invading organisms. Theor. Popul. Biol. 43, 141–158. Lewis, M.A., Pacala, S., 2000. Modeling and analysis of stochastic invasion processes. J. Math. Biol. 41, 387–429. Kantz, H., Schreiber, T., 1997. Nonlinear Time Series Analysis. Cambridge University Press, Cambridge. Kawasaki, K., Mochizuki, A., Matsushita, M., Umeda, T., Shigesada, N., 1997. Modelling spatio-temporal patterns generated by Bacillus subtilis. J. Theor. Biol. 188, 177–185. Kolb, A., Alpert, P., Enters, D., Holzapfel, C., 2004. Environmental stress and plant community invasibility in a coastal grassland in California. ESA Annual Conference, Portland, Oregon, July 30–August 6, 2004. Kolmogorov, A.N., Petrovsky, I.G., Piskunov, N.S., 1937. Investigation of the equation of diffusion combined with increasing of the substance and its application to a biology problem. Bull. Moscow State Univ. Ser. A: Math. Mech. 1 (6), 1–25. Malchow, H., Petrovskii, S.V., 2002. Dynamical stabilization of an unstable equilibrium in biological and chemical systems. Math. Comput. Modelling 36, 307–319. Malchow, H., Petrovskii, S.V., Medvinsky, A.B., 2002. Numerical study of plankton-fish dynamics in a spatially structured and noisy environment. Ecol. Modelling 149, 247–255. Malchow, H., Hilker, F., Petrovskii, S.V., 2004. Noise and productivity dependence of spatiotemporal pattern formation in a preypredator system. Discrete Continuous Dynamical Systems 4, 707–713. Medvinsky, A.B., Petrovskii, S.V., Tikhonova, I.A., Malchow, H., Li, B.-L., 2002. Spatiotemporal complexity of plankton and fish dynamics. SIAM Rev 44, 311–370. Mikhailov, A.S., 1990. Foundations of Synergetics I. Springer, Berlin. Mimura, M., Sakaguchi, H., Matsushita, M., 2000. Reaction–diffusion modeling of bacterial colony patterns. Physica A 282, 283–303. Morozov, A.Y., 2003. Mathematical modelling of ecological factors affecting stability and pattern formation in marine plankton communities. Ph.D. Dissertation. Shirshov Institute of Oceanology, Moscow. Morozov, A.Y., Petrovskii, S.V., Li, B.-L., 2004. Bifurcation and chaos in a predator–prey system with the Allee effect. Proc. R. Soc. London B 271, 1407–1414. Murray, J.D., 1989. Mathematical Biology. Springer, Berlin. Nayfeh, A.H., Balachandran, B., 1995. Applied Nonlinear Dynamics. Wiley, New York. Nisbet, R.M., Gurney, W.S.C., 1982. Modeling Fluctuating Populations. Wiley, Chichester. Owen, M.R., Lewis, M.A., 2001. How predation can slow, stop or reverse a prey invasion. Bull. Math. Biol. 63, 655–684. Pascual, M., 1993. Diffusion-induced chaos in a spatial predator–prey system. Proc. R. Soc. London B 251, 1–7.

ARTICLE IN PRESS A. Morozov et al. / Journal of Theoretical Biology 238 (2006) 18–35 Petrovskii, S.V., Malchow, H., 1999. A minimal model of pattern formation in a prey–predator system. Math. Comput. Modelling 59, 157–174. Petrovskii, S.V., Malchow, H., 2000. Critical phenomena in plankton communities: KISS model revisited. Nonlinear Anal.: Real World Appl. 1, 37–51. Petrovskii, S.V., Malchow, H., 2001. Wave of chaos: new mechanism of pattern formation in spatiotemporal population dynamics. Theor. Popul. Biol. 59, 157–174. Petrovskii, S.V., Venturino, E., 2004. Patterns of patchy spatial spread in models of epidemiology and population dynamics. University of Torino, preprint # 20/2004 (available online: http://www.dm.unito.it/quadernidipartimento/quaderni.php). Petrovskii, S.V., Morozov, A.Y., Venturino, E., 2002a. Allee effect makes possible patchy invasion in a predator–prey system. Ecol. Lett. 5, 345–352. Petrovskii, S.V., Vinogradov, M.E., Morozov, A.Y., 2002b. Formation of the patchiness in the plankton horizontal distribution due to biological invasion in a two-species model with account for the Allee effect. Oceanology 42, 363–372. Petrovskii, S.V., Li, B.-L., Malchow, H., 2003. Quantification of the spatial aspect of chaotic dynamics in biological and chemical systems. Bull. Math. Biol. 65, 425–446. Petrovskii, S.V., Morozov, A.Y., Li, B.-L., 2005a. Regimes of biological invasion in a predator–prey system with the Allee effect. Bull. Math. Biol. 67, 637–661. Petrovskii, S.V., Malchow, H., Li, B.-L., 2005b. An exact solution of a diffusive predator–prey system. Proc. R. Soc. London A 461, 1029–1053. Petrovskii, S.V., Malchow, H., Hilker, F., Venturino, E., 2005c. Patterns of patchy spread in deterministic and stochastic models of biological invasion and biological control. Biol. Invasions, in press.

35

Sakai, A.K., Allendorf, F.W., Holt, J.S., Lodge, D.M., Molofsky, J., With, K.A., Baughman, S., Cabin, R.J., Cohen, J.E., Ellstrand, N.C., McCauley, D.E., O’Neil, P., Parker, I.M., Thompson, J.N., Weller, S.G., 2001. The population biology of invasive species. Annu. Rev. Ecol. Syst. 32, 305–332. Segel, L.A., Jackson, J.L., 1972. Dissipative structure: an explanation and an ecological example. J. Theor. Biol. 37, 545–559. Sherratt, J.A., 2001. Periodic travelling waves in cyclic predator–prey systems. Ecol. Lett. 4, 30–37. Sherratt, J.A., Lewis, M.A., Fowler, A.C., 1995. Ecological chaos in the wake of invasion. Proc. Natl Acad. Sci. USA 92, 2524–2528. Sherratt, J.A., Lambin, X., Thomas, C.J., Sherratt, T.N., 2002. Generation of periodic waves by landscape features in cyclic predator–prey systems. Proc. R. Soc. London B 269, 327–334. Sherratt, J.A., Lambin, X., Sherratt, T.N., 2003. The effects of the size and shape of landscape features on the formation of traveling waves in cyclic populations. Am. Nat. 162, 503–513. Shigesada, N., Kawasaki, K., 1997. Biological Invasions: Theory and Practice. Oxford University Press, Oxford. Swope, S., Carruthers, R., Anderson, G., Bubenheim, D., 2004. USDA and NASA collaborate to use hyperspectral imagery to detect invasive species through time and space. ESA Annual Conference, Portland, Oregon, July 30–August 6, 2004. Thomas, J., 1995. Numerical Partial Differential Equations: Finite Difference Methods. Texts in Applied Mathematics, Vol. 22. Springer, New York. Volpert, A.I., Volpert, V.A., Volpert, V.A., 1994. Travelling Wave Solutions of Parabolic Systems. American Mathematical Society, Providence. Wolf, A., Swift, H., Swinney, H., Vastano, J., 1985. Determining Lyapunov exponent from a time series. Physica D 16, 285–317.