## Abstract

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),(2.1)(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):(2.3)where *K* is the prey carrying capacity and *U*_{0} is a certain measure of the Allee effect. Depending on the sign of *U*_{0}, the Allee effect can be either ‘strong’ (*U*_{0}>0) or ‘weak’ (−*K*<*U*_{0}≤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*<*U*_{0}.

Regarding the *per capita* predator mortality, we assume that it is described by the following function:(2.4)where *M*, *d*_{0} 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 form(2.5)(2.6)Introducing dimensionless variables(2.7)we obtain from equations (2.5) and (2.6)(2.8)(2.9)where *β*=*U*_{0}*K*^{−1}, *k*=*κη*(*α**K*)^{−1}, *m*=*M*(*α**K*^{2})^{−1} and *δ*=*d*_{0}*αK*^{2}*η*^{−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)=*u*_{0}(*x*) and *v*(*x*, 0)=*v*_{0}(*x*). The form of *u*_{0}(*x*), *v*_{0}(*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:(2.10a)(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 system(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 it is somewhat less obvious. The following theorem addresses this issue.

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

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

The corresponding stationary values 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 which, after transformations, can be written in the form(2.18)Thus, the stationary values are given by(2.19)and(2.20)where plus and minus correspond to and , respectively.

Under the assumptions of theorem 1.1, *v*_{1,2}>0 and . Thus, the steady states (,), (,) lie inside and (,) lies outside The only remaining question is about the sign of and the corresponding steady state (,).

Obviously, the inequality(2.21)is equivalent to(2.22)Thus, the steady state (,) lies inside if the inequality (2.22) is true, and outside 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* *if* , *and it has three steady* *states in* *if* .

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 and, taking into account (2.10b), to(2.23)By virtue of equation (2.10b), *m*+1−*k*>0; thus, from (2.23), we arrive at(2.24)and then at(2.25)Taking into account that 0<*β*<1, from equation (2.10a), we obtain 0<*m*<1. Then, inequality (2.25) holds when(2.26)Considering (2.26) together with , we finally obtain:(2.27)Now, when *h*>0 and the system (2.11) has three steady states, from , we first arrive at(2.28)which is equivalent to(2.29)Considering it together with , we obtain that(2.30)However, for 0<*m*<1, , and inequality (2.30) is equivalent to(2.31)The proof is complete. ▪

Another important point is steady states' stability. Taking into account that the stationary values , 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 (,), 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 (,), a sufficient condition of its stability is given by the following theorem.

*Let conditions* (2.10) *be satisfied and*(2.32)*then the steady state* (,) *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 (,) 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 with the position of the hump. Under conditions (2.10), 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 (,) is stable for(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* , *the* *system* (2.8) *and* (2.9) *has the following exact solution*:(2.34)*where* *and* *are the steady states of the system* (2.11) *given by equations* (2.16) *and* (2.19), *respectively*, , ,(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 *n*_{1} and *n*_{2}, 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 that(2.36a)(2.36b)Obviously, *n*_{1}>*n*_{2} and *n*_{1} is always positive, whereas *n*_{2} can be either positive or negative depending on the parameter values:(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,(2.38)Thus, the wave propagating with the speed *n*_{1} is a travelling front connecting the steady states (0,0) and (,).

Owing to *n*_{1}>*n*_{2}, 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:(2.39)Therefore, in the large-time limit, the solution (2.34) describes a travelling front connecting the states (0,0) and (,) and propagating with the speed *n*_{2}. The direction of propagation depends on parameter values (see equation (2.37), and see also figures 3 and 4).

Curiously, the actual dynamics described by solution (2.34) is not exhausted by the two travelling fronts propagating with speeds *n*_{1} and *n*_{2}. 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):(2.40)Since(2.41)where , the solution in the crossover region apparently behaves as a travelling wave connecting the states (,) and (,) and propagating with the speed .

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 (,) 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 (,) 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 (,), 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 (,) 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 (,) 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 (,) to the stable homogeneous distributions at (0,0) and (,).

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 , for *x*>0.

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 (,) to the unstable state (,), 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:(3.1)(3.2)(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 6*a* shows sensitivity/robustness of the density profiles with respect to variation of the diffusion coefficients so that *D*_{pred}/*D*_{prey}≠1, other parameters are *m*=0.2, *β*=0.2, *k*=0.911, *δ*=12. It is readily seen that variation of *D*_{pred}/*D*_{prey} 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 *D*_{pred}/*D*_{prey} 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.

Although robustness of the wave shape with respect to reasonably small parameter variations seems to be a general property (cf. figure 6*a*–*c*), the wave speed appears to be rather sensitive to variations of the interaction parameters. Figure 6*b* 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 6*c* 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).

The exact solution (2.34), especially in its large-time asymptotical form (2.4), has a clear ecological/biological meaning. In the case *n*_{2}<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 *n*_{2}<0 the population front propagates to the region with low population density, which corresponds to species invasion, in the case *n*_{2}>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)*u*^{2}. 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, The system (2.8) and (2.9) is then virtually reduced to one equation:(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):(3.5)In dimensionless variables, *s*_{0}≤*s*_{1}≤*s*_{2} are the roots of the polynomial and minus corresponds to our actual choice of the conditions at infinity, that *u*(−∞,*t*)=*s*_{0}, *u*(+∞,*t*)=*s*_{2}. In the prey-only limit of the system (2.8) and (2.9), that is, in the case *v*≡0, we obtain:(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:(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:(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, *D*_{1}≪*D*_{2}. By introducing a small parameter *ϵ*=*D*_{1}/*D*_{2}, 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.

## Acknowledgments

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 equations(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 or for ∀*x*, *t*), then *σ* can have an almost arbitrary value with the only restriction or , respectively.

Substitution of equation (A 1) into (2.8) leads to the following equation:(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 system(A3)(A4)(A5)assuming that *z*_{x}≠0.

Similarly, substituting equation (A 1) into (2.9) and assuming *w*_{x}≠0, we obtain the following equations:(A6)(A7)(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 to(A9a)and(A9b)If we choose plus in both of the above equations (other possibilities are mentioned in §3), then equation (A 9) takes the form(A10a)and(A10b)which is consistent only in the case of the following relation between *μ* and *γ*:(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 equation(A12)Obviously, equation (A 12) becomes trivial under the following constraints on the parameter values:(A13)(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:(A15)The system (A 3)–(A 8) is now reduced to only three equations:(A16)(A17)(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*):(A19)Equation (A 19) is linear and its general solution has the form(A20)where the functions *f*_{0}, *f*_{1}, *f*_{2} still need to be determined and the eigenvalues *λ*_{1,2} are the solutions of the following equation:(A21)so that(A22)Taking into account equations (A 10a) and (A 20), we obtain(A23)where *g*_{0}(*t*) is a certain function.

To find the functions *f*_{0}, *f*_{1}, *f*_{2}, *g*_{0}, we substitute equations (A 20) and (A 23) into equation (A 16):(A24)Since *λ*_{1}≠*λ*_{2} (based on the assumption that ), then and are linear independent functions of *x*. Thus, equation (A 24) holds for any *x* only if functions *f*_{0}, *f*_{1}, *f*_{2}, *g*_{0} give a solution of the following system:(A25)From equation (A 25), we obtain(A26)where and the constants *C*_{1,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):(A27)(A28)where , *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, .

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 . Introducing new constants as(A29)and taking into account that , (cf. equations (2.16), (2.19) and (A22)), the solution (A 27) and (A 28) takes the following form:whereand *λ*_{1,2} are given by equation (A 22). That completes the proof of theorem 1.3.

## Footnotes

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

- © 2005 The Royal Society