An exact solution of a diffusive predator–prey system

Sergei Petrovskii, Horst Malchow, Bai-Lian Li


We consider a system of two nonlinear partial differential equations describing the spatio-temporal dynamics of a predator–prey community where the prey per capita growth rate is damped by the Allee effect. Using an appropriate change of variables, we obtain an exact solution of the system, which appears to be related to the issue of biological invasion. In the large-time limit, or for appropriate parameter values, this solution describes the propagation of a travelling population front. We show that the properties of the solution exhibit biologically reasonable dependence on the parameter values; in particular, it predicts that the travelling front of invasive species can be stopped or reversed owing to the impact of predation.

1. Introduction

Mathematical problems of population dynamics have been attracting considerable attention for years. Initiated by the classical papers by Lotka (1925), Volterra (1926), Fisher (1937) and Kolmogorov et al. (1937), this subject has now expanded into a vast scientific field at the interface between mathematics and ecology (Aronson & Weinberger 1978; Murray 1993; Shigesada & Kawasaki 1997). Currently, different mathematical techniques are used to describe the dynamics of populations (cf. Matsuda et al. 1992; Sato et al. 1994; Kot et al. 1996; Cannas et al. 2003). However, the classic approach based on differential equations still seems to be the most common (Holmes et al. 1994; Owen & Lewis 2001; Sherratt 2001; Wang & Kot 2001).

Among others, the issues of growing interest and increasing importance are those involving the spatial aspects of the population dynamics, particularly the problems inspired by biological invasions (Shigesada & Kawasaki, 1997; Frantzen & van den Bosch 2000; Fagan et al. 2002). Relevant mathematical studies led to the consideration of travelling population fronts (Hadeler & Rothe 1975; Aronson & Weinberger 1978; Dunbar 1983, 1986; Berezovskaya & Karev 1999; Leach et al. 2002) and related issues such as calculation of the speed of the fronts (Owen & Lewis 2001; Petrovskii et al. 2001; Wang & Kot 2001), the convergence of transient solutions to a stationary travelling wave (Hadeler & Rothe 1975; Larson 1978; also see Volpert et al. 1994 for a thorough review of more recent results), the processes preceding the propagation of stationary travelling front (Kawahara & Tanaka 1983; Ognev et al. 1995; Petrovskii & Shigesada 2001) and pattern formation induced by the front propagation (Yachi et al. 1989; Sherratt et al. 1995; Petrovskii & Malchow 2000; Petrovskii et al. 2001; Malchow & Petrovskii 2002).

In general, considerable progress has been made in the mathematics of bioinvasion during the decades after the pioneering works by Fisher (1937) and Kolmogorov et al. (1937). However, only a few exact solutions have been found, and these relate to the models consisting of a single equation (Ablowitz & Zeppetella 1979; Kawahara & Tanaka 1983; Newman 1980; Petrovskii & Shigesada 2001; Petrovskii & Li 2003). This situation is easily understood if one takes into account that population dynamics are usually nonlinear and most of the standard mathematical tools are not applicable. As for the case when the dynamics are described by a system of nonlinear partial differential equations (PDEs), we are not aware of any known exact solution with a clear biological meaning. However, exact solutions are of considerable use in theoretical studies because they provide a good opportunity to understand how the properties of a given model depend on parameter values, and they provide a convenient test for checking numerical simulations.

In spite of impressive recent advances made in the theory of nonlinear PDEs (Newell 1985), an appropriate change of variables still remains one of the main mathematical tools commonly used to obtain exact solutions (Calogero & Xiaoda 1991). According to Hopf (1950), a successful change of variables may sometimes linearize the equations under study (Kudryashov 1993; Ognev et al. 1995; Petrovskii 1999). In this paper, we extend this approach to a system of diffusive predator–prey equations and obtain an exact solution describing the propagation of travelling population fronts. We demonstrate that the solution possesses biologically meaningful properties. In particular, we show that a wave of invasion can be slowed down or reversed by an increased predation. This prediction is in good agreement with results of recent ecological studies (Fagan & Bishop 2000).

2. Main equations and theorems

According to a widely accepted approach, the spatio-temporal dynamics of a predator–prey system can be described by the equations (cf. Murray 1993; Holmes et al. 1994; Shigesada & Kawasaki 1997),Embedded Image(2.1)Embedded Image(2.2)where U and V are the densities of prey and predator, respectively, at position X and time T; the function f(U), is the per capita growth rate of the prey; the term r(U)V stands for predation; κ is the coefficient of food utilization; and g(V) is the per capita mortality rate of predator. Here, the first term on the right-hand side of equations (2.1) and (2.2) describes the spatial mixing caused either by self-motion of individuals (Skellam 1951; Okubo 1980) or by properties of the environment, for example, for plankton communities the mixing is attributed to turbulent diffusion (Dubois 1975; Okubo 1980). D is the diffusion coefficient, which we assume to be the same for prey and predator; the ecological relevance of this assumption will be discussed in §3.

For different species, functions f, r and g can be of different types. We assume that the prey dynamics are subjected to the Allee effect (Allee 1938; Dennis 1989), so that its per capita growth rate is not a monotonically decreasing function of the prey density, but possesses a local maximum. In this case, the standard parameterization is as follows (Lewis & Kareiva 1993):Embedded Image(2.3)where K is the prey carrying capacity and U0 is a certain measure of the Allee effect. Depending on the sign of U0, the Allee effect can be either ‘strong’ (U0>0) or ‘weak’ (−K<U0≤0; cf. Owen & Lewis 2001; Wang & Kot 2001). In this paper, we focus on the ‘strong’ Allee effect when the prey growth rate becomes negative for U<U0.

Regarding the per capita predator mortality, we assume that it is described by the following function:Embedded Image(2.4)where M, d0 and n are positive parameters. Function g(V) gives the so-called closure term because it is supposed not only to describe the processes taking place inside the predator population (such as natural or linear mortality, competition, possibly cannibalism, etc.), but also, virtually to take into account the impact of higher predators that are not included into the model explicitly (Steele & Henderson 1992). Different authors consider various functional forms for the closure term, particularly different n (cf. Edwards & Yool 2000 and references therein). It should be mentioned, however, that the accuracy of ecological observations is usually rather low, so that it is not often possible to give a reliable estimate of n. In this paper, for the sake of analytical tractability we restrict our consideration to the case n=2.

Finally, we assume that the predator shows a linear response to prey according to the classical Lotka–Volterra model, that is, r(U)=ηU. Then, equations (2.1) and (2.2) with (2.3) and (2.4) take the formEmbedded Image(2.5)Embedded Image(2.6)Introducing dimensionless variablesEmbedded Image(2.7)we obtain from equations (2.5) and (2.6)Embedded Image(2.8)Embedded Image(2.9)where β=U0K−1, k=κη(αK)−1, m=M(αK2)−1 and δ=d0αK2η−2 are positive dimensionless parameters, subscripts x and t stand for the partial derivatives with respect to dimensionless space and time, respectively. We consider equations (2.8) and (2.9) in an infinite space, −∞<x<∞, and for t>0. Functions u(x,t), v(x,t) are assumed to be bounded for x→±∞. To complete the problem, equations (2.8) and (2.9) should be provided with the initial conditions u(x, 0)=u0(x) and v(x, 0)=v0(x). The form of u0(x), v0(x) will be specified later.

Exact solutions of nonlinear PDEs are often ad hoc solutions obtained either for a specific form of nonlinearity (cf. Newell & Whitehead, 1969; Petrovskii & Shigesada 2001) or for certain restrictions on parameter values (Ablowitz & Zeppetella 1979), although it is difficult to say whether that stems from immaturity of contemporary theory of nonlinear PDEs or reflects certain intrinsic symmetries of given PDEs. In the rest of this paper, for the sake of analytical tractability, we assume the following relations between the equation parameters:Embedded Image(2.10a)Embedded Image(2.10b)Apparently, constraints (2.10) reduce the full four-dimensional parameter space of system (2.8) and (2.9) to a nonempty two-dimensional domain. The origin and the ecological meaning of the constraints (2.10) will become clear later.

Before proceeding to the analysis of system (2.8) and (2.9), we are going to give a brief illustration of the properties of the corresponding spatially homogeneous systemEmbedded Image(2.11)The question of primary importance is the existence of steady states, which as is usual in the case of autonomous systems of second order, are given by the intersection points of the zero-isoclines of the system. It is readily seen that equation (2.11) always has the trivial steady state (0,0) corresponding to species extinction and the two ‘prey only’ states (β,0), (1,0). As for the existence and the number of the steady states in the interior of the first quadrant, that is, in Embedded Image it is somewhat less obvious. The following theorem addresses this issue.

Let conditions (2.10) be satisfied and Embedded Image. Then, the number of steady states of system (2.11) inside Embedded Image depends on the sign of the following quantity:Embedded Image(2.12)Namely, there are two steady states in the case h<0 and three steady states in the case h>0.

Obviously, inside Embedded Image the isoclines of the system (2.11) are given by the equations (see figure 1)Embedded Image(2.13a)Embedded Image(2.13b)For any steady state (Embedded Image,Embedded Image), Embedded Image and Embedded Image are a solution of system (2.13). Having substituted equation (2.13a) to (2.13b), we obtainEmbedded Image(2.14)which, taking into account conditions (2.10), after standard although tedious algebraic transformations, can be written in the formEmbedded Image(2.15)Thus, possible stationary values Embedded Image are given byEmbedded Image(2.16)andEmbedded Image(2.17)It is readily seen that under the assumptions of theorem 1.1, all Embedded Image, i=1,…,4.

Figure 1

Example of the phase plane of the system (2.11) obtained for m=0.2, δ=14, k=0.933. Solid and dashed lines show the isoclines for prey and predator, respectively. Circles show the steady states and numbering of the states corresponds to the numbering in the text.

The corresponding stationary values Embedded Image can be obtained substituting equations (2.16) and (2.17) into (2.13a). However, we use another approach that appears to be less cumbersome. Having substituted equation (2.13b) to (2.13a), and taking into account (2.10), we obtain the equation for Embedded Image which, after transformations, can be written in the formEmbedded Image(2.18)Thus, the stationary values Embedded Image are given byEmbedded Image(2.19)andEmbedded Image(2.20)where plus and minus correspond to Embedded Image and Embedded Image, respectively.

Under the assumptions of theorem 1.1, v1,2>0 and Embedded Image. Thus, the steady states (Embedded Image,Embedded Image), (Embedded Image,Embedded Image) lie inside Embedded Image and (Embedded Image,Embedded Image) lies outside Embedded Image The only remaining question is about the sign of Embedded Image and the corresponding steady state (Embedded Image,Embedded Image).

Obviously, the inequalityEmbedded Image(2.21)is equivalent toEmbedded Image(2.22)Thus, the steady state (Embedded Image,Embedded Image) lies inside Embedded Image if the inequality (2.22) is true, and outside Embedded Image otherwise. This proves theorem 1.1. ▪

The above proof, however, leaves open the question of whether the conditions of theorem 1.1 can be met simultaneously; that is, can h actually be of different sign, taking into account that, by virtue of equation (2.10) in (2.12), only two parameters are independent? A close inspection of equation (2.12) along with equation (2.10) shows that the answer to this question is positive; the details are given by the following corollary.

The system (2.11) has two steady states in Embedded Image if Embedded Image, and it has three steady states in Embedded Image if Embedded Image.

Let us start with the case where h<0, and thus the system (2.11) has two steady states. Apparently, h<0 is equivalent to Embedded Image and, taking into account (2.10b), toEmbedded Image(2.23)By virtue of equation (2.10b), m+1−k>0; thus, from (2.23), we arrive atEmbedded Image(2.24)and then atEmbedded Image(2.25)Taking into account that 0<β<1, from equation (2.10a), we obtain 0<m<1. Then, inequality (2.25) holds whenEmbedded Image(2.26)Considering (2.26) together with Embedded Image, we finally obtain:Embedded Image(2.27)Now, when h>0 and the system (2.11) has three steady states, from Embedded Image, we first arrive atEmbedded Image(2.28)which is equivalent toEmbedded Image(2.29)Considering it together with Embedded Image, we obtain thatEmbedded Image(2.30)However, for 0<m<1, Embedded Image, and inequality (2.30) is equivalent toEmbedded Image(2.31)The proof is complete. ▪

Another important point is steady states' stability. Taking into account that the stationary values Embedded Image,Embedded Image appear as the solutions of fourth-order algebraic equations, a thorough investigation of this issue is very difficult. However, for the goals of this paper, it seems possible to restrict our consideration to the stability of only two steady states, that is, (0,0) and (Embedded Image,Embedded Image), where subscript ‘2’ corresponds to plus in equations (2.16) and (2.19). It is easy to see that (0,0) is always stable. As for (Embedded Image,Embedded Image), a sufficient condition of its stability is given by the following theorem.

Let conditions (2.10) be satisfied andEmbedded Image(2.32)then the steady state (Embedded Image,Embedded Image) is stable.

The conclusion of the theorem almost immediately follows from a more general fact that a steady state is stable when it arises as an intersection point of decreasing zero-isocline for prey (i.e. originated from the equation for prey) and increasing zero-isocline for predator. In terms of system (2.11), it means that (Embedded Image,Embedded Image) is stable when it is situated on the right of the hump of curve 1 (cf. figure 1). Thus, all we need is to compare Embedded Image with the position of the hump. Under conditions (2.10), Embedded Image is given by equation (2.16), and the maximum of the isocline for predator is situated at (β+1)/2. Then, taking into account (2.10b), the steady state (Embedded Image,Embedded Image) is stable forEmbedded Image(2.33)which is equivalent to assumption (2.32). That proves the theorem. ▪

Now we are ready to proceed to the main result of this paper concerning an exact solution of equations (2.8) and (2.9).

Under constraints (2.10) and Embedded Image, the system (2.8) and (2.9) has the following exact solution:Embedded Image(2.34)where Embedded Image and Embedded Image are the steady states of the system (2.11) given by equations (2.16) and (2.19), respectively, Embedded Image, Embedded Image,Embedded Image(2.35)and ϕ1,2 are arbitrary constants.

Theorem 1.2 can be proved by substituting (2.34) into equations (2.8) and (2.9). A simple way to arrive at solution (2.34) will be shown in §3. A somewhat more formal approach, based on an appropriate change of variables, that leads to solution (2.34) and also helps to provide an understanding of the origin of the relations (2.10), is given in Appendix A.

Now, since the solution (2.34) is likely to have a variety of ecological and biological applications, we are going to have a closer look at its properties. The form of solution (2.34) suggests that it describes propagation of two waves travelling with the speeds n1 and n2, correspondingly (see figure 2). Assuming, without any loss of generality, that λ1<λ2 (i.e. choosing minus for λ1 and plus for λ2 in equation (2.35)), we immediately obtain thatEmbedded Image(2.36a)Embedded Image(2.36b)Obviously, n1>n2 and n1 is always positive, whereas n2 can be either positive or negative depending on the parameter values:Embedded Image(2.37)In order to gain insight into the nature of the waves, let us consider the case when x and t, λ1ξ1≃1 or even λ1ξ1≪1 while λ2ξ2 is negative and large. For sufficiently small t, it can always be achieved by a proper choice of ϕ1 and ϕ2. Then, in this domain,Embedded Image(2.38)Thus, the wave propagating with the speed n1 is a travelling front connecting the steady states (0,0) and (Embedded Image,Embedded Image).

Figure 2

The densities of prey (solid line) and predator (dashed line) as given by the exact solution (2.34) shown for different time for parameters m=0.2, δ=34, k=1.029, ϕ1=30, ϕ2=−20. The homogeneous species distribution in the right-hand side of the domain corresponds to the stable steady state (Embedded Image,Embedded Image), and the homogeneous species distribution in the middle of the domain (cf. the upper two panels) corresponds to the unstable equilibrium (Embedded Image,Embedded Image). Transient dynamics at small times thus correspond to a decay of the unstable homogeneous distribution via propagation of the ‘partial’ fronts connecting (Embedded Image,Embedded Image) to (0,0) and to (Embedded Image,Embedded Image); see equations (2.38) and (2.40), respectively. Interaction of the partial fronts leads to formation of the travelling front connecting (0,0) and (Embedded Image,Embedded Image); cf. the panel at the bottom.

Owing to n1>n2, at any fixed point x variable ξ1 decreases at a higher rate than ξ2 and, regardless of the values of ϕ1 and ϕ2, for a sufficiently large time, the solution (2.34) can be approximated as follows:Embedded Image(2.39)Therefore, in the large-time limit, the solution (2.34) describes a travelling front connecting the states (0,0) and (Embedded Image,Embedded Image) and propagating with the speed n2. The direction of propagation depends on parameter values (see equation (2.37), and see also figures 3 and 4).

Figure 3

The densities of prey (solid line) and predator (dashed line) as given by the exact solution (2.34) shown for different time for ϕ1=−20, ϕ2=−10. Other parameters are the same as in figure 2. The homogeneous species distribution in the right-hand side of the domain corresponds to the stable steady state (Embedded Image,Embedded Image). The front propagation corresponds to species invasion.

Figure 4

The densities of prey (solid line) and predator (dashed line), as given by equation (2.34), shown for the case when the travelling front is reversed owing to the impact of predation (see comments in §3a). Parameters are m=0.2, δ=12, k=0.911, ϕ1=−20, ϕ2=−10. The homogeneous species distribution in the right-hand side of the domain corresponds to the stable steady state (Embedded Image,Embedded Image). The front propagation corresponds to species retreat.

Curiously, the actual dynamics described by solution (2.34) is not exhausted by the two travelling fronts propagating with speeds n1 and n2. Considering λ1ξ1≫1, in the crossover region where λ1ξ1 and λ2ξ2 are of the same magnitude we obtain from equation (2.34):Embedded Image(2.40)SinceEmbedded Image(2.41)where Embedded Image, the solution in the crossover region apparently behaves as a travelling wave connecting the states (Embedded Image,Embedded Image) and (Embedded Image,Embedded Image) and propagating with the speed Embedded Image.

Thus, in general, the propagation of the travelling front (2.39) can be preceded by the propagation of the ‘partial’ fronts (2.38) and (2.40) (see figure 2). They propagate toward each other, merge and create the single front (2.39). However, in cases where ϕ2 is much larger than ϕ1, the travelling fronts (2.38) and (2.40) may never be seen explicitly, and the exact solution (2.34) is well approximated by equation (2.38) for any t>0 (figures 3 and 4).

3. Discussion and conclusions

In this paper, we have considered the spatio-temporal dynamics of a predator–prey system described by two coupled nonlinear PDEs of ‘diffusion-reaction’ type. Using an appropriate change of variables, which arises as a generalization of an approach earlier successfully applied to single equations (Hopf 1950; Kawahara & Tanaka 1983; Kudryashov 1993; Ognev et al. 1995; Petrovskii 1999; Petrovskii & Li 2003), we have obtained an exact solution (2.34) describing propagation and interaction of travelling population fronts.

In the large-time limit, the solution describes a travelling front connecting the extinction state (0,0), which is always stable to the ‘upper’ coexistence state (Embedded Image,Embedded Image) which is stable under assumptions of theorem 1.2. The direction of the wave propagation depends on parameter values so that the front propagates toward the steady state that has a lower stability limit (see Murray 1989, pp. 297–302, for more details). Briefly, this situation becomes possible owing to the existence of the ‘intermediate’ unstable steady state (Embedded Image,Embedded Image) In our model, the intermediate state appears as a result of the strong Allee effect. In contrast, in a system without Allee effect (e.g. when the prey population exhibit logistic growth), the intermediate state does not exist, the extinction state (0,0) is always unstable and the direction of wave propagation does not depend on parameters. The front always propagates toward the region where the species are absent.

Regarding stability of the steady state (Embedded Image,Embedded Image), it should be mentioned that theorem 1.2 gives only a sufficient condition of stability, not a necessary one. Although we cannot prove the stability of (Embedded Image,Embedded Image) under less restrictive assumptions, the results of the numerical solution of equation (2.11) show that the actual range of its stability is wider than that given by theorem 1.2. In addition, the results of numerical simulations show that the state (Embedded Image,Embedded Image) is unstable in a wide range of parameter values. These results on the steady states' stability help us to understand better the transient dynamics described by the solution (2.34) at small times. Propagating toward each other, ‘partial’ travelling fronts (2.38) and (2.40) correspond to a ‘decay’ of the quasi-homogeneous species distribution at the unstable level (Embedded Image,Embedded Image) to the stable homogeneous distributions at (0,0) and (Embedded Image,Embedded Image).

Note that, since equations (2.8) and (2.9) are invariant with respect to the transformation x→(−x), one can expect the existence of the solution symmetrical to (2.34), that is, with similar travelling waves but propagating in opposite directions. Indeed, the procedure described in Appendix A leads to the solution with these properties if we choose minus in both equations (A 9a) and (A 9b). It is readily seen that other options (i.e. plus in equation (A 9a) and minus in equation (A 9b) or vice versa), lead to solutions of (2.8) and (2.9) which are not nonnegative and, thus, do not seem to have a clear ecological meaning.

The initial conditions corresponding to the exact solution (2.34) can be immediately obtained from equation (2.34) letting t=0. An important point, however, is that the meaning of solution (2.34) is not restricted to this specific case. Numerical simulations of system (2.8) and (2.9) show that the profile described by equation (2.34) actually appears, owing to solution convergence, for initial conditions from a wide class. Figure 5 gives an example of such convergence obtained in the circumstance that the initial conditions for equations (2.8) and (2.9) are chosen as piecewise-constant functions, that is, u(x,0)=0, v(x,0)=0 for x>0 and Embedded Image, Embedded Image for x>0.

Figure 5

Convergence of piecewise-constant initial conditions to the exact solution (2.34). Dashed curves show prey density u(x, t) as given by equation (2.34) for parameters m=0.2, δ=12, k=0.911, ϕ1=0, ϕ2=0.3. Solid curves show u(x, t) obtained in numerical simulations, equations (2.8) and (2.9) being solved by finite-difference method. Predator density v(x, t) exhibits similar behaviour.

In general, the rate of convergence to the solution (2.34) and the form of the transients can depend on the initial conditions. In particular, propagation of the partial waves connecting the stable steady states (0,0) and (Embedded Image,Embedded Image) to the unstable state (Embedded Image,Embedded Image), see equations (2.38) and (2.40), can be observed if the initial conditions include a sufficiently large domain where the species densities are at the unstable equilibrium, for example, for the following profile:Embedded Image(3.1)Embedded Image(3.2)Embedded Image(3.3)Apparently, the duration of the stage related to partial waves propagation depends on Δ. Discussion of ecological processes and circumstances that could lead to the initial conditions of that type is beyond the scope of this paper. For initial conditions of an arbitrary form, the partial waves are unlikely to be seen.

A point of interest is how the character of the travelling fronts (e.g. their speed and shape) depends on variations of the parameters, in particular, on variations violating conditions (2.10) and/or the assumption of equal species diffusivity (the meaning and ecological relevance of these constraints will be explained below). Figure 6 shows the profiles of the prey density obtained by means of numerical solution of the system (2.8) and (2.9) for different parameter values (dashed curves) as well as the exact solution (thick curves) at t=200. The thick curves on the left show the initial condition given by the exact solution for t=0. Since the predator density exhibits similar properties, only prey density is shown for the sake of brevity. Figure 6a shows sensitivity/robustness of the density profiles with respect to variation of the diffusion coefficients so that Dpred/Dprey≠1, other parameters are m=0.2, β=0.2, k=0.911, δ=12. It is readily seen that variation of Dpred/Dprey between 0.1 (the dashed curve on the right) and 1.0 (the thick curve in the middle) leads to less than 30% variation in the value of the wave speed. Variations of Dpred/Dprey between 1.0 and 3.0 (the dashed curve on the left) have somewhat larger impact and lead to about 40% variation in the wave speed. In both cases, the shape of the wave remains virtually unchanged.

Figure 6

Profiles of the prey density showing sensitivity/robustness of the solution to parameter variation. For each panel, the dashed curves show the solution of equations (2.8) and (2.9) obtained numerically at t=200. The thick curve in the middle shows the corresponding exact solution and the thick curve at the left-hand side shows the initial conditions: (a) variation of the species diffusivity, dashed curves from left to right, Dpred/Dprey=3.0, 2.0, 0.5, 0.1. The thick curve shows exact solution corresponding to Dpred/Dprey=1, other parameters m=0.2, β=0.21, k=0.911, δ=12; (b) variation of β, dashed curves from left to right, β=0.215, 0.210, 0.210, 0.205, 0.195, 0.190. The thick curve shows exact solution corresponding to β=0.2, other parameters as in case (a); (c) variation of k dashed curves from left of right: 0.895, 0.920, 0.935. The thick curve shows exact solution corresponding to k=0.911; other parameters as above.

Although robustness of the wave shape with respect to reasonably small parameter variations seems to be a general property (cf. figure 6ac), the wave speed appears to be rather sensitive to variations of the interaction parameters. Figure 6b shows the solution of equations (2.8) and (2.9) obtained numerically (dashed curves) for parameter β varied between 0.190 (right) and 0.215 (left), the thick curve in the middle shows the exact solution obtained for β=0.2. Other parameters are the same as above. Figure 6c shows the wave profiles obtained numerically for k varied from 0.895 (left) to 0.935 (right). The thick curve in the middle shows the exact solution corresponding to k=0.911 (other parameters are the same as above). It must be mentioned that larger parameter variations can change the type of the system dynamics at all, for example, change the regime of front propagation to a more complicated pattern of spread combined with spatio-temporal pattern formation (cf. Morozov et al. 2004). As a whole, a diffusive predator–prey system with the strong Allee effect for prey exhibits very rich dynamics and its sensitivity to parameter variations seems to be a typical property (Morozov et al. 2004; Petrovskii et al. 2004).

(a) Ecological relevance of the exact solution

Exact solution (2.34) was obtained under the assumption that the diffusion coefficient has the same value for the prey and predator. This assumption, however, is not very restrictive. Although predator often has higher diffusivity, a closer inspection of population communities reveals many trophical relations where diffusivity of prey and predator is of the same magnitude. One example is given by a plankton community where spatial mixing takes place mainly owing to turbulence (Okubo 1980), which has the same impact on phyto (prey) and zooplankton (predator). In terrestrial ecosystems, examples include lynx and hare, wolf and deer, and so on. In fact, predator's success is often reached not due to a faster motion but due to an optimal foraging strategy while its diffusivity must not be necessarily higher than that of prey (figure 7).

Figure 7

The map in the (β, δ) parameter plane. The domains above and below curve 1 correspond to species invasion and retreat, respectively, as it is given by relations (2.37). Curve 2 shows the boundary of the domain (above the curve) where exact solution (2.34) exists. The domain on the left of the straight line 3 corresponds to species retreat owing to increased Allee effect. The interior of the curvilinear triangle made by the lines 1, 2 and 3 corresponds to the species retreat because of predation (see comments in the text).

The exact solution (2.34), especially in its large-time asymptotical form (2.4), has a clear ecological/biological meaning. In the case n2<0, it describes the third ‘geographical’ stage of biological invasion when introduced exotic species, having established themselves in the new environment, start invading new areas (Shigesada & Kawasaki 1997). Remarkably, the direction of propagation of the population front can be different (see equation (2.37) and figures 3 and 4). While in the case n2<0 the population front propagates to the region with low population density, which corresponds to species invasion, in the case n2>0 the front propagates to the region with high population density and thus it corresponds to species retreat. Below, we will show that relations (2.37) actually describe an interplay between two different mechanisms, each of them can reverse the travelling front. One of these mechanisms is associated with the Allee effect and the other relates to the impact of predation.

In order to distinguish between the two mechanisms, let us start with a closer look at the relations (2.10) between the interaction parameters. In spite of their rather special form, they can be given a clear ecological interpretation. Consider first (2.10a). Both β and m have the same meaning. They give linear per capita mortality rate of prey and predator, respectively, in the case that all density-dependence phenomena can be neglected, for example, in the case of low population density.

To reveal the meaning of equation (2.10b), let us look at the reaction terms in the right-hand side of equations (2.8) and (2.9) from the point of (bio)mass vertical flow through the food web. Parameter k then quantifies the mass flow from the lower level (prey) to the upper level (predator) while δ quantifies the mass flow from predator to the higher trophical levels that are virtually taken into account by the last term in equation (2.9) (see the lines after equation (2.4)). Apparently, the larger k is, the larger the mass inflow, and the smaller δ is, the smaller the mass outflow. It means that the mass kept at the predator level can be quantified by the expression (k+δν) where ν is a positive parameter. On the other hand, the actual biomass production at the prey level is described by (β+1)u2. Constraint (2.10b) thus means that, up to the exact value of ν, the rate of biomass production at the lower level must be consistent with the rate of biomass assimilation at the upper level of this simple trophic web. It may mean that there should be certain similarity between the populations of prey and predator regarding their response to vertical biomass flow.

Altogether, along with the assumption of equal diffusion coefficients and equal mortality rates, it means that the exact solution (2.34) describes the dynamics of a population system where prey and predator are in some sense similar, for example, belong to the same taxonomic group. Indeed, one can observe that under constraints (2.10) variables u and v become proportional to each other, Embedded Image The system (2.8) and (2.9) is then virtually reduced to one equation:Embedded Image(3.4)Thus, solution (2.34) gives an extension of the Kawahara & Tanaka (1983) solution to the case of a system of two interacting species. We want to emphasize, however, that this extension is nontrivial. The point is that, in the single-species model with the population growth described by a cubic polynomial (see equation (3.4)), the speed of the front propagation is given by the equation (Murray 1993; Volpert et al. 1994):Embedded Image(3.5)In dimensionless variables, s0s1s2 are the roots of the polynomial and minus corresponds to our actual choice of the conditions at infinity, that u(−∞,t)=s0, u(+∞,t)=s2. In the prey-only limit of the system (2.8) and (2.9), that is, in the case v≡0, we obtain:Embedded Image(3.6)so that the front propagates toward the region where u≈0 for β<0.5 (species invasion) and toward the region where u≈0 for β<0.5 (species retreat); β=0.5 corresponds to the front with zero speed. Note that β is a dimensionless measure of the Allee effect; thus, an increased Allee effect can turn invasion to retreat.

Now, what can change in the presence of predator, that is, in case v≠0 identically? Having applied the comparison theorem for nonlinear parabolic equations (Volpert et al. 1994) to equation (2.8), it is readily seen that predation cannot turn retreat back to invasion. However, the impact of predation can turn invasion to retreat even when the invasion would be successful in the absence of predator. In the predator–prey system (2.8) and (2.9), the ‘turning’ relation between the parameters is:Embedded Image(3.7)(see equation (2.37)), so that the values of k greater than the one given by equation (3.7) correspond to species invasion, smaller k corresponds to species retreat (see figures 3 and 4).

Allowing for relations (2.10), the condition of species retreat takes the following form:Embedded Image(3.8)In the case β>0.5, inequality (3.8) describes species retreat resulting from a joint effect of two factors, that is, increased Allee effect (large β) and the impact of predation. As can be seen from comparison between equations (2.36b) and (3.6), the speed of retreat appears to be higher in the presence of predator. In the case β<0.5, however, inequality (3.8) describes species retreat, which must be essentially attributed to the impact of predation because in the absence of predator, the condition β<0.5 always leads to species invasion. Figure 6 shows the solution of inequality (3.8) in (β,δ) plane; the area inside the curvilinear triangle corresponds to the predator-caused retreat. Thus, based on the properties of exact solution (2.34), we have shown that when the population front is reversed in the predator–prey system (meaning that inequality (3.8) holds), the invasion could be successful in the absence of predator.

Factors that can either speed up or slow down the spread of invasive species have recently been a subject of considerable interest (cf. Kot et al. 1996; Sherratt & Marchant 1996; Savill & Hogeweg 1999; Frantzen & van den Bosch 2000; Campos et al. 2002; Fagan et al. 2002). In particular, the impact of predation was considered by Owen & Lewis (2001). However, these authors were more interested in the case when prey is virtually immobile, that is, D1D2. By introducing a small parameter ϵ=D1/D2, it was then possible to apply a singular perturbation analysis and to obtain an approximate condition of wave blocking. In contrast, in this paper, we are more interested in the case when the diffusivity of prey and predator is of the same order of magnitude. Under some additional constraints, it appears possible to solve the equations analytically and to obtain an exact condition of invasion wave blocking owing to predation. Predictions of our exactly solvable model (2.8) and (2.9) are in good agreement with the results of qualitative analysis by Owen & Lewis (2001) and with field observations by Fagan & Bishop (2000), both indicating that invasion can be stopped or reversed owing to predation.


This work was partially supported by the Russian Foundation for Basic Research under grants 03-04-48018 and 04-04-49649; by the Deutsche Forschungsgemeinschaft (DFG) under grant 436 RUS 113/631; by the US National Science Foundation under grants DEB-0080529 and DEB-0409984; by the University of California Agricultural Experiment Station; and by UCR Center for Conservation Biology.

Formal proof of theorem 1.3

Let us introduce new variables z(x, t) and w(x, t) by means of the equationsEmbedded Image(A1)where μ and γ are parameters and σ is a constant included into the denominator of (A 1) in order to avoid singularities because, for biological reasons, we are primarily interested in bounded solutions of system (2.8) and (2.9). If we restrict our analysis to the case where z+w is semi-bounded (i.e. either Embedded Image or Embedded Image for ∀x, t), then σ can have an almost arbitrary value with the only restriction Embedded Image or Embedded Image, respectively.

Substitution of equation (A 1) into (2.8) leads to the following equation:Embedded Image(A2)Because σ is (nearly) arbitrary and functions (z+w)j are linearly independent for different j (except for the trivial case z+w≡const.), equation (A 2) can hold only if the expressions in the square brackets equal zero identically. Thus, after some obvious transformations we arrive at the systemEmbedded Image(A3)Embedded Image(A4)Embedded Image(A5)assuming that zx≠0.

Similarly, substituting equation (A 1) into (2.9) and assuming wx≠0, we obtain the following equations:Embedded Image(A6)Embedded Image(A7)Embedded Image(A8)The idea of the further analysis is that the system (A 3)–(A 8) is over-determined and its consistency may only take place under certain constraints on the parameter values. In particular, equations (A 3) and (A 6) are equivalent toEmbedded Image(A9a)andEmbedded Image(A9b)If we choose plus in both of the above equations (other possibilities are mentioned in §3), then equation (A 9) takes the formEmbedded Image(A10a)andEmbedded Image(A10b)which is consistent only in the case of the following relation between μ and γ:Embedded Image(A11)Next, equations (A 4) and (A 7) become equivalent when their right-hand sides coincide. Taking into account (A 10), that leads to the equationEmbedded Image(A12)Obviously, equation (A 12) becomes trivial under the following constraints on the parameter values:Embedded Image(A13)Embedded Image(A14)It is readily seen that under condition (2.10b) each of equations (A 11), (A 13) and (A 14) is equivalent to the following:Embedded Image(A15)The system (A 3)–(A 8) is now reduced to only three equations:Embedded Image(A16)Embedded Image(A17)Embedded Image(A18)The number of equations still exceeds the number of variables. However, under condition (2.10a), equations (A 17) and (A 18) become identical and the system is reduced to only two equations.

Differentiating equation (A 16) with respect to x, substituting (A 17) and (A 18) and taking into account (A 10) and (A 15), we obtain the following equation for z(x, t):Embedded Image(A19)Equation (A 19) is linear and its general solution has the formEmbedded Image(A20)where the functions f0, f1, f2 still need to be determined and the eigenvalues λ1,2 are the solutions of the following equation:Embedded Image(A21)so thatEmbedded Image(A22)Taking into account equations (A 10a) and (A 20), we obtainEmbedded Image(A23)where g0(t) is a certain function.

To find the functions f0, f1, f2, g0, we substitute equations (A 20) and (A 23) into equation (A 16):Embedded Image(A24)Since λ1λ2 (based on the assumption that Embedded Image), then Embedded Image and Embedded Image are linear independent functions of x. Thus, equation (A 24) holds for any x only if functions f0, f1, f2, g0 give a solution of the following system:Embedded Image(A25)From equation (A 25), we obtainEmbedded Image(A26)where Embedded Image and the constants C1,2 are determined by the initial conditions.

Note that function z+w, as defined by equations (A 20), (A 23) and (A 26), is semi-bounded, which agrees with our earlier assumption (see the lines below equation (A 1)). Thus, our analysis has been consistent. Substituting equations (A 20), (A 23) and (A 26) into equation (A 1), we obtain the following solution of the diffusive predator–prey system (2.8) and (2.9):Embedded Image(A27)Embedded Image(A28)where Embedded Image, I=1,2. Thus, relations (2.10) allow the integration of the system at the cost of u and v being proportional to each other, Embedded Image.

It is not difficult to see that to avoid singularities in the right-hand sides of equations (A 27) and (A 28), it is necessary that Embedded Image. Introducing new constants asEmbedded Image(A29)and taking into account that Embedded Image, Embedded Image (cf. equations (2.16), (2.19) and (A22)), the solution (A 27) and (A 28) takes the following form:Embedded ImagewhereEmbedded Imageand λ1,2 are given by equation (A 22). That completes the proof of theorem 1.3.


    • Received November 27, 2003.
    • Accepted September 14, 2004.


View Abstract