Bistability induced by generalist natural enemies can reverse pest invasions
BBistability induced by generalist natural enemies can reverse pestinvasions
Sten Madec a , J´erˆome Casas b,c , Guy Barles a , Christelle Suppo ba LMPT, UMR CNRS 7350, Universit´e Fran¸cois-Rabelais de Tours, 37200 Tours - France b IRBI, UMR CNRS 7261, Universit´e Fran¸cois-Rabelais de Tours, 37200 Tours, France c Institut Universitaire de France
Key words : Reaction Diffusion system, long time dynamics, traveling wave, invasion process, biologi-cal control, prey-predator interaction, generalist predator : 35B40, 35K57, 92D25, 92D40, 92B99
Abstract.
Reaction-diffusion analytical modeling of predator-prey systems has shown that specialist natural ene-mies can slow, stop and even reverse pest invasions, assuming that the prey population displays a strong Alleeeffect in its growth. Few additional analytical results have been obtained for other spatially distributed predator-prey systems, as traveling waves of non-monotonous systems are notoriously difficult to obtain. Traveling waveshave indeed recently been shown to exist in predator-prey systems, but the direction of the wave, an essential itemof information in the context of the control of biological invasions, is generally unknown. Preliminary numericalexplorations have hinted that control by generalist predators might be possible for prey populations displayinglogistic growth. We aimed to formalize the conditions in which spatial biological control can be achieved by gene-ralists, through an analytical approach based on reaction-diffusion equations.The population of the focal prey — the invader — is assumed to grow according to a logistic function. Thepredator has a type II functional response and is present everywhere in the domain, at its carrying capacity,on alternative hosts. Control, defined as the invader becoming extinct in the domain, may result from spatiallyindependent demographic dynamics or from a spatial extinction wave. Using comparison principles, we obtainsufficient conditions for control and for invasion, based on scalar bistable partial differential equations (PDEs).The searching efficiency and functional response plateau of the predator are identified as the main parametersdefining the parameter space for prey extinction and invasion. Numerical explorations are carried out in the regionof those control parameters space between the super- and subsolutions, in which no conclusion about controllabi-lity can be drawn on the basis of analytical solutions.The ability of generalist predators to control prey populations with logistic growth lies in the bistable dynamicsof the coupled system, rather than in the bistability of prey-only dynamics as observed for specialist predatorsattacking prey populations displaying Allee effects. Analysis of the ordinary differential equations (ODEs) systemidentifies parameter regions with monostable (extinction) and bistable (extinction or invasion) dynamics. Bycontrast, analysis of the associated PDE system distinguishes different and additional regions of invasion andextinction. Depending on the relative positions of these different zones, four patterns of spatial dynamics can beidentified : traveling waves of extinction and invasion, pulse waves of extinction and heterogeneous stationarypositive solutions of the Turing type. As a consequence, prey control is predicted to be possible when space isconsidered in additional situations other than those identified without considering space. The reverse situation isalso possible. None of these considerations apply to spatial predator-prey systems with specialist natural enemies.The consideration of space in predator-prey systems involving generalist predators with a parabolic functionalresponse is thus crucial. a r X i v : . [ m a t h . A P ] J u l Introduction
Biological invasions are a major contemporary problem (Pimentel 2011; Garnier et al. 2012; Mistroet al. 2012; Potapov & Rajakaryne 2013; Wang et al. 2013; Savage & Renton 2013) for which fewsolutions are available, all of which are very costly. The use of natural enemies for the biological controlof invading insects is one of the most promising possibilities (Moffat et al. 2013; Li et al. 2014; Ye etal. 2014; Basnet & Mukhopadhyay 2014). As invasion is essentially a spatial process, the potential ofnatural enemies to stop or even reverse an invasion is of particular interest. The fundamental analyticalwork of Owen & Lewis 2001 showed that specialist predators could potentially slow, stop or reverse thespread of invasive pests. The reversal of pest spread by specialist predation requires a strong Allee effectfor the pest-only dynamics, defined as a negative growth rate for the prey population at low density. Inthe presence of a weak Allee effect, the predator can stop, but not reverse the wave of invasion. Theseconclusions have been confirmed in several other theoretical studies (Cai et al. 2014; Boukal et al.2007; Morozov & Petrovskii 2009).Generalist predators can also control prey effectively (Erbach et al. 2014; Chakraborty 2015). Theiruse could be promoted through conservation biological control programs without the need for exogenousspecialist natural enemies. Unfortunately, the role of generalist predators in the spatial control of theirprey has been much less studied than that of specialist predators, due to the intrinsic difficulties of havingto work with a system of equations rather than with a single scalar equation. However, two importantstudies have been carried out in this area : the analytical and comprehensive study of Du & Shi 2007,and the preliminary simulation study of Fagan et al. 2002. Both used the same model structure as we dohere, with logistic growth for both prey and predator populations, and a type II functional response forpredators. The convergence of these models was strengthened further by the in-depth analysis of Magalet al. 2008 in which space was not considered. It is difficult to use these models in a spatial context : thework of Du & Shi 2007 cannot deal with invasion and traveling waves, because it deals with a boundedspace. The numerical simulations of Fagan et al. 2002 are restricted to a few parameter values. Theyare, however, valuable, because they suggest conditions in which a generalist predator might be able tostop, and even reverse the invasion wave of a pest population displaying logistic growth. Fagan et al . alsoreported the results of field studies indicating that predators with diffusion coefficients higher than thoseof their prey are poor control organisms. The authors provided an explanation for this finding foundedon logical arguments, but without a firm mathematical foundation. This result has been confirmedby a few numerical simulations including space, as reported by Magal et al. 2008, revealing a strongdependence of system dynamics on the relative rates of diffusion of the prey and the predator. It isthus important to take space into account, by considering diffusion coefficients of both predator andprey. This conclusion accounts for the interest of scientists in questions of this type (Lewis et al. 2013;Hastings 2000; De Roos et al. 1991, 1998).There are therefore hopes that it might be possible to extend the conditions for the control of inva-sive prey organisms to (i) generalist predators and (ii) prey populations displaying growth patterns notdependent on the restrictive assumption of Allee effects. Such control approaches would have a majorimpact in the field, given the high degree of generalism obtained. The aim of this study was, thereforeto formalize the conditions in which spatial biological control can be achieved by generalists, throughan analytical approach based on traveling waves solutions of reaction-diffusion equations.Traveling wave solution describes a constant profile U moving through space at a speed c . Such wavesare often observed in nonlinear reaction-diffusion systems modeling various phenomena. They are par-ticularly suitable for describing the propagation of invasive fronts. In systems modeling a single species,described by a scalar equation, this type of solution is very well understood (Fischer 1937; Kolmogorovet al. 1937 and Volpert et al. 1994 for a complete theory). Two particular classes of equations canbe distinguished : monostable equations (like the Fisher-KPP equations) and bistable equations (oftenmodeling the Allee effect). In monostable equations, there is a minimal wave speed c ∗ such that, for2ny c ≥ c ∗ , a wave solution with speed c exists. In bistable equations traveling waves exist for a unique speed c ∗ . The sign of this speed c ∗ distinguishes between invasion or extinction of prey, which is a keyproperty for our purposes.For interactions of several species (described by a multidimensional system), the situation is muchmore complex. However, for some type of interaction, cooperation for instance, the system possessesa strong structural property, namely monotonicity. Essentially, this monotony makes it possible to usethe comparison principle, which is always possible for one-dimensional systems, and the theory is thencomplete (see Volpert et al. 1994). Unfortunately, our system, and prey-predator systems in general,do not have such a monotonous structure. This method is then unsuitable for monotonous systems andonly a few results have been published. One of the key reasons for this is as follows : when we search fortraveling wave solutions for a system with N equations, we obtain a system of N second order ordinaryequations that can be reformulated as a system of 2 N first order ordinary equations. In the scalar case( N = 1) , it is therefore possible to study trajectories in a plane, available using classical tools fortwo-dimensional dynamical systems. For several species ( N > N -dimensional space, which may be very difficult.Hence, the first rigorous results demonstrating the existence of traveling waves in prey-predatorsystems were based on a generalization, to the fourth dimension, of the classical shooting method in thephase plane (Dunbar 1984a,b). This approach has since been generalised (Huang et al. 2003; Xu &Weng 2012). However, all these studies simply investigate the mere existence of traveling waves. Theydo not determine the direction of the wave or the global dynamics for general initial conditions. Othermethods have recently been developed in similar models (Huang & Weng 2013; Ducrot & Langlais2012), but they are subject to the same limitations. A last approach is to use the degree theory (seee.g. Giovangigli 1990; Volpert et al. 1994) to obtain the existence of traveling waves. These homotopymethods may occasionally give some information on the speed c . Unfortunately, this needs additionalestimates which are very difficult to obtain here. We therefore required another method.The analysis provided in Magal et al. 2008 gives conditions for preys’ control by predators, butthis analysis was carried out largely without reference to space. Thus, we have extended the system ofMagal et al. 2008 by adding spatial diffusion. We find that the conditions for control are very differentfrom those for the system in which space is not considered. The conditions for prey extinction andinvasion are discussed in terms of two essential parameters : the encouter rate E and the handling time h . Increasing E clearly increases predator pressure. Conversely, increasing h decreases predator pressure.The paper is organized as follows. In section 2 we present the mathematical model and the mainresult of this work : theorem 2.1 describes invasion conditions for the ODE system and the theorems2.4 and 2.6 the invasion conditions for the PDE system. The mathematical results are completed bynumerical simulations in section 3. The results are discussed in section 4. The final section 5 is devotedto the mathematical proofs. We analyze a system of partial differential equations for a prey population with logistic growth, anda generalist predator population with logistic growth on alternative prey in the absence of the invadinghost. The functional response is of Holling type II. The prey-predator interactions are modeled by thefollowing partial differential equation system : ∂ t u = D u ∆ x u + r u (cid:16) − uK (cid:17) − Euv Ehu ,∂ t v = D v ∆ x v + r v (cid:16) − vK (cid:17) + γ Euv Ehu , x ∈ R , t ∈ R + u ( x,
0) = u ( x ) , v ( x,
0) = 1 (1)3ith : u ( t, x ) = prey density at time t and at point x. v ( t, x ) = predator density at time t and at point x. D u = diffusion rate of prey D v = diffusion rate of predators r = growth rate of prey r = growth rate of predators K = carrying capacity of prey K = carrying capacity of predators in absence of focal prey E = encounter rate h = handling time γ = conversion efficiency u ≥ D u , D v , r , r , K , K , E, h and γ are positive constant parameters.We carried out the following adimensionalization : t (cid:48) = r t ; x (cid:48) = x (cid:113) r D u ; u (cid:48) ( x (cid:48) , t (cid:48) ) = u ( t, x ) /K ; v (cid:48) ( x (cid:48) , t (cid:48) ) = v ( t, x ) /K d (cid:48) = D v /D u ; r (cid:48) = r /r ; E (cid:48) = EK /r ; h (cid:48) = r hK /K ; γ (cid:48) = γK /K ; α = γ (cid:48) r (cid:48) .Removing the sign ’ to simplify the notation, the system reads (cid:40) ∂ t u = ∆ x u + u (1 − u ) − Euv Ehu , x ∈ R , t ∈ R + ∂ t v = d ∆ x v + r (cid:16) v (1 − v ) + α Euv Ehu (cid:17) (2)with the initial conditions (cid:40) u (0 , x ) = u ( x ) ∈ [0 ,
1] ; lim x →−∞ u ( x ) = 1 ; lim x → + ∞ u ( x ) = 0 v (0 , x ) = 1 . (3) We distinguish two ways in which a predator can control the prey, one taking space into accountand the other not considering this factor (mathematical definitions are provided in definition 2.2.1).– The spatially uniform extinction results exclusively from local demographic processes and is inde-pendent of space.– The extinction wave is due to both demographic and diffusive processes and may take variousforms, from a traveling front to a pulse.Conversly, invasion is defined as prey survival and we distinguish two ways in which the prey can invade.– The spatially uniform invasion , which is independent of space.– The non-uniform invasion , described by various spatial dynamics, from Turing phenomena toinvasion waves.
Definition 2.2.1.
Let ( u ( x ) , v ( x )) be an initial condition verifying (3) and ( u ( t, x ) , v ( t, x )) be thecorresponding solution of (2) .– Extinction of prey occurs if ∀ x ∈ R , lim t → + ∞ u ( t, x ) = 0 . – Prey extinction is uniform if it is uniform with respect to x ∈ R , that is, if there exists a map φ ( t ) verifying ∀ x ∈ R , ∀ t > t , ≤ u ( t, x ) ≤ φ ( t ) and lim t → + ∞ φ ( t ) = 0 . – Prey extinction is non uniform if there is extinction but no uniform extinction.– Invasion of prey occurs if there is no extinction, that is if ∃ x ∈ R , lim sup t → + ∞ u ( t, x ) >
1. All our results remain true for various different initial conditions. The essential condition is that the solutions of thescalar systems we consider converge to traveling wave solutions. In particular, compact support may be allowed for u . SeeFife 1979 for a detailed discussion. .2.1 Analysis of the associated ODE system If space is not taken into account, system (2) may be rewritten as follows. ddt u = u (1 − u ) − Euv Ehu , t ∈ R + ddt v = r (cid:16) v (1 − v ) + α Euv Ehu (cid:17) . < u (0) = u ≤ ≤ v (0) = v (4)System (4) is well understood (Magal et al. 2008). Indeed, it is clear that there are always threetrivial stationary states : (0 ,
0) and (1 , , E >
1. Moreover, there are no more than three non-trivial positive steady states. Weare interested principally in the case
E >
1. In this case, there are either no or two stationary positivesteady states. If the two steady states exist, denoted ( (cid:98) u, (cid:98) v ) and ( u ∗ , v ∗ ) with (cid:98) u < u ∗ and (cid:98) v < v ∗ , then( (cid:98) u, (cid:98) v ) is always unstable and ( u ∗ , v ∗ ) is most often stable. In this case, there are two stable nonnegativesolutions, (0 ,
1) and ( u ∗ , v ∗ ) and the system is bistable.We are interested principally in the conditions for prey extinction. If E <
1, then (0 ,
1) is unstableand no extinction occurs. We are therefore interested only in the case
E >
1. Now, if
E >
1, there aretwo possibilities. In the non bistable case, (0 ,
1) is globally stable and there is extinction. In the bistablecase, provided that u is initially small enough, say u < µ for some 0 < µ <
1, then u ( t ) → t → + ∞ . Conversely, if u is initially large enough, say 0 < µ < u , then u ( t ) → u ∗ as t → + ∞ . Thus,in this case, the outcome — extinction or invasion — depends on the initial conditions. The followingresult provides an explicit statement of the above in the parameter space ( E ; h ) and is proven insection 5.1. Theorem 2.1.
Let
E > and α ≥ be fixed.(i) Existence of positive solutions.
There exists a unique h ∗ = h ∗ ( E, α ) such that– If h < h ∗ , then there is no positive stationary solution and there is extinction of prey for theODE system.– If h > h ∗ ( E, α ) , then there exist two positive solutions for the ODE system ( (cid:98) u, (cid:98) v ) and ( u ∗ , v ∗ ) with (cid:98) u < u ∗ .(ii) Stability of the solutions.
Let h > h ∗ ( E, α ) . The solution ( (cid:98) u, (cid:98) v ) is always unstable.Moreover, there exists a unique h ∗∗ ( E, α ) > h ∗ ( E, α ) such that– If h > h ∗∗ ( E, α ) then ( u ∗ , v ∗ ) is stable.– If h ∗ ( E, α ) < h < h ∗∗ ( E, α ) , the stability of ( u ∗ , v ∗ ) depends on r . It is unstable if r is smallenough and stable otherwise. Remark 1.
Our calculations show that the gap between h ∗ and h ∗∗ is very small, so that, roughlyspeaking, ( u ∗ , v ∗ ) is stable whenever it exists. However, if h belongs to the conditional stability zone, i.e. h ∈ ( h ∗ , h ∗∗ ) , and if r is very small, stability is lost and the system becomes excitable. This explains, inparticular, the presence of pulses for small values of r when dealing with spatial interactions (see section3). The map (
E, α ) (cid:55)→ h ∗ ( E, α ) has the following properties, as proved in section 5.1.
Properties 2.2.
Let α ≥ be fixed. The map E (cid:55)→ h ∗ ( E, α ) is increasing and one has the explicitlimits : lim E → h ∗ ( E, α ) = (cid:26) α if α < √ α if α ≥ E → + ∞ h ∗ ( E, α ) = 2 + 2 √ α. Let E ≥ be fixed. The map α (cid:55)→ h ∗ ( E, α ) is increasing and one has the explicit limits : h ∗ ( E,
0) = 1 E (cid:16) E − (cid:112) E ( E − (cid:17) := h ( E ) and lim α → + ∞ h ∗ ( E, α ) = + ∞ . Figure 1 illustrates the maps h ∗ and h ∗∗ and the possible outcomes for system (4).5 ∗ h ∗∗ encounter rate ( E ) h a nd li n g t i m e ( h ) (0 , Conditional Bistable Zone
Figure
1: Description of the dynamics of the ODE system (4) in the E − h plane. If E <
1, the control solution(0 ,
1) is unstable. In this zone there exists at least one positive steady state and the prey neverdisappears entirely. If
E >
1, then the control solution (0 ,
1) is always (locally) stable. Moreover, the
E > h ∗ curve, (0 ,
1) is the only non-negative steadystate and is a global attractor : this is a monostable zone. Above the h ∗ curve, there are two additionalpositive steady states, one of which is always unstable while the second, denoted ( u ∗ , v ∗ ), may be stableor unstable. Above the h ∗∗ curve, ( u ∗ , v ∗ ) is always stable. In this subzone, the asymptotic behaviordepends only on the initial conditions : this is a bistable zone. Between the h ∗ and the h ∗∗ curves, thestability of ( u ∗ , v ∗ ) depends on other parameters : this is a conditional bistable zone. For illustrativepurpose, the size of this last subzone has been considerably increased. We wish to identify the parameter conditions required to obtain prey extinction in the PDE system(2). A simple stability analysis shows that, if
E <
1, then invasion occurs in the PDE system. If
E > u ∗ , v ∗ ) and thatthe initial condition u ( x,
0) is close to u ∗ at some places x and close to 0 at other places. Since both u ∗ and 0 are stable, the demographic phenomena lead to an agregation near u ∗ and an agregation near0. However, diffusion allows individuals to move around in space, so one of 0 or u ∗ may be the finalglobal attractor. In other words, there may be a (stable) traveling wave joining u ∗ to 0. The directionof this wave, given by the sign of the speed of the wave, indicates whether extinction or invasion occurs.However, there are difficulties associated with this argument.– There can be no homogeneous stationary solution of (2), only stable heterogeneous positive sta-tionary solutions. In other words, it is possible that h < h ∗ ( E, α ) without control occurring.– Even in the case of bistability ( h > h ∗ ( E, α )), the bistable system (2) is neither competitive norcooperative. Little theoretical knowledge is available concerning the occurrence of traveling wavesin such systems, with even less known about the stability and direction of the wave.Using super and subsolutions, we show here how to obtain the conditions sufficient (but not neces-sary) for extinction and for invasion, based on well known scalar bistable PDEs. Roughly speaking, let( u ( t, x ) , v ( t, x )) be the solution of (2). If we find a positive constant v such that, for any ( t, x ) ∈ R + × R , v ( t, x ) ≥ v then it comes ∂ t u ( t, x ) − ∆ x u ( t, x ) ≤ u ( t, x )(1 − u ( t, x )) − Eu ( t, x )1 + Ehu ( t, x ) v, t > , x ∈ R .
2. It suffices that this condition occurs for t > t for some t > u be the solution of (cid:40) ∂ t u ( t, x ) − ∆ x u ( t, x ) = u ( t, x )(1 − u ( t, x )) − Eu ( t,x )1+ Ehu ( t,x ) v,u (0 , x ) ≥ u (0 , x ) . (5)The comparison principle implies that u ( t, x ) ≥ u ( t, x ). Now, if u ( x, t ) → t → + ∞ then u ( t, x ) → t → + ∞ and extinction of prey occurs (see figure 2-(a)). Moreover, if u ( x, t ) = φ ( t )does not depend on x , then there is aspatial control (see figure 2-(b)).Conversely, if we can identify a positive constant v such that v ( x, t ) ≤ v , then it comes ∂ t u ( t, x ) − ∆ x u ( t, x ) ≥ u ( t, x )(1 − u ( t, x )) − Eu ( t, x )1 + Ehu ( t, x ) v, t > , x ∈ R . Now, define u ( t, x ) as the solution of (cid:40) ∂ t u ( t, x ) − ∆ x u ( t, x ) = u ( t, x )(1 − u ( t, x )) − Eu ( t,x )1+ Ehu ( t,x ) v,u (0 , x ) ≤ u (0 , x ) . (6)The comparison principle implies that u ( t, x ) ≤ u ( t, x ). It follows that if for some x ∈ R , lim sup t → + ∞ u ( x, t ) > t → + ∞ u ( t, x ) > E , α and h , for extinction or invasion to occur. All theorems are proven in section 5. Webegin with a sufficient condition for uniform extinction. u ( x, t ) u ( x, t ) space ( x ) p r e y d e n s i t y ( u ) u ( x, t ) φ ( t ) space ( x ) p r e y d e n s i t y ( u ) (a) Extinction (b) Uniform extinction Figure
2: (a) Sufficient condition for extinction. The solution u ( t, x ) is majored by a supersolution u ( t, x ). Iflim t → + ∞ u ( t, x ) = 0, then lim t → + ∞ u ( t, x ) = 0 and there is extinction.(b) Sufficient condition for uniform extinction. The solution u ( t, x ) is majored by a supersolution φ ( t )which does not depend on x . If lim t → + ∞ φ ( t ) = 0, then lim t → + ∞ u ( t, x ) = 0 uniformly in x and there isuniform extinction. Theorem 2.3.
Let
E > and define h ( E ) = 1 E (cid:16) E − (cid:112) E ( E − (cid:17) . If h < h ( E ) then there is uniform extinction. In other words, for any initial condition verifying (3) ,there exists φ ( t ) ≥ such that any solution ( u ( t, x ) , v ( t, x )) of (2) verifies ∀ x ∈ R , ≤ u ( t, x ) ≤ φ ( t ) and lim t → + ∞ φ ( t ) = 0 . ( x, t ) u ( x, t ) space ( x ) p r e y d e n s i t y ( u ) Invasion
Figure
3: Sufficient condition for invasion. The solution u ( t, x ) is minored by a subsolution u ( t, x ). Iflim sup t → + ∞ u ( t, x ) >
0, then lim sup t → + ∞ u ( t, x ) > If h > h ( E ), then there can be invasion or extinction. The following theorem gives a sufficientcondition for extinction to occur. Theorem 2.4.
Let
E > be fixed and let ( u, v ) be the solution of (2) - (3) . Define v = 1 and let u be asolution of (6) together with u (0 , x ) = u (0 , x ) . There exists a unique h − = h − ( E ) > h ( E ) , such that– If h < h − ( E ) , then ∀ x ∈ R , lim t → + ∞ u ( t, x ) = 0 – If h > h − ( E ) , then ∀ x ∈ R , lim t → + ∞ u ( t, x ) = µ where µ = µ ( E, h ) is a positive scalar.As a consequence, if h < h − ( E ) there is extinction of prey. The map E (cid:55)→ h − ( E ) verifies the following properties proved in section 5.3. Properties 2.5.
The map E (cid:55)→ h − ( E ) is increasing and admits the following explicit limits : lim E → h − ( E ) = 1 ; lim E → + ∞ h − ( E ) = 163 . Our last result gives a sufficient condition for invasion to occur.
Theorem 2.6.
Let
E > and α ≥ be fixed and let ( u, v ) be the solution of (2) - (3) . Define v =1 + α E Eh and let u be a solution of (6) together with u (0 , x ) = u (0 , x ) . There exists a unique h + = h + ( E, α ) > h ∗ ( E, α ) such that– If h < h + ( E, α ) , then ∀ x ∈ R , lim t → + ∞ u ( t, x ) = 0 – If h > h + ( E, α ) , then ∀ x ∈ R , lim t → + ∞ u ( t, x ) = µ where µ = µ ( E, h ) is a positive scalar.As a consequence, if h > h + ( E, α ) there is invasion of prey. Finally, the following result specifies the behavior of the map h + . Properties 2.7.
The maps E (cid:55)→ h + ( E, α ) and α (cid:55)→ h + ( E, α ) are increasing. For any E > , h + ( E,
0) = h − ( E ) and lim α → + ∞ h + ( E, α ) = + ∞ . Finally, α ≥ being fixed, one has the explicit limit lim E → + ∞ h + ( E, α ) = 83 (cid:32) (cid:114) α (cid:33) . Remark 2.
The limit lim E → h + ( E, α ) remains unknown. However, it can be proved that this limit existsand is greater than lim E → h ∗ ( E, α ) . { ( E, h ) , E > , h − ( E ) < h < h + ( E, α ) } ,which we will refer to as the ‘transition zone’, it is not possible to draw any conclusions concerningwhether prey invasion or extinction is likely to occur. Indeed, this zone can be separated into twosubzones, according to the parameters values :Zone I = { ( E, h ) , E > , max ( h ∗ ( E, α ) , h − ( E )) < h < h + ( E, α ) } , Zone II = { ( E, h ) , E > , h − ( E ) < h < h ∗ ( E, α ) } . In Zone I, our numerical simulations show non-monotonous traveling waves. In zone II, simulations showvarious types of behavior, including pulse and even heterogeneous positive stationary solutions. Thisphenomena are discussed in the section 3. h ∗ h − h + h encounter rate ( E ) h a nd li n g t i m e ( h ) Uniform extinctionExtinctionInvasion U n i f o r m i n v a s i o n Transition Zone (I)Transition Zone (II)
Figure
4: Description of the dynamic of the PDE system (2) in the E − H plan. If E <
E > h < h ( E ), extinction for h ( E ) < h < h − ( E ) and invasionfor h + ( E, α ) < h . The zone between h − and h + is called the transition zone. This transition zone issplitted into two subzones : Zone I and Zone II, separated by the h ∗ curve. In these two zones, bothextinction or invasion of prey may occur due to various spatial phenomena. Numerical study of the transition zone α . The mathematical results above demonstrate the influence of the parameters E and h , and, indirectly,that of the conversion rate α , on the long-term behavior of the system. More precisely, when E > h . Indeed, we can define two values h − = h − ( E ) < h + = h + ( E, α ) (see theorems 2.4 and 2.6). Extinction occurs if h < h − and invasionoccurs if h > h + . When h ∈ ( h − , h + ) we observe richer dynamics, which may depend on other factors.We refer to this zone as the transition zone . Note that, as h − is not dependent on α and h + is anincreasing function of α (proposition 2.7), the size of this transition zone increases with increasing α . Figure
5: Computation of h ∗ ( E, α ), h − ( E ) and h + ( E, α ) for four values of α . The so-called transition zone liesbetween the red curve h + and the blue curve h − . The transition zone increases with increasing α . Thistransition zone can be split into two subzones separated by h ∗ (black line). A first clue to the possible dynamics in the transition zone is provided by an understanding of thedynamics of the ODE system (4) described in the theorem 2.1. The dynamic of (4) is essentially de-pendent on the position of h relative to h ∗ = h ∗ ( E, α ). When h < h ∗ , there is no positive stationarysolution, whereas for h > h ∗ there are two positive stationary solutions, one of which, the larger of thetwo, is (nearly always) stable.The position of h ∗ relative to h − and h + provides a first description of the transition zone. Byvirtue of proposition 2.2, one gets the following. We always have h ∗ < h + but the position of h ∗ relativeto h − is dependent on α . On the one hand, from the facts that h ∗ (1 , α ) > h − (1) for α > h ∗ ( E,
0) = h ( E ) < h − ( E ) for E >
1, we deduce that h − > h ∗ for large enough values of E and small
3. The dynamics generally also depends on a quantity h ∗∗ , defined in the theorem 2.1, slightly greater than h ∗ that isnot taken into account here for the sake of simplicity. α . On the other hand, h ∗ is an increasing function of α tending to + ∞ . We obtainthat h ∗ ( E, α ) > h − ( E ) for large values of α and any E > Remark 3.
When h − > h ∗ , which may occur for sufficiently small values of α , we see that taking spaceinto account automatically increases the potential of extinction of preys. The transition zone can thus be separated into two subzones : one in which h < h ∗ (Zone II) and onein which h > h ∗ (Zone I). Figure 5 sums up this discussion. As we will see below, both extinction andinvasion are possible in each of these zones, but the phenomena at work differ considerably, accordingto whether h < h ∗ or h > h ∗ . These phenomena are studied in more detail below, using a numericalapproach. Figure
6: Numerical computation of h crit in the E − h space for the fixed values of α = 4, d = 1 and r = 1.The transition zone h ∈ ( h − , h + ) is split into two subzones separated by h crit . Extinction occurs below h crit and invasion occurs above. r and d . Our numerical analysis shows that both invasion and extinction are possible in the transition zone.When space is taken into account we see that h ∗ does not separate the zone of invasion from that ofextinction. These two zones are, indeed, separated by a new critical value of the handling time denoted h crit ∈ ( h − , h + ), which is dependent on E and α , of course, but also on the relative rates of growth ( r )and diffusion ( d ) of the predator population : h crit = h crit ( E, α, r, d ) . As expected, when h > h crit , prey invasion is observed, whereas extinction is observed when h < h crit .Thus, higher values of h crit are associated with more effective predation and thus with less effectiveinvasion by the prey.Figure 6 completes the theoretical scheme represented in figure 4, by presenting an example of thecurve ( E, h crit ) for a particular selection of values for the parameters α , r and d in the E − h plan.11ike h ∗ and h + , h crit increases with both E and α . This naturally translates into the fact that,higher values of E increases the chance of meeting between predators and prey and that at higher α values, the predator is able to make greater use of the prey and can therefore eliminate it.Given the multiple dependence of h crit on different parameters, figure 6 can only represent a parti-cular case, chosen for its simplicity. Figure 7 shows the relationships between h crit and d for fixed valuesof E and α and for various values of r . We see that h crit increases with r . This translates into the factthat for small values of r , the predators growth rate is small and their effectiveness reduced. Conversly, h crit (essentially) decreases with increasing d . This is due to the fact that for large d , predators spreadinto a zone in which prey are not present, resulting in a weakening predation. Figure
7: Numerical computation of h crit with respect to d for four values of r with E = 2 and α = 4. Fora given value of the parameters, there is invasion of prey if h > h crit and extinction if h < h crit .Therefore, the larger h crit is, the greater is the potential for extinction of prey. We see that h crit isincreasing in r and (essentially) decreasing in d . Thus, an increase of d or a decrease of r decrease theimpact of the predation on invasive prey. Remark 4.
For small values of r and intermediate values of d , predators may increase their effectivenessby increasing d . In that case, the predator growth rate being small, predators density remains high for along time even if prey are absent. Now, if d is large enough but not large, predators may spread into azone in which prey are not present and remain there at a high density and long enough to stop the preyto invade. This phenomenon may enable predators to form a barrier to prey’s movement, preventingthereby prey propagation. For too large values of d , the loss in predators effectiveness due to movementis too strong and the above phenomenon does not hold any more. This explains why, for small values of r , h crit first increases and then decreases with increasing d . < h crit : Extinction h > h crit : Invasion h < h ∗ : Zone II PULSE ( r (cid:28)
1) TURING ( d (cid:29) h = 5 .
35 ; d = 1 ; r = 0 . h = 5 .
35 ; d = 100 ; r = 1 h > h ∗ : Zone I ETW ITW h = 5 . d = 1 ; r = 1 h = 6 ; d = 1 ; r = 0 . Table
1: Summary of the four different situations involving space in the transition zone. Here, E = 2 and α = 4which yields h ∗ ≈ .
4. The X axis represents space and the Y axis the concentration of species. Theblue curve (top of each graph) represents the concentration of predators in space while the black curveis the concentration of prey. Extinction and Invasion Traveling Waves are abbreviated ETW and ITWrespectively.
We know that extinction occurs when h < h crit , whereas invasion occurs when h > h crit . It should beborne in mind that the existence or absence of non-trivial solutions that are homogeneous over space aredependent on the position of h with respect to h ∗ . Consequently, the processes at work during extinctionor invasion are highly dependent on these position.We will now describe the different dynamics occurring in the transition zone, summarized in table1. See section 2.2 for a precise definition of the various quantities described here. A] Invasion ( h > h crit ) – Turing instabilities : h < h ∗ and d (cid:29) . As d increases, predators spread out, moving into areasfrom which the prey is absent, leading to a decrease in the size of the predator population. Thisphenomenon leads to a decrease in predator density throughout the space occupied by the predator,allowing the prey to survive in certain zones. We thus obtain a periodic distribution in space anda constant distribution over time of the densities of the prey and predator. Mathematically, thisphenomenon is described by a Turing bifurcation.– Invasion traveling waves (ITW) : h > h ∗ . The invasion is described simply by an invasiontraveling wave : a wave of propagation linking the two stable solutions ( u ∗ , v ∗ ) and (0 ,
1) in thedirection of the positive solution ( u ∗ , v ∗ ). When r is large this traveling wave is monotone. By13ontrast, when r (cid:28) r (cid:28) d is not too large, we observethat h crit is approximately equal to h ∗ (see figure 7). This indicates that there is an invasiontraveling wave if a positive solution exists. In this case, ahead of the front, v = 1 and, since r (cid:28) h > h − , the prey invades the space whenthe predator is at concentration 1. This leads to front advancing. Behind the front, the predatorhas had sufficient time to increase the size of its population and, therefore, to decrease the size ofthe prey population. As h > h ∗ , the population of the prey decreases towards the positive solution u ∗ and we observe a non-monotonous invasive traveling wave. B] Extinction ( h < h crit ) – Pulse : h ∈ [ h − , h ∗ ] and r (cid:28) . Since r (cid:28)
1, the front of the wave is similar to the ITWdescribed above. However, as h < h ∗ , here is no homogeneous positive solution u ∗ . The preypopulation therefore decreases to zero behind the front, whereas it continues to advance in aheadof the front. We thus obtain a pulse.– Extinction traveling wave (ETW) : h > h ∗ . This corresponds to the simplest case describedabove. We observe a propagation wave linking the two stable solutions ( u ∗ , v ∗ ) and (0 ,
1) in thedirection of the control solution (0 , E − h plane for fixed values of α , d and r . It furthermorespecifies the possible dynamics in each zone. Figure
8: Different spatial dynamics in the transition zone. To ensure that all possible situations are represented,we choose α = 0 . r = 0 .
01 and d = 10. Note that traveling waves of extinction (ETW) and travelingwaves of invasion (ITW) may be obtained outside the transition zone. We have shown that invasion occurs if
E <
1. If predators do not encounter their prey they cannotcontrol them. If
E >
1, then extinction or invasion can occur, depending on the parameters h , E , α r . Uniform extinction occurs for h < h ( E ) and extinction occurs for h ( E ) < h < h − ( E ). Invasionoccurs for h + ( E, α ) < h . Thus, if h increases, we move from a zone of extinction without a considera-tion of space to a zone of extinction requiring a consideration of spatial aspects and then to a zone ofinvasion. For intermediate values of E , the zones of extinction increase with increasing E resulting in ahigher potential of extinction, as h , h ∗ , h − are increasing functions of E . When E is large, the zones ofcontrol do not depend on E any more because h , h − , h + have finite limits when E → + ∞ . Thus, E canonly play a role in prey control if it takes intermediate values and if h is not too large. In summary, forlow values of E ( E <
1) or high values of E or h , the outcome of the interaction (extinction or invasion)is independent of E .There is furthermore a transition zone splitted in two subzones, with various spatio-temporal phe-nomena and wherein both extinction and invasion can occur. The size of this transition zone greatlyincreases when the conversion rate α increases. Depending on the relative positions of these two zoneswith regard to the zones of extinction and of invasion, four spatial dynamics were identified : extinctionand invasion traveling waves, extinction pulse waves and heterogeneous stationary positive solutions ofthe Turing type. We have shown that an increase in E increases the potential of extinction while an increase in h increases the potential of invasion. This translates the fact that a highly effective predator does havea high encounter rate and a small handling time. Furthermore, since h + is an increasing function of α , an increase in α decreases the potential of invasion and increases the size of the transition zone,which in turn increases the potential for the system to have complex dynamics. Finally, an increase ofthe diffusion rate d and a decrease of the amplitude of the predators growth rate r both increase thepotential of invasion of prey. Thus, a generalist predator loses its effectiveness to exterminate invasiveprey if it diffuses too fast or if it has a too slow dynamics.The above results are stated in term of adimensionalized parameters (see section 2.1). By choosing theappropriate spatio-temporal variables, we may define D u = r = 1. Thus, the biological interpretationsof d and r are accurate. Conversely, the definitions of the searching efficiency, E = EK , the handlingtime h = h K K and the conversion rate α = γr complicate the biological interpretation of these threeparameters. Thus, in addition to the above discussion about the influence of E and h , we now discussour results in terms of the other biological variables : K , K , γ and r . The parameter h being increasingin the carrying capacity K of prey, prey with high carrying capacity show a high risk of being invasive.Conversly, h and E are respectively decreasing and increasing in the carrying capacity of predators K .Thus, predators with high carrying capacity have a high potential to control prey invasion. Otherwise γ < γ > γ and a decrease in theamplitude of the predator growth rate r yield an increase of α . Therefore, predators with a preference forthe invasive prey or predators with a slow dynamics might display a complex dynamics. In particular,the likelihood of the system to exhibit a pulse wave is then important. The model analyzed here was studied without taking space into account, except for a numericalexploration in the discussion, in an article by Magal et al. 2008. As explained in the introduction,models of identical structure have been proposed independently by Fagan et al. 2002 and Chakraborty2015. We will now discuss our results in the context of these previous studies. Adding a spatial componentto predator-prey systems makes any prediction about the controllability of the system difficult, as itthen depends on the values of several parameters. The comparison between situations with and without15he consideration of space is epitomized by the distinction between h ∗ , separating parameter regionsof mono- and bistability in the ODE system, and h crit , separating parameter regions of invasion andextinction in the PDE system. We will now focus on the case of E >
1, as values of
E < h < h ∗ , as 0 is a global attractor. This is stilltrue in situations in which space is taken into account, if h < h where h is smaller than h ∗ . When h isbetween h and h ∗ , 0 is only a local attractor, so it is not possible to state that control is always attained.In this respect, adding consideration of space decreases the potential for control. Furthermore, whenspace is not taken into account, there is either extinction or invasion when h > h ∗ ( E, α ), dependingon initial conditions. Incorporating consideration of space changes the region where invasion occurs,for any values of the other parameters and for appropriate initial conditions, into h > h + ( E, α ), with h + > h ∗ . Thus, the consideration of space reduces the size of the zone wherein the invasion is certainand is detrimental to the invading prey. Finally, the relative levels of predator and prey diffusion alsodetermine the potential for control. Our model shows that control is increased by predators being lessmobile than prey. If predator mobility levels are too high, the predators become to thinly spread on theground. For similar reasons, too high a level of prey mobility leaves the prey vulnerable to predators. Thisis entirely consistent with the experimental findings of Fagan et al. 2002. In conclusion, taking spaceinto account can lead to an increase or a decrease in the controllability of invading prey by predators ; theaddition of space to the model has no generic implication for considerations of predator-prey dynamics(see also Lam & Ni 2012; Braverman et al. 2015). The originality of this study lies in its consideration of a generalist predator in a spatial context.When studying generalist predators, it is common practice to assume that the functional response isof type III, due to switching between prey species (Erbach et al. 2014; van Leuven et al. 2007, 2013;Morozov & Petrovskii 2013). However, this approach is not mandatory, and other works (Basnet &Mukhopadhyay 2014; Krivan & Eisner 2006; Hoyle & Bowers 2007) have considered a type II functionalresponse. Altering our model to include a type III functional response would be very costly in terms ofunderstanding, because such responses lead to a loss of bistability. Its derivative would be null withoutprey, so some of our demonstration would fail and the analytical complexity would be greatly increased.However, traveling waves for specialist predator with type III functional response are known to exist (Li& Wu 2008) which indicates that our result may be extended to this case.The complexity of analytical studies of spatial predator-prey interactions lies in the reaction termsbeing of alternative signs in the equations, making the study of the systems of equations essential(Dunbar 1984b; Huang et al. 2003; Huang & Weng 2013). Other interactions, such as competitionof two species (all negative) and symbiosis (all positive), are simpler, as their studies are similar tothe study of a single equation (Volpert et al. 1994; Alzahrani et al. 2012). This accounts for theslow scientific progress in this otherwise highly relevant topic. However, several major results have beenobtained in recent decades, including those of the fundamental work of Owen & Lewis 2001. The findingof Owen and Lewis that predators can slow, stop, and even reverse invasion by their prey was basedon the bistability of the prey-only dynamics of systems consisting of specialist predators attacking preypopulations displaying Allee effects. By contrast, our work shows that the ability of generalist predatorsto control prey populations with logistic growth lies in the bistable dynamics of the coupled system.We also observe pseudo-Allee effects in our system, but their physics is quite different. An analysis ofthe ODE system identified parameter regions of monostable (extinction) and bistable (extinction orinvasion) dynamics, but analysis of the associated PDE was able to distinguish different and additionalregions of invasion and extinction. As a consequence, prey control was predicted to be possible whenspace was considered in additional situations other than those identified without considering space. Thereverse situation was also possible. None of these considerations apply to spatial predator-prey systemswith specialist natural enemies. 16
Proofs (2.1)
Let
E > α ≥ h ≥
0, the system (4) can be rewriten as (cid:26) ddt u = Θ h ( u )( f h ( u ) − v ) ddt v = rv ( g h ( u ) − v ) (7)where Θ h ( u ) = Eu Ehu , f h ( u ) = E (1 − u )(1 + Ehu ) and g h ( u ) = 1 + α Eu Ehu . (i) Proof of the existence. Define H ( h, u ) = f h ( u ) − g h ( u ). For a given h ≥
0, a couple ( u, v ) is a positivestationary solution of (7) if and only if v = f h ( u ) and u ∈ ]0 ,
1[ is a solution of H ( h, u ) = 0 . (8)Now, fix u ∈ ]0 , ∂ h H ( h, u ) = u (1 − u ) + α (cid:16) Eu Ehu (cid:17) >
0, one sees that the map h (cid:55)→ H ( h, u ) isincreasing. From E >
1, we get H (0 , u ) < u ∈ ]0 ,
1[ we get lim h → + ∞ H ( h, u ) = + ∞ . The map h → H ( h, u ) being continuous, this implies that for any u ∈ ]0 , h ( u ) > H ( h, u ) < h < h ( u ) ,H ( h, u ) = 0 if h = h ( u ) ,H ( h, u ) > h > h ( u ) . The smooth function u (cid:55)→ h ( u ) may be computed explicitly by noting that for any h ≥
0, the equation(8) is equivalent to the algebraic equation u ∈ ]0 ,
1[ is a solution of P h ( u ) = 0 (9)wherein we have set P h ( u ) = (1 − u )(1 + Ehu ) − E (1 + Ehu + Eαu ) . This yields the explicit formula h ( u ) = 1 Eu (1 − u ) (cid:20) E (cid:16) (cid:112) αu (1 − u ) (cid:17) + u − (cid:21) . (10)In particular lim u → + h ( u ) = lim u → − h ( u ) = + ∞ . (11)This implies that the minimum of h ( u ) is obtained for some u crit ∈ ]0 , h ∗ = inf u ∈ (0 , h ( u ) = min u ∈ (0 , h ( u ) = h ( u crit ) . (12)The definition of h ∗ shows that if h < h ∗ , then (9) has zero solution. This also implies, together withthe limits (11) and the continuity of u (cid:55)→ h ( u ), that for any fixed h > h ∗ the equation (9) admits atleast two solutions (see the figure 9). In addition to this, for any fixed h >
0, one has deg( P h ) = 3 and P h always admits a negative roots for P (0) = 1 − E < x →−∞ P h ( x ) = + ∞ . This implies that (9)admits at most two solutions . In conclusion, (9) has exactly two positive solutions if h > h ∗ and zeropositive solution if h < h ∗ . This ends the proof of ( i ).
4. Remark that h ∗ > E , because h ≤ E , implies that f (cid:48) h < ,
1) and (8) has no solution.5. These arguments also show that u (cid:55)→ h ( u ) is decreasing on (0 , u ∗ ) and increasing on ( u ∗ ,
1) ; for otherwise it ispossible to choose h > u ∈ (0 ,
1) such that h = h ( u ), which is equivalent to P h having at least 4 roots. See the figure 9. crit uh 0 10 h ∗ f h ( u ) > g h ( u ) f h ( u ) < g h ( u ) h > h ∗ g h ( u ) E uv f h ( u ) h = h ∗ g h ( u ) E uv f h ( u ) h < h ∗ g h ( u ) E uv f h ( u ) Figure
9: The four figures are computed for E = α = 2, which gives h ∗ ≈ .
36. The figure on the left representsthe curve u (cid:55)→ h ( u ) (in bold). Above the curve f h ( u ) > g h ( u ) while below the curve f h ( u ) < g h ( u ). Fora given h , the ordinate u of a point of this curve verifies f h ( u ) = g h ( u ) and corresponds to the positivestationary solution ( u, f h ( u )) of (7). The figures on the right represent the isoclines v = f h ( u ) and v = g h ( u ) for the three fixed values h = 3, h = 4 .
36 and h = 5. A positive stationary solution of thesystem (7) corresponds to an intersection of these two isoclines. The system (7) admits two positivesolutions for h > h ∗ , one (double) solution for the critical case h = h ∗ and zero positive solution for h < h ∗ . (ii) Proof of the stability . Let h > h ∗ be fixed. Let ( u, v ) be a positive stationary solution of (7). Since( u, v ) verifies f h ( u ) = g h ( u ) = v , the Jacobian matrix at ( u, v ) reads J ( u, v ) = (cid:20) Θ h ( u ) f (cid:48) h ( u ) − Θ h ( u ) rvg (cid:48) h ( u ) − rv (cid:21) hence det (cid:0) J ( u, v ) (cid:1) = r Θ h ( u ) v (cid:0) g (cid:48) h ( u ) − f (cid:48) h ( u ) (cid:1) . From the proof of (i), we know that the system (7) admits exactly two positive solutions denotedrespectively as ( (cid:98) u, (cid:98) v ) and ( u ∗ , v ∗ ) with (cid:98) u < u crit < u ∗ and such that h = h ( (cid:98) u ) = h ( u ∗ ) , (13)where u (cid:55)→ h ( u ) is given by (10). In particular, (cid:98) u and u ∗ are the solution of H ( h ( u ) , u ) = 0 . (14)Differentiating the equation (14) with respect to u gives ∂ u H ( h ( u ) , u ) = − h (cid:48) ( u ) ∂ h H ( h ( u ) , u ) . ∂ h H ( h ( u ) , u ) >
0, the identity (13) and the footnote 5, one gets ∂ u H ( h, (cid:98) u ) < ∂ u H ( h, u ∗ ) >
0. Since ∂ u H ( h, u ) = − ( g (cid:48) h ( u ) − f (cid:48) h ( u )), this shows that det (cid:0) J ( (cid:98) u, (cid:98) v ) (cid:1) < (cid:98) u, (cid:98) v ) follows.By contrast, one has det (cid:0) J ( u ∗ , v ∗ ) (cid:1) > u ∗ , v ∗ ) is given by the signof tr (cid:0) J ( u ∗ , v ∗ ) (cid:1) = Θ h ( u ∗ ) f (cid:48) h ( u ∗ ) − rv ∗ . (15)In order to highlight the dependence on h , for any h > h ∗ , we note u ∗ = u ∗ ( h ) and we also define µ ( h ) = Eh − Eh . From f (cid:48) h ( u ) = E h ( µ ( h ) − u ) and (15), we infer the following :– If µ ( h ) ≤ u ∗ ( h ) then ( u ∗ , v ∗ ) is asymptotically stable.– If µ ( h ) > u ∗ ( h ) then the stability of ( u ∗ , v ∗ ) depends on r . More precisely, define r crit = Θ h ( u ∗ ) f (cid:48) h ( u ∗ ) v ∗ > . (It is easy to show that r crit ≤ r > r crit , then ( u ∗ , v ∗ ) is asymptotically stable.– If r < r crit , then ( u ∗ , v ∗ ) is unstable.The sign of µ ( h ) − u ∗ ( h ) with respect to the parameter h remains to be found.On a first hand, the explicit expression of f (cid:48) h shows that f (cid:48) h ∗ is decreasing and that f (cid:48) h ∗ ( µ ( h ∗ )) = 0.Moreover, the definition (12) of h ∗ yields u ∗ ( h ∗ ) = u crit and f (cid:48) h ∗ ( u crit ) = g (cid:48) h ∗ ( u crit ) >
0. Hence, µ ( h ∗ ) E (see footnote 4), this implies that µ ( h ) = u ∗ ( h ) has exactly onesolution in ( h ∗ , + ∞ ). We note this unique solution h ∗∗ = h ∗∗ ( E, α ).Finally, it is clear from the above arguments that µ ( h ) < u ∗ ( h ) if h ∈ ( h ∗ , h ∗∗ ) and that µ ( h ) > u ∗ ( h )if h > h ∗∗ . This ends the proof of the Theorem. Proofs of properties 2.2
Similarly to the proof of the theorem 2.1, and to highlight the role of theparameters E and α , let us define H ( E, α, h, u ) = 1 E (1 − u )(1 + Ehu ) − (cid:18) α Eu Ehu (cid:19) . From the proof of the theorem 2.1, we know that the quantity h ∗ = h ∗ ( E, α ) and the corresponding u crit = u crit ( E, α ) are characterized by the two equations H ( E, α, h ∗ , u crit ) = 0 (16a) ∂ u H ( E, α, h ∗ , u crit ) = 0 . (16b)The implicit function theorem immediatly shows that the maps ( E, α ) (cid:55)→ ( h ∗ , u crit ) belongs to C (cid:0) (1 , + ∞ ) × [0 , + ∞ ); R (cid:1) .– Proofs of the growth of E (cid:55)→ h ∗ ( E, α ) and of α (cid:55)→ h ∗ ( E, α ) . Let α ≥ E and use (16b) yields ∂ E H ( E, α, h ∗ , u crit ) + ∂ E h ∗ ( E, α ) · ∂ h H ( E, α, h ∗ , u crit ) = 0 . We already know that ∂ h H < ∂ E H ( E, α, h ∗ , u crit ) = 1 E (1 + Eh ∗ u crit ) > . It follows that ∂ E h ∗ ( E, α ) >
0. Similar arguments show that ∂ α h ∗ ( E, α ) > Computation of lim E → h ∗ ( E, α ) . Let α ≥ h ∗ ( · , α ) is increasing and positive on (1 , + ∞ ), there exists a nonnegativescalar h ∗ ( α ) such that h ∗ ( E, α ) → h ∗ ( α ) as E →
1. To compute this limit, denote P ( E, α, h, u ) =(1 +
Ehu ) H ( E, α, h, u ). From (16) and the definition of h ∗ , one see that h ∗ is the minimal valueof h such that ∃ u ∈ [0 , , P ( E, α, h, u ) = ∂ u P ( E, α, h, u ) = 0 . (17)In other words, h ∗ is the minimal value of h such that P ( E, α, h, · ) admits a multiple root in [0 , E → h ∗ ( α ) is the minimal value of h such that ∃ u ∈ [0 , , P (1 , h, α, u ) = ∂ u P (1 , h, α, u ) = 0 . (18)Explicit computations give P (1 , α, h, u ) = u (cid:0) h u + h (2 − h ) u + α + 1 − h (cid:1) . The multiplicity of the roots of P (1 , α, h, · ) need now to be discussed. • If h < √ α then 0 is the only root of P (1 , α, h, · ) and the multiplicity is one. • If h ≥ √ α then P (1 , α, h, · ) has two other real roots and may be explicitly written as P (1 , α, h, u ) = h u ( u − u − )( u − u + )where u ± = 12 h (cid:16) h − ± (cid:112) (2 − h ) − α + 1 − h ) (cid:17) . – If h = 2 √ α then u − = u + = 1 − √ α and this multiple root belongs to [0 ,
1) if and only if α ≥ √ α < h < α + 1 then the three roots 0, u − and u + are distinct and there is no multiple root.– If h = α + 1 then either 0 = u − or 0 = u + , depending on the sign of α −
1. In both cases 0 is amultiple root.– Finally, if h > α + 1 then u − < < u + and all the roots of P (1 , α, h, · ) have multiplicity one.The above discussion shows, using the characterization of h ∗ ( α ), that h ∗ ( α ) = (cid:26) α if α ≤ √ α if α > . (19)Note that α (cid:55)→ h ∗ ( α ) belongs to C ([0 , + ∞ ) , R + ).– Computation of lim E → + ∞ h ∗ ( E, α ) . Since h ∗ ( · , α ) is positive and increasing, one haslim E → + ∞ h ∗ ( E, α ) = 1 (cid:96) α for some nonnegative number (cid:96) α (wherein we have set = + ∞ ).However, since u crit ( E, α ) is bounded, up to a subsequence again denoted by E , one has u crit ( E, α ) → µ for some positive number µ (eventually depending on α ). By passing to the limit in (16a), wededuce µ > H ( E, α, h ∗ , u crit ) → + ∞ . It follows that taking E → + ∞ in (16b), one obtains µ = 1 /
2. Thus, by passing to the limit in (16a), one gets1 (cid:96) α = 2 (cid:96) α + 2 α (cid:96) α being nonnegative) 1 (cid:96) α = 2 + 2 √ α ;– Computation of h ∗ ( E, . If α = 0 one gets P ( E, h, , u ) = 1 E (1 + Ehu ) (cid:0) − Ehu + u ( Eh −
1) + 1 − E (cid:1) . (20)Standard computations show that P ( E, h, , · ) admits two nonnegative roots if and only if h > E (2 E − (cid:112) E ( E −
1) := h ( E ) and one double nonnegative root if h = h ( E ). This shows that h ∗ ( E,
0) = h ( E ).– Computation of lim α → + ∞ h ∗ ( E, α ) . Since E (cid:55)→ h ∗ ( E, α ) is increasing one gets for any E ≥ h ∗ ( E, α ) ≥ h ∗ ( α ). The explicit expression (19) of h ∗ ( α ) implies lim α → + ∞ h ∗ ( E, α ) = + ∞ . Let
E > u ( t, x ) , v ( t, x )) be a solution of (2) with initial condition ( u ( x ) , v ( x ))verifying (3). Step 1.
It is clear that u ( t, x ) ≥ t > ∂ t v ( t, x ) − d ∆ x v ( t, x ) ≥ rv ( t, x )(1 − v ( t, x )) , t > , x ∈ R . Thus, using the comparison principle and the hypothesis that v (0 , x ) ≥
1, we deduce v ( t, x ) ≥ t ≥ x ∈ R . It follows that ∂ t u ( t, x ) − ∆ x u ( t, x ) ≤ u ( t, x )(1 − u ( t, x )) − Eu ( t, x )1 + Ehu ( t, x ) , t > , x ∈ R . By the comparison principle, we infer that any solution u of ∂ t u ( t, x ) − ∆ x u ( t, x ) = u ( t, x )(1 − u ( t, x )) − Eu ( t, x )1 + Ehu ( t, x ) , t > , x ∈ R (21)such that u (0 , x ) ≥ u ( x ) for any x ∈ R satisfies ∀ t > , ∀ x ∈ R , u ( t, x ) ≤ u ( t, x ) . In particular, let φ ( t ) be the solution of the ordinary differential equation ddt φ ( t ) = φ ( t )(1 − φ ( t )) − Eφ ( t )1 + Ehφ ( t ) , t > φ (0) = 1. φ is a homogeneous solution of (21) and from u (0 , x ) ≤ φ (0) we deduce ∀ t > , ∀ x ∈ R , ≤ u ( t, x ) ≤ φ ( t ) . tep 2. Behavior of φ ( t ) as t → + ∞ .It is clear that 0 is always a steady state of (22) and is asymptotically stable for E > φ > φ is a root of the polynomial P ( E, h, , · ) whichis studied in the proof of the property 2.2 (see (20)). Hence, if h ≥ h ( E ) then (22) has two positivesteady states that we denote as u − ( E, h ) ≤ u + ( E, h ) ≤ u ± ( E, h ) = 12 − Eh ± (cid:115)(cid:18) − Eh (cid:19) − E − Eh . (23)A linear analysis shows that for h > h ( E ), u − ( E, h ) is unstable and u + ( E, h ) is (asymptotically) stable.Finally, classical arguments show that φ ( t ) → u + ( E, h ) if h > h ( E ) and φ ( t ) → h < h ( E ). Thisends the proof of the theorem. Let
E > u ( t, x ) , v ( t, x )) be a solution to (2) with initial condition ( u ( x ) , v ( x )) verifying(3).By the argument of the step 1 of the proof of 5.2, one already knows that : u ( t, x ) ≤ u ( t, x ) (24)where u ( t, x ) is a solution of (21) with u (0 , x ) = u ( x ). Moreover, from the step 2 of the proof of 5.2,one knows that if h > h ( E ) then the equation (21) is bistable since (21) admits two stable nonnegativesteady states : u = 0 and µ = u + ( E, h ) <
1. We prove here that, in that case, there exists a trave-ling wave connecting µ to 0 at a negative speed if and only if h < h − ( E ) for some (implicit) number h − ( E ) > h ( E ). This implies that for any x ∈ R , u ( t, x ) → h < h − ( E ) and u ( t, x ) → µ if h > h − ( E ),which proves the theorem.Let h > h ( E ) be fixed. It can be proven (see Fife 1979, Volpert et al. 1994 and the referencetherein) that there exists a unique speed c = c ( E, h ) such that (21) admits a traveling solution of speed c which connects µ to 0. More precisely, there exists a profile U ( ξ ) verifying U ( −∞ ) = µ , U (+ ∞ ) = 0and U (cid:48) ( ±∞ ) = 0, such that U ( x − ct ) = u ( x, t ) is a solution of (21). Moreover (see Fife 1979) thistraveling wave describes the asymptotic behavior of all solutions provided the initial condition (3) areverified. The theorem is proven by showing that there exists h − ( E ) such that c ( E, h ) < h < h − ( E )and c ( E, h ) > h > h − ( E ). This result is a direct consequence of the two following lemma. The firstlemma gives a characterization of the sign of c by an explicit function (see figure 10). Lemma 5.1.
Define W ( E, h, u )) = (cid:90) u (cid:18) s (1 − s ) − Es Ehs (cid:19) ds.
Then sign ( c ( E, h )) = sign ( W ( E, h, µ )) where µ := u + ( E, h ) .Proof. Denote c = c ( E, h ) and let ξ = x − ct and U ( ξ ) = u ( x, t ) with U ( −∞ ) = µ , U (+ ∞ ) = 0 and U (cid:48) ( ±∞ ) = 0. We have − cU (cid:48) = U (cid:48)(cid:48) + U (1 − U ) − EU EhU = U (cid:48)(cid:48) + ∂W∂U ( E, h, U ) . U (cid:48) and integrating over R one gets − c (cid:90) + ∞−∞ ( U (cid:48) ) dZ = (cid:90) + ∞−∞ (cid:18) U (1 − U ) − EU EhU (cid:19) U (cid:48) dZ = (cid:90) µ ∂W∂U ( E, h, U ) dU = − W ( E, h, u + )and sign ( c ) = sign ( W ( E, h, u + ( E, h )) follows. u u − u + − W ( u + ) > c < u u − u + − W ( u + ) < c > Figure
10: Graph of u (cid:55)→ − W ( E, h, u ). The steady states 0, u − and u + correspond to critical points of thepotential − W . A stable steady state is a local minimum of W and an unstable steady state is a localmaximum. One sees that 0 and u + are two stable steady states. The sign of c characterizes which ofthem is the final global attractor : if c <
0, then 0 is the global attractor. If c >
0, then u + is theglobal attractor. The second lemma gives the sign of W with respect to h . Lemma 5.2.
For any
E > , there exists (a unique) h − ( E ) > h ( E ) such that sign ( W ( E, h, u + ( E, h ))) = sign (cid:0) h − h − ( E ) (cid:1) . Lemma 5.1 and 5.2 show that– There is invasion of prey for (21) ( c ( E, h ) >
0) if h ∈ ( h − ( E ) , + ∞ ). In that case, for any x ∈ R , u ( t, x ) → u + ( E, h ) as t → + ∞ .– There is extinction of prey for (21) ( c ( E, h ) <
0) if h ∈ [ h ( E ) , h − ( E )). In that case, for any x ∈ R , u ( t, x ) → t → + ∞ .In particular, we infer from the inequality (24), that if h ∈ [ h ( E ) , h − ( E )), then for any x ∈ R , u ( t, x ) → t → + ∞ . This ends the proof of the theorem.It remains to prove the lemma 5.2. Proof of lemma 5.2.
Define W ( E, h ) = W ( E, h, u + ( E, h )). From the lemma 5.1, we know that sign ( c ) = sign ( W ( E, h )). We show here that there exists h − ( E ) > h ( E ) such that sign ( W ( E, h )) = sign ( h − h − ( E )). Step 1.
Differentiate W with respect to h gives ∂ h W ( E, h ) = ∂ h W ( E, h, u + ( E, h )) + ∂ u W ( E, h, u + ( E, h )) ∂ h u + ( E, h ) . By the very definition of W and u + , one has ∂ u W ( E, h, u + ( E, h )) = 0, which yields ∂ h W ( E, h ) = ∂ h W ( E, h, u + ( E, h )) . W ( E, h, u ) = u − u − uh + 1 Eh ln(1 + Ehu ) . (25)Differentiate this expression with respect to h and denoting z = Ehu + ( E, h ) provides ∂ h W ( E, h, u + ( E, h )) = 1 Eh (cid:18) z − z ) + z z (cid:19) . A standard analysis shows that the map z (cid:55)→ z − z ) + z z takes positive values for z >
0. Itfollows that the map h (cid:55)→ W ( E, h ) is increasing.
Step 2.
Recalling that u + is the largest roots of (20), one verifies that u + ( E, h ) → h → + ∞ and W ( E, + ∞ ) = 1 / − / / > Step 3.
This step consists in proving that for any
E > W ( E, h ( E )) = W ( E, h ( E ) , u + ( E, h ( E ))) := g ( E ) is negative.From the explicit expression (23) of u + ( E, h ) and the definition of h ( E ), we get u + ( E, h ( E )) = 12 (cid:18) − Eh ( E ) (cid:19) = 1 − E + (cid:112) E ( E − . On a first hand, we have g (cid:48) ( E ) = ∂ E W ( E, h ( E ) , u + ( E, h ( E ))) + ∂ h W ( E, h ( E ) , u + ( E, h ( E ))) · h (cid:48) ( E ) . Straightforward calculations give g (cid:48) ( E ) = (cid:18) Eh ( E ) (cid:19) (cid:20) − ln(1 + z ) (cid:18) z (cid:19)(cid:21) . wherein we have set z = Eh ( E ) u + ( E, h ( E )) = √ E − (cid:16) √ E + √ E − (cid:17) = 12 ( Eh ( E ) − . A standard analysis shows that 2 zz +2 < ln(1 + z ) for any z >
0. Thus g (cid:48) ( E ) < E > h (1) = 1 and u + (1 ,
1) = 0, so that g (1) = 0. It follows that g ( E ) < E > Conclusion
For any
E >
1, the map h (cid:55)→ W ( E, h ) is increasing and verifies W ( E, h ( E )) < h → + ∞ W ( E, h ) = 1 /
6. By continuity, for any
E >
1, there exists a unique h − = h − ( E ) ∈ ( h ( E ) , + ∞ )such that sign ( W ( E, h, u + )) = sign ( h − h − ( E )). This ends the proof of lemma 5.2. Proofs of properties 2.5.
Let
E >
Growth of h − ( · ) . Recall that h − ( E ) is characterized by W ( E, h − ( E ) , u + ( E, h − ( E ))) = 0 (26)where the expressions of W and u + are respectively given in (25) and (23). By definition of W ,one has ∂ u W (cid:0) E, h − ( E ) , u + ( E, h − ( E )) (cid:1) = 0 . ∂ E W (cid:0) E, h − ( E ) , u + ( E, h − ( E )) (cid:1) + dh − dE ( E ) · ∂ h W (cid:0) E, h − ( E ) , u + ( E, h − ( E )) (cid:1) = 0 . The explicit computation of ∂ E W and ∂ h W are done in proof 5.3 and one gets ∂ E W (cid:0) E, h − ( E ) , u + ( E, h − ( E )) (cid:1) < ∂ h W (cid:0) E, h − ( E ) , u + ( E, h − ( E )) (cid:1) > dh − dE ( E ) = − ∂ E W ( E, h − ( E ) , u + ( E, h − ( E ))) ∂ h W ( E, h − ( E ) , u + ( E, h − ( E ))) > Limit of h − ( E ) as E → . Let
E >
1. One knows that h − ( · ) is increasing on (1 , + ∞ ) and minored by h ( E ) >
0. Hence h − ( E ) admits a limit (cid:96) − as E →
1, and from h (1) = 1, one obtains (cid:96) − ≥ . The explicit expressions (23) of u ± yield u − (1 , (cid:96) − ) = 0 and u + (1 , (cid:96) − ) = 1 − (cid:96) − (cid:96) − . Assume by contradiction that u + (1 , (cid:96) − ) >
0. Then, for any h ≥ < u < u + ( h, ∂ u W (1 , h, u ) = u (1 − u ) −
11 + hu > W (1 , h,
0) = 0 . This implies W (1 , (cid:96) − , u + (1 , (cid:96) − )) <
0, which contradicts (26).It follows that u + (1 , (cid:96) − ) = 0. Thus (cid:96) − = 1, which reads lim E → h − ( E ) = 1 . – Limit of h − ( E ) as E → + ∞ . One knowns that h − ( E ) is increasing so that there is some (cid:96) ≥ E → + ∞ h − ( E ) = 1 (cid:96) (wherein we have set ∞ = 0). By taking the limit E → + ∞ in (23), we obtainlim E → + ∞ u + ( E, h − ( E )) = u (cid:96) := 12 (1 + √ − (cid:96) ) > . and then, by taking the limit in (26), u (cid:96) − u (cid:96) − u (cid:96) (cid:96) = 0 . Thus, (cid:96) = which reads, lim E → + ∞ h − ( E ) = . As in the proof 5.3, we start here by showing that 0 ≤ u ( t, x ) ≤ u ( t, x ) for some function u verifyinga scalar reaction-diffusion equation (27) depending on E , h and α . Next, we show that there exists h + ( E, α ) such that lim t → + ∞ u ( t, x ) = µ > h > h + ( E, α ). This step is done using the proof ofTheorem 2.4 together with an appropriate change of variables.25 tep 1.
Let
E > ≤ u ≤
1, one has 1 ≤ v ≤ v with v := v ( E, h, α ) = 1 + α E
Eh .
From the estimate v ≤ v, we get u ( t, x ) ≤ u ( t, x ) where u verifies ∂ t u = ∆ x u + u (1 − u ) − Evu
Ehu , t > , x ∈ R . (27)Denoting (cid:101) E = Ev and (cid:101) h = hv , this equation reads simply ∂ t u = ∆ x u + u (1 − u ) − (cid:101) Eu (cid:101) E (cid:101) hu , t > , x ∈ R . (28)Note that, removing the (cid:101) , this equation is nothing but (21). As a consequence, the proof of the theorem2.4 implies the following result on (28). Lemma 5.3.
Let (cid:101) u ( t, x ) be a solution of (28) verifying the initial condition (3) . There exists h − ( (cid:101) E ) >h ( (cid:101) E ) such that– if (cid:101) h < h − ( (cid:101) E ) then for any x ∈ R , lim t → + ∞ (cid:101) u ( t, x ) = 0 (extinction),– if (cid:101) h > h − ( (cid:101) E ) then for any x ∈ R , lim t → + ∞ (cid:101) u ( t, x ) = µ ( (cid:101) E, (cid:101) h ) > (invasion). Step 2.
Recalling that (cid:101) E = Ev ( E, h, α ) and (cid:101) h = hv ( E,h,α ) do depend on h , we set µ = µ ( E, h ) = µ ( (cid:101) E, (cid:101) h )and the previous lemma gives the following implicit condition on h for invasion to occur.For any x ∈ R , lim t → + ∞ u ( t, x ) = µ , provided : h > v ( E, h, α ) h − ( Ev ( E, h, α )) . (29)The following lemma gives an equivalent condition for this implicit condition to occur. This ends theproof of theorem 2.4. Lemma 5.4.
For any
E > and α ≥ , there exists (a unique) h + ( E, α ) ≥ h − ( E ) such that : (29) holds true if and only if h > h + ( E, α ) . Proof of lemma 5.4 :
Let
E > α ≥ F E,α ( h ) = h − v ( E, h, α ) h − ( Ev ( E, h, α )) . (30)By the construction of h − via the implicit function theorem, one knows that F E,α ( · ) is a C function.Moreover, since h − and v are bounded, we have F E,α (+ ∞ ) = + ∞ and F E,α (0) = − v ( E, , α ) h − ( Ev ( E, , α )) = − (1 + αE ) h − ( E (1 + αE )) < . Therefore, it suffices to show that F E,α is increasing.One has F (cid:48) E,α ( h ) = 1 − ∂ h v ( E, h, α ) (cid:18) h − ( Ev ) + v dh − dE ( Ev ) (cid:19) . From the expression of v , we infer ∂ h v < dh − dE >
0. It follows that F (cid:48) E,α ( h ) >
0. This ends the proof of the lemma.26 roofs of properties 2.7. – Proof of the growth of h + ( E, · ) and h + ( · , α ) . h + ( E, α ) is characterized by an equality in (29), that is, F E,α ( h + ( E, α )) = 0 (31)where F E,α is defined in (30). Differentiating (31) with respect to E gives ∂ E h + ( E, α ) · F (cid:48) E,α ( h + ( E, α )) + ∂ E F E,α ( h + ( E, α )) = 0and then ∂ E h + ( E, α ) · F (cid:48) E,α ( h + ( E, α )) = (cid:2) ( Ev + h − ) ∂ E v + v (cid:3) · dh − dE ( Ev ) . wherein we have set h − = h − ( Ev ) , v = v ( E, h + ( E, α ) , α ) and ∂ E v = ∂ E v ( E, h + ( E, α ) , α ) . One already knows that F (cid:48) E,α ( h + ( E, α )) > dh − dE >
0. A direct computation shows that ∂ E v >
0, which leads to ∂ E h + ( E, α ) > α gives, with obvious notations, ∂ α h + ( E, α ) · F (cid:48) E,α ( h + ( E, α )) = ∂ α v · (cid:18) h − + E dh − dE (cid:19) , and since ∂ α v >
0, one obtains ∂ α h + ( E, α ) > Limits of h + ( E, α ) as α → and α → + ∞ . From v ( E, h,
0) = 1 we deduce h + ( E,
0) = h − ( E ). Now, it is clear from the construction of h + ,that h + ( E, α ) ≥ h ∗ ( E, α ). From the properties 2.2, we obtain h + ( E, α ) → + ∞ as α → + ∞ . – Limit of h + ( E, α ) as E → + ∞ . First, recall that h + ( E, α ) is characterized by (31), which reads h + ( E, α ) = (cid:18) α E Eh + ( E, α ) (cid:19) · h − (cid:18) E + α E Eh + ( E, α ) (cid:19) . (32)Since h + ( · , α ) is increasing, there exists (cid:96) α ≥ E → + ∞ h + ( E, α ) = (cid:96) α . Taking thelimit E → + ∞ in (32) and using h − ( E ) → , one obtains (cid:96) α = · (1 + α(cid:96) α ). The resolution ofthis equation ends the proof. Acknowledgements
The authors would like to thank the two anonymous reviewers for their valuablecomments and suggestions to improve the quality of this manuscript.
R´ef´erences
Alzahrani EO, Davidson FA and Dodds N (2012) Reversing invasion in bistable systems. Journal of MathematicalBiology, 65 :1101–1124.Basnet K and Mukhopadhyay A (2014) Biocontrol potential of the lynx spider
Oxyopes javanus (Araneae :Oxyopi-dae) against the tea mosquito bug,
Helopeltis theivora (Heteroptera :Miridae). International Journal of TropicalInsect Science 34(4) :232-238.Braverman E, Kamrujjaman M and Korobenko L (2015) Competitive spatially distributed population dynamicsmodels : Does diversity in diffusion strategies promote coexistence ? Mathematical biosciences, 264 :63-73. oukal DS, Sabelis MW and Berec L (2007) How predator functional responses and Allee effects in prey affectthe paradox of enrichment and population collapses. Theoretical Population Biology, 72 :136-147.Cai Y, Banerjee M, Kang Y and Wang W (2014) Spatiotemporal Complexity in a predator-prey model with weakallee Effects. Mathematical Biosciences and Engineering, 11 :1247-1274Chakraborty S (2015) The influence of generalist predators in spatially extended predator–prey systems. EcologicalComplexity, 23 :50-60.De Roos AM, Mccauley E and Wilson WG (1991) Mobility versus density-limited predator–prey dynamics ondifferent spatial scales. Proceedings of the Royal Society of London B : Biological Sciences, 246(1316) :117-122.De Roos AM, Mccauley E and Wilson WG (1998) Pattern formation and the spatial scale of interaction betweenpredators and their prey. Theoretical Population Biology, 53(2) :108-130.Du Y and Shi J (2007) Allee effect and bistability in a spatially heterogeneous predator-prey model. Transactionsof the American Mathematical Society, 359(9) :4557-4593.Ducrot A and Langlais M (2012) A singular reaction-diffusion system modelling prey-predator interactions :Invasion and co-extinction waves. Journal of Differential Equations, 253 :502-532.Dunbar SR (1984) Traveling Wave Solutions of Diffusive Lotka-Volterra Equations. Journal of MathematicalBiology, 17 :11-32.Dunbar SR (1984) Traveling Wave Solutions of Diffusive Lotka-Volterra Equations : A heteroclinic connection inR4. Transaction of the American Mathematical Society, 286 :557-594.Erbach A, Lutscher F and Seo G (2014) Bistability and limit cycles in generalist predator–prey dynamics. Eco-logical Complexity, 14 :48–55.Fagan WF, Lewis MA, Neurbert MG, Driessche Pvd (2002), Invasion theory and biological control. EcologyLetters, 5 :148-157.Fife PC (1979), Long Time Behavior of Solutions of Bistable Nonlinear Diffusion Equations. Archive for RationalMechanics and Analysis, 70 :31–46.Fischer RA (1937), The wave of advance of an advantageous gene. Annals of Eugenics, 7 :353-369.Garnier J, Roques L and Hamel F (2012), Success rate of a biological invasion in terms of the spatial distributionof the founding population. Bulletin of Mathematical Biology, 74(2) :453-473.Giovangigli V (1990), Nonadiabatic Plane Laminar Flames and Their Singular Limits SIAM J. Math. Anal.,21(5) :1305-1325.Hastings A (2000), Parasitoid spread : lessons for and from invasion biology. Parasitoids Population Biology, (M.E. Hochberg and A. R. Ives eds). Princeton, NJ : Princeton University Press, 70-82.Hoyle A and R.G. Bowers RG (2007), When is evolutionary branching in predator–prey systems possible with anexplicit carrying capacity ? Mathematical Biosciences, 210 :1–16.Huang J, Lu G and Ruan S (2003), Existence of traveling wave solutions in a diffusive predator-prey model.Journal of Mathematical Biology, 46 :132-152.Huang Y and Weng P (2013), Traveling waves for a diffusive predator-prey system with a general functionalresponse. Nonlinear Analysis : Real World Applications, 14 :940-959.Krivan V and Eisner J (2006), The effect of the Holling type II functional response on apparent competition.Theoretical Population Biology, 70 :421–430.Kolmogorov AN, Petrowskii I and Piscounov N (1937), Etude de l’´equation de la diffusion avec croissance de laquantit´e de mati´ere et son application `a un probl`eme biologique. Moscow University Mathematics Bulletin,1 :1-25. am KY and Ni WM (2012), Uniqueness and complete dynamics in heterogeneous competition-diffusion systems.SIAM Journal on Applied Mathematics, 72 :1695-1712.van Leeuwen E, Jansen VAA and Bright PW (2007), How population dynamics shape the functional response ina one-predator–two-prey system. Ecology, 88(6) :1571-1581.van Leeuwen R, Brannsrom A, Jansen VAA, Dieckmann U and Rossberg AG (2013), Journal of TheoreticalBiology, 328 :89-98.Li WT and Wu SL (2008), Traveling waves in a diffusive predator–prey model with holling type-III functionalresponse. Chaos, Solitons and Fractals, 37 :476-486.Lewis MA, Maini PK and Petrovskii SV (2013), Dispersal, Individual Movement and Spatial Ecology, Berlin,Germany : SpringerLi DS, Liao C, Zhang BX and Song ZW (2014), Biological control of insect pests in litchi orchards in China.Biological Control, 68 :23-36.Magal C, Cosner C, Ruan S and Casas J (2008), Control of invasive hosts by generalist parasitoids. MathematicalMedicine and Biology, 25 :1-20.Mistro DC, Rodrigues LAD and Petrovskii S (2012), Spatiotemporal complexity of biological invasion in a space-and time-discrete predator-prey system with the strong Allee effect. Ecological Complexity, 9 :16-32.Moffat CE, Lalonde RG, Ensing DJ, De Clerck-Floate RA, Grosskopf-Lachat G and Pither J (2013), Frequency-dependent host species use by a candidate biological control insect within its native European range. BiologicalControl, 67 :498-508.Morozov A and Petrovskii S (2009), Excitable population dynamics, biological control failure, and spatiotemporalpattern formation in a model ecosystem. Bulletin of Mathematical Biology,71, 863-887.Morozov A and Petrovskii S (2013), Feeding on Multiple Sources : Towards a Universal Parameterizationof the Functional Response of a Generalist Predator Allowing for Switching. PLoS ONE 8(9) : e74586.doi :10.1371/journal.pone.0074586Owen MR, Lewis MA (2001), How Predation can Slow, Stop or Reverse a Prey Invasion. Bulletin of MathematicalBiology, 63 : 655-684.Pimentel D (2011), Biological Invasions : Economic and Environmental Costs of Alien Plant, Animal, and MicrobeSpecies. Second Edition. CRC Press, New York.Potapov A and Rajakaruna H (2013), Allee threshold and stochasticity in biological invasions : Colonization timeat low propagule pressure. Journal of Theoretical Biology, 337 :1-14.Savage D and Renton M (2013), Requirements, design and implementation of a general model of biological invasion.Ecological Model, 272 :394-409.Smoller J (1983), Shock waves and reaction-diffusion equations. Springer-Verlag, New York.Turchin P (2003), Complex Population Dynamics : A Theoretical/Empirical Synthesis, Monographs in PopulationBiology. Princeton University Press, Princeton, NJ.Volpert AI, Volpert VA and Volpert VA (1994), Traveling Wave Solutions of Parabolic Systems. MathematicalMonographs 140.Wang W, Feng X and Chen X (2013), Biological Invasion and Coexistence in Intraguild Predation. Journal ofApplied Mathematics, 2013 :12pp.Z. Xu Z and P. Weng P (2012), Traveling Waves in a Diffusive Predator-Prey Model with General FunctionalResponse. Electronic Journal of Differential Equations, 197 :1-13.Ye GY, Xiao Q, Chen M, Chen X, Yuan Z, Stanley DW and Hu C (2014), Tea : Biological control of insect andmite pests in China. Biological Control, 68 :73-91.am KY and Ni WM (2012), Uniqueness and complete dynamics in heterogeneous competition-diffusion systems.SIAM Journal on Applied Mathematics, 72 :1695-1712.van Leeuwen E, Jansen VAA and Bright PW (2007), How population dynamics shape the functional response ina one-predator–two-prey system. Ecology, 88(6) :1571-1581.van Leeuwen R, Brannsrom A, Jansen VAA, Dieckmann U and Rossberg AG (2013), Journal of TheoreticalBiology, 328 :89-98.Li WT and Wu SL (2008), Traveling waves in a diffusive predator–prey model with holling type-III functionalresponse. Chaos, Solitons and Fractals, 37 :476-486.Lewis MA, Maini PK and Petrovskii SV (2013), Dispersal, Individual Movement and Spatial Ecology, Berlin,Germany : SpringerLi DS, Liao C, Zhang BX and Song ZW (2014), Biological control of insect pests in litchi orchards in China.Biological Control, 68 :23-36.Magal C, Cosner C, Ruan S and Casas J (2008), Control of invasive hosts by generalist parasitoids. MathematicalMedicine and Biology, 25 :1-20.Mistro DC, Rodrigues LAD and Petrovskii S (2012), Spatiotemporal complexity of biological invasion in a space-and time-discrete predator-prey system with the strong Allee effect. Ecological Complexity, 9 :16-32.Moffat CE, Lalonde RG, Ensing DJ, De Clerck-Floate RA, Grosskopf-Lachat G and Pither J (2013), Frequency-dependent host species use by a candidate biological control insect within its native European range. BiologicalControl, 67 :498-508.Morozov A and Petrovskii S (2009), Excitable population dynamics, biological control failure, and spatiotemporalpattern formation in a model ecosystem. Bulletin of Mathematical Biology,71, 863-887.Morozov A and Petrovskii S (2013), Feeding on Multiple Sources : Towards a Universal Parameterizationof the Functional Response of a Generalist Predator Allowing for Switching. PLoS ONE 8(9) : e74586.doi :10.1371/journal.pone.0074586Owen MR, Lewis MA (2001), How Predation can Slow, Stop or Reverse a Prey Invasion. Bulletin of MathematicalBiology, 63 : 655-684.Pimentel D (2011), Biological Invasions : Economic and Environmental Costs of Alien Plant, Animal, and MicrobeSpecies. Second Edition. CRC Press, New York.Potapov A and Rajakaruna H (2013), Allee threshold and stochasticity in biological invasions : Colonization timeat low propagule pressure. Journal of Theoretical Biology, 337 :1-14.Savage D and Renton M (2013), Requirements, design and implementation of a general model of biological invasion.Ecological Model, 272 :394-409.Smoller J (1983), Shock waves and reaction-diffusion equations. Springer-Verlag, New York.Turchin P (2003), Complex Population Dynamics : A Theoretical/Empirical Synthesis, Monographs in PopulationBiology. Princeton University Press, Princeton, NJ.Volpert AI, Volpert VA and Volpert VA (1994), Traveling Wave Solutions of Parabolic Systems. MathematicalMonographs 140.Wang W, Feng X and Chen X (2013), Biological Invasion and Coexistence in Intraguild Predation. Journal ofApplied Mathematics, 2013 :12pp.Z. Xu Z and P. Weng P (2012), Traveling Waves in a Diffusive Predator-Prey Model with General FunctionalResponse. Electronic Journal of Differential Equations, 197 :1-13.Ye GY, Xiao Q, Chen M, Chen X, Yuan Z, Stanley DW and Hu C (2014), Tea : Biological control of insect andmite pests in China. Biological Control, 68 :73-91.