Multiple wave solutions in a diffusive predator-prey model with strong Allee effect on prey and ratio-dependent functional response
Edgardo Villar-Sepúlveda, Pablo Aguirre, Víctor F. Breña-Medina
MMULTIPLE WAVE SOLUTIONS IN A DIFFUSIVE PREDATOR-PREY MODELWITH STRONG ALLEE EFFECT ON PREY AND RATIO-DEPENDENTFUNCTIONAL RESPONSE
EDGARDO VILLAR–SEP ´ULVEDA (cid:63) , PABLO AGUIRRE (cid:63) , AND V´ICTOR F. BRE ˜NA–MEDINA † Abstract.
A thorough analysis is performed in a predator-prey reaction-diffusion model which in-cludes three relevant complex dynamical ingredients: (a) a strong Allee effect; (b) ratio-dependentfunctional responses; and (c) transport attributes given by a diffusion process. As is well-known inthe specialized literature, these aspects capture adverse survival conditions for the prey, predationsearch features and non-homogeneous spatial dynamical distribution of both populations. We lookfor traveling-wave solutions and provide rigorous results coming from a standard local analysis, nu-merical bifurcation analysis, and relevant computations of invariant manifolds to exhibit homoclinicand heteroclinic connections and periodic orbits in the associated dynamical system in R . In sodoing, we present and describe a diverse zoo of traveling wave solutions; and we relate their occur-rence to the Allee effect, the spreading rates and propagation speed. In addition, homoclinic chaos ismanifested via both saddle-focus and focus-focus bifurcations as well as a Belyakov point. An actualcomputation of global invariant manifolds near a focus-focus homoclinic bifurcation is also presentedto enravel a multiplicity of wave solutions in the model. A deep understanding of such ecologicaldynamics is therefore highlighted. Introduction
One of the main objectives of ecology is to understand the interactions of individual organisms withthose of the same species and with other actors in their environment, and how populations distributewhen forming communities. In particular, predation —and its associated mathematical models— hasbeen part of the dominant topics in ecology from the beginning of this discipline [12, 51]. Sincethe early works of Volterra [62] onwards, predator-prey models have been continuously studied bymathematicians and biologists thanks to their adequate generalization of practical scenarios and forbeing an essential guide in biological populations. In all these cases, predator-prey models may shedlight on necessary conditions for species to coexist and not to be doomed to extinction [12, 51].In this work, we are interested in the interaction between predator and prey when this is driven byspatial distribution and propagation features —related to a transport process that favors food searchingevents—, and the critical densities that the prey population must achieve in order to survive in thelong term. The starting point of our work is the following Gause type deterministic predator-preymodel first studied in [4]: dNdt = rN (cid:18) − NK (cid:19) (cid:18) NM − (cid:19) − aN PP + ahN ; dPdt = e aN PP + ahN − qP. (1) Mathematics Subject Classification.
Key words and phrases.
Homoclinic and heteroclinic orbits, traveling waves, invariant manifolds, bifurcation analysis,Allee effect.EVS was partially funded by Programa de Incentivo a la Iniciaci´on Cient´ıfica PIIC DGIIP-UTFSM. EVS and PAwere partially funded by Proyecto Interno UTFSM PI-LI-19-06. PA also thanks Proyecto Basal CMM Universidad deChile. VFBM thanks the financial support by Asociaci´on Mexicana de Cultura A.C.. a r X i v : . [ n li n . C D ] A ug E. Villar–Sep´ulveda, P. Aguirre and V. Bre˜na–Medina
Here N = N ( t ) and P = P ( t ) are functions of time that represent the population sizes of prey andpredators, respectively, provided that P + ahN >
0. The main ingredients of (1) are a strong Alleeeffect on prey and a ratio-dependent functional response. Let us give a brief overview of these twofeatures. The Allee effect [56] is a phenomenon that may affect a species at low population densities,sometimes even related to the emergence of extinction risks. At low densities, predators can inducea phenomenon called the
Allee effect by predation affecting the prey species [28, 30]. Under theseunfavorable conditions, some species can have difficulties defending themselves from predators, hidingor taking anti-predation measures when their population density is low [64, 66]. Another manifestationof the Allee effect is the difficulty of a species to use cooperation and aggregation tactics in otherstrategic activities such as food searching, to resist unfavorable environmental conditions, or to findmating partners. Hence, while at a low population density, each individual may take advantage of ahigher amount of resources, sometimes this benefit is not as significant as the losses associated with itssmall density. The study of the Allee effect can have a massive impact on the dynamics of populationsof many plants and animals [20]. According to the literature, a species can be affected by two maintypes of Allee effect. The strong Allee effect is characterized by the presence of a threshold (known asthe
Allee threshold ). Below this threshold, the population density growth becomes negative, so thatthe population becomes prone to extinction [20]. In model (1) we assume that the prey is affectedby a strong Allee effect. This is manifested in the presence of parameter M . Indeed, in the absenceof predators, if the prey population is below the threshold M , we have dN/dt < < M (cid:28) K , where K is the carrying capacity of the prey.On the other hand, the so-called weak Allee effect describes a scenario in which low-level populationsshow a reduced (but positive) per capita growth rate [20]. Hence, a subtle relevant distinction betweenthese two density-dependent growth rates lies on the fact that the strong Allee effect captures whethera population critical size is reached. In other words, if the prey population size is exactly the Allee-threshold parameter value, the prey growth will only be dependent on interaction, as in (1) for instance,with the predators.Since the Allee effect is related to population extinction, an understanding of this phenomenon isvery relevant for ecological administration and conservation [56]. Many studies have been carried out inthis regard, as can be seen in [4, 18, 21, 56]. Different models have shown that the Allee effect increasesthe risk of population extinction as the Allee threshold rises [56]. When predator-prey interaction ispresent —as is the case in (1)—, this imposes a two-dimensional state space for the variables x and y . In such case, the Allee threshold is typically a one-dimensional object in the ( x, y )-plane. Recentfindings reveal that this threshold can be the basin boundary of a positive (coexistence) equilibrium,or a limit cycle, or a homoclinic orbit [4, 18].The second feature in (1) is a ratio-dependent functional response. The motivation for this comesfrom two well-known paradoxes in Gause type predator-prey models. First, enriching a predator-prey system causes an increase in the equilibrium density of the predator but not in that of theprey, destabilizing the community equilibrium; this is known as the paradox of enrichment [57, 42].Also, in Gause type models there can not exist a stable low prey density equilibrium, also knownas the biological control paradox [14, 48]. One way to deal with these paradoxes is by consideringratio-dependent functional responses [42]. A generic ratio-dependent functional response depends onthe population density ratio between prey and predators N/P . Some examples of ratio-dependentfunctional responses can be formulated by the direct substitution N (cid:55)→ N/P in the classic Hollingfamily of functional responses [22, 36, 37, 61]. In particular, in (1), prey consumption by predatorsfollows a modified Holling II ratio-dependent functional response. Here, e > h > a > ultiple wave solutions in a diffusive predator-prey model 3 consumption rate of predators. The remaining parameters in (1) are the prey population naturalintrinsic growth rate r >
0, and the predator natural mortality rate q >
N, P ) = (0 , u, v, τ ) := (cid:18) NK , PahK , rK (1 + ah ) M ( P + ahN ) t (cid:19) , ( α, β, γ, m ) := (cid:18) aMrK , eMrhK , qMrK , MK (cid:19) . In this way we obtain dudt = u ( u − m )(1 − u )( u + v ) − αuv,dvdt = βuv − γv ( u + v ) , (2)where u = u ( t ) and v = v ( t ) represent rescaled population sizes of prey and predators, respectively, andwe have renamed τ (cid:55)→ t once again for notation convenience. System (2) represents a C ∞ -qualitativelyequivalent polynomial extension of (1): it is C ∞ -qualitatively equivalent to (1) whenever ( P, N ) (cid:54) =(0 ,
0) and, it is well-defined at the origin. In [4], the authors sought to explain the consequencesof the (dis)appearance of limit cycles, homoclinic orbits, and heteroclinic connections in the globalarrangement of the phase plane near a Bogdanov-Takens bifurcation. In particular, they point out thatthe Allee threshold in the two-dimensional system is given as the boundary of the basin of attractionof an attracting positive equilibrium, and they determine conditions for mutual extinction or survivalof the populations.On the other hand, the spatial component of ecological interactions has been identified as a signifi-cant factor in the behavior of communities [31, 50, 51]. For instance, breaking symmetry features areusually found in ratio-dependent models that take into account certain heterogeneous spatial popula-tion distributions, (see e.g. [9]). We consider that even though individuals may follow a simple randommotion, the population however may behave as a whole in a deterministic way. Such a property yieldsto a transport process that can be understood from a mean-field viewpoint; that is, the populationtransport dynamics is captured by the Fokker-Planck equation, which in the simplest version particu-larly accounts for a standard diffusion process; see [52, 53]. Recent studies have considered differentspatial phenomena in population dynamics, such as propagation and migration [17, 44]. These haveled to obtaining a solid argument that highlights the relevance of their study by finding traveling wavesolutions that move in space, preserving their shape as time passes [46, 47, 59]. This kind of solutionsemerges as a suitable mathematical approach to describe wave-like spatial movement of populations.Evidence of this type of behavior has been found in several models used to model pest outbreaks,chemical concentration, colonization of a space region by a population, spatial dispersal of epidemics,etc.; see, for instance, the textbooks [50, 51] and the references therein. Typically, traveling wavesrepresent spatiotemporal transitions from one homogeneous steady state to another one, or to itself;see [38, 46, 51, 65]. Moreover, recent evidence indicates that multi-year periodic traveling waves canalso be found [59]. Following this motivation, we consider an extension of (2) by taking into account adiffusion transport process on both populations. In so doing, we obtain the following reaction-diffusion
E. Villar–Sep´ulveda, P. Aguirre and V. Bre˜na–Medina system: (cid:40) u t = D u xx + u ( u − m )(1 − u )( u + v ) − αuv ,v t = D v xx + βuv − γv ( u + v ) , (3)where partial derivatives short notation is used: w t = ∂w/∂t and w xx = ∂ w/∂x . In (3), x ∈ [ − L, L ] isthe spatial variable and the characteristic length of the predator-prey interaction domain is assumed tobe such that L → ∞ as this sort of solutions are known to arise in systems as (3) (see e.g. [15, 49, 45]); t > D , D > Auto [25] (with its extension
HomCont [16]) and
Matcont [23] with high accuracy. Inaddition, one can locate homoclinic and heteroclinic connections as intersections of global invariantmanifolds. This can be achieved by direct computation and inspection of the manifolds [2, 7] or bysetting additional techniques such as Lin’s method [39]. While one-dimensional invariant manifoldscan be approximated using straightforward integration from a given initial condition, the computationof higher-dimensional manifolds of equilibria and periodic orbits requires advanced numerical tech-niques; see for example [1, 32, 40, 41]. Moreover, while some works have dealt with traveling wavesassociated with three-dimensional vector fields in population models [47, 60, 67], our problem involvesa four-dimensional phase space, which is as a major challenge. Although the human brain is effi-cient when capturing depth in flat images of 3-dimensional objects, this ability is not as effective inhigher dimensions [19, 54]. When trying to visualize objects in a 4D phase space —such as invariantmanifolds—, standard projections may give rise to false intersections between them. These artifactsdue to projections must be detected and differentiated from real intersections to ensure or discard theexistence of homoclinic or heteroclinic connections. To do so, we make extensive use of dynamicalsystems theory and topological arguments to justify our findings.With a strategic combination of invariant manifold analysis and bifurcation theory we find thetraveling wave solutions identifying each of them as a specific heteroclinic/homoclinic connection ora periodic orbit in the four-dimensional phase space. We classify these solutions into 12 differentclasses depending on the topological type of the associated orbit. We also determine conditions on themodel parameters so that there is such a particular kind of solution and identify homoclinic chaoticdynamics as one of the sources of complicated behaviour. All in all, this collection of dynamical objectsrepresent desired wave fronts, wave pulses and wave trains, corresponding to bounded solutions of (3)which ensure survival of both populations in the long term.This paper is organized as follows: section 2 presents some preliminary results and background.Local analysis near steady states is included in section 3, while the bifurcation analysis is presentedin section 4. Wave pulses, wave trains and wave fronts are discussed in sections 5, 6 and 7, respec-tively. Section 8 presents a description of the multiple wave fronts found near a focus-focus homoclinicbifurcation. Section 9 discussed the influence of the propagation speed and the diffusion ratio on the ultiple wave solutions in a diffusive predator-prey model 5 occurrence of the different wave pulses. Section 10 analyzes the existence of traveling fronts for theprey population in the absence of predators. Finally, section 11 presents a summary and discussion ofthe main results. In addition, in the Appendix A, we briefly review the method we used to obtain theglobal two-dimensional manifolds.2.
Model reduction, definitions and notation
For future convenience, the first step to study traveling wave solutions in (3) is to make a timerescaling and a change of parameters given, respectively, by t → D t, ( d, s, b, g, a, m ) = (cid:18) D D , D , sβ, sγ, sα, m (cid:19) ∈ R × ]0 , . (4)Thus, we can write system (3), equivalently, as (cid:40) u t = du xx + su ( u − m )(1 − u )( u + v ) − auv ; v t = v xx + buv − gv ( u + v ) . (5)Here, d = D /D represents the ratio of diffusion rates of prey and predators, respectively, and appearsas an explicit parameter in (5). If d > d < z := x + ct , where c > U ( z ) = u ( x, t ), V ( z ) = v ( x, t ). Applying the chain rule and substitutingthis into (5), we obtain the following set of second order ordinary differential equations c dUdz = d d Udz + sU ( U − m )(1 − U )( U + V ) − aU V ; c dVdz = d Vdz + bU V − gV ( U + V ) . (6)Naming the auxiliary variables W = dU/dz and R = dV /dz , system (6) can be expressed equiva-lently as the vector field X : dUdz = W,dVdz = R,dWdz = 1 d ( cW − sU ( U − m )(1 − U )( U + V ) + aU V ) ,dRdz = cR − bU V + gV ( U + V ) . (7)The biologically relevant region of the four-dimensional phase space of (7) is the set Ω = { ( U, V, W, R ) ∈ R : U ≥ , V ≥ } . A traveling wave of (5) is any bounded solution of the system (7) contained inthe realistic domain Ω. These solutions are functions that “travel” in space, preserving their form astime goes by. A sketch of a traveling wave moving to the left with speed c > E. Villar–Sep´ulveda, P. Aguirre and V. Bre˜na–Medina + Figure 1.
Sketch of the spatial propagation of a traveling wave solution, with c >
Pulses and fronts as homoclinic and heteroclinc orbits.
An orbit of the vector field (7) iscalled homoclinic if it connects an equilibrium to itself converging to this point as the orbit parameter z → ±∞ . Figure 2(a) shows an actual homoclinic orbit of (7) found with the method presented laterin section 5. The homoclinic orbit connects the equilibrium point q = ( q u , q v , ,
0) (given explicitlyin section 3) to itself. The time series of U and V associated with this connecting orbit are shownin figure 2(b) in blue and orange, respectively. The profile of this solution is characterized by a largedeviation (or pulse) in the amplitudes of U and V followed by a convergence back to the resting state.This corresponds to a wave pulse in the original system (5) that travels from a spatially homoge-neous stationary solution to itself. Therefore, homoclinic orbits of (7) correspond to wave pulses of(5).Formally, this wave pulse satisfies for all x ∈ R in (5), and respectively in (7), thatlim t →±∞ ( u ( t, x ) , v ( t, x )) = ( q u , q v ) ⇔ lim z →±∞ ( U ( z ) , V ( z ) , W ( z ) , R ( z )) = q . (8)On the other hand, an orbit is said to be heteroclinic if it converges to different equilibra at eachend z → ∞ or z → −∞ . Heteroclinic orbits of (7) correspond to wave fronts of (5). Figure 3(a) showsa heteroclinic orbit connecting equilibrium q to equilibrium p = ( p u , p v , ,
0) (given explicitly below insection 3). This corresponds to a wave of (5) that makes the transition from one spatially homogeneousstationary solution to another as is shown in figure 3(b). Formally, this wave front satisfieslim t →−∞ ( u ( t, x ) , v ( t, x )) = ( q u , q v ) , ∀ x ∈ R , lim t → + ∞ ( u ( t, x ) , v ( t, x )) = ( p u , p v ) , ∀ x ∈ R , in (5) ⇔ lim z →−∞ ( U ( z ) , V ( z ) , W ( z ) , R ( z )) = q , lim z → + ∞ ( U ( z ) , V ( z ) , W ( z ) , R ( z )) = p , in (7). (9) ultiple wave solutions in a diffusive predator-prey model 7 q • V UW (a) (b)
U q u Vq v z Figure 2.
Profile of a wave pulse. Panel (a) shows a homoclinic orbit which joins the equilibrium q to itself in the long term, projected onto the U V W space, while panel (b) shows the time series of U and V associated with this solution in blue and orange, respectively; here, the range of z values in thehorizontal axis is restricted to the interval where the variables U and V develop the pulse. Parametervalues are a = 24, b = 19, g = 1, c = 1, s = 100, m = 0 . d = 1 . z ∈ ] − ∞ , ∞ [. However, for computational procedures, this indepen-dent variable is rescaled to z ∈ ]0 ,
1[ in all our results; see the Appendix A for more details. Moreover,in figure 2 (and for every other wave pulse shown throughout this paper) we restrict the values of z tothose compact subintervals of ]0 ,
1[ where the pulses are easier to see.2.2.
Wave trains and periodic orbits.
The third kind of wave solution of (5) we are interested inis wave trains. These solutions correspond to periodic orbits of system (7), as is illustrated in figure 4.A cycle of (7) is associated with a traveling wave of (5) with a periodic profile in both space and time.More formally, for every periodic orbit in the domain Ω of (7) there exists
T > period )such that (cid:0) u (cid:0) ˆ x, ˆ t (cid:1) , v (cid:0) ˆ x, ˆ t (cid:1)(cid:1) = ( u ( x, t ) , v ( x, t )) in (5) , ⇔ (10) ( U (ˆ z ) , V (ˆ z ) , W (ˆ z ) , R (ˆ z )) = ( U ( z ) , V ( z ) , W ( z ) , R ( z )) in (7) , where ˆ x = x + αT, ˆ t = t + (1 − α ) c T, ˆ z = z + T, and α ∈ [0 , T = 1; see theAppendix A. In particular, in figure 4 (and for every other wave train shown throughout this paper)we restrict the values of z to one such period of the cycle. E. Villar–Sep´ulveda, P. Aguirre and V. Bre˜na–Medina p ❄ •• q U VW (a) (b)
U q u p u Vp v q v z Figure 3.
Profile of a wave front. Panel (a) shows a heteroclinic orbit from q to p , projected ontothe U V W space, while panel (b) shows the time series of U and V associated with the same solution.Parameter values are the same as in figure 2 except for d = 2 . c = 0 . -8-6-4-20 W -4 V U z U V V UW (a) (b)
U Vz
Figure 4.
Profile of a wave train. Panel (a) shows a periodic orbit projected onto the
U V W space,while panel (b) shows the time series of U and V contained in a (rescaled) period of length 1. Parametervalues are the same as in figure 2 except for d = 1 . ultiple wave solutions in a diffusive predator-prey model 9 Invariant manifolds.
Homoclinic and heteroclinic orbits in ordinary differential equations —such as (7)— are intersections of invariant manifolds of equilibria. If all the eigenvalues of the associatedJacobian matrix at an equilibrium, say x ∗ , have non-zero real part, we say that x ∗ is hyperbolic ;otherwise, it is called non-hyperbolic . In the hyperbolic case, if all the eigenvalues have negative (resp.positive) real part, the equilibrium is called an attractor (resp. repeller ). If x ∗ has eigenvalues withpositive and negative real parts, the equilibrium is a saddle . In such case, the stable manifold theorem[33, 43] ensures the existence of the local stable and unstable manifolds of x ∗ defined, respectively, as W s loc ( x ∗ ) = { y ∈ N : ϕ z ( y ) → x ∗ when z → ∞ , and ϕ z ( y ) ∈ N, ∀ z ≥ } ,W u loc ( x ∗ ) = { y ∈ N : ϕ z ( y ) → x ∗ when z → −∞ , and ϕ z ( y ) ∈ N, ∀ z ≤ } , where N is a neighborhood of x ∗ . In our particular case, ϕ z is the flow of the vector field (7), where theorbits are parameterized by the independent variable z . The sets W s loc ( x ∗ ) and W u loc ( x ∗ ) are inmersed,smooth manifolds tangent at x ∗ to the eigenspaces E s ( x ∗ ) and E u ( x ∗ ), respectively. Moreover, thedimension of W s loc ( x ∗ ) (resp. W u loc ( x ∗ )) is the same as that of E s ( x ∗ ) (resp. E u ( x ∗ )). That is,dim W s loc ( x ∗ ) (resp. dim W u loc ( x ∗ )) is the number of eigenvalues with negative (resp. positive) real part.By extending W s loc ( x ∗ ) and W u loc ( x ∗ ) by the flow ϕ z we obtain the global stable and unstable manifolds given, respectively, by W s ( x ∗ ) = { y ∈ R : ϕ z ( y ) → x ∗ as z → ∞} ,W u ( x ∗ ) = { y ∈ R : ϕ z ( y ) → x ∗ as z → −∞} . In the example of figure 2, the homoclinic connection is an orbit in W u ( q ) which comes back to q along the stable manifold W s ( q ), i.e., it is in the intersection W u ( q ) ∩ W s ( q ). Furthermore, in figure 3,the heteroclinic connection is an orbit in W u ( q ) which moves away from q , but it is also contained in W s ( p ) and, hence, heads toward p . That is, this heteroclinic orbit lies in W u ( q ) ∩ W s ( p ).On the other hand, if x ∗ is a non-hyperbolic equilibrium, there is a locally invariant manifold,denoted as W c ( x ∗ ), and known as the center manifold . The manifold W c ( x ∗ ) is tangent at x ∗ to theeigenspace E c ( x ∗ ) and dim W c ( x ∗ ) = dim E c ( x ∗ ), and captures the essential long term dynamics ina neighbourhood of x ∗ ; we refer to [33] for more details.3. Local analysis
System (7) has at most five equilibrium points in the domain Ω, which are given by p = (0 , , , p m = ( m, , , p = (1 , , , p = ( p u , p v , , q = ( q u , q v , , q u = bs (1 + m ) + √ bs ∆2 bs , (11a) q v = ( b − g ) g q u , (11b) p u = bs (1 + m ) − √ bs ∆2 bs , (11c) p v = ( b − g ) g p u , (11d)provided ∆ = bs ( m − − a ( b − g ) ≥ q u , q v , p u , p v ≥
0. Under these conditions, we have thefollowing result on the stability of p , p m and p : Proposition 1.
Let us consider the quantities ∆ m := c + 4( g − b ) m, ∆ m := c − d (1 − m ) m s. Then, system (7) satisfies the following statements:(1) p is an unstable non-hyperbolic equilibrium, with dim( W u ( p )) = 2 and dim( W c ( p )) = 2 .(2) If b < g , then p m and p are hyperbolic saddles, with dim( W s ( p m )) = 1 , dim( W u ( p m )) = 3 , dim( W s ( p )) = 2 and dim( W u ( p )) = 2 . (3) If b > g , then p m is a hyperbolic repeller and p is a hyperbolic saddle with dim( W s ( p )) = 1 and dim( W u ( p )) = 3 . In addition, if ∆ m > and ∆ m > , then p m is a repelling node.Proof. If we denote the vector field (7) as X , its Jacobian matrix evaluated at the points p , p m and p is given, respectively, by: DX ( p ) = cd
00 0 0 c , DX ( p m ) = − (1 − m ) m sd amd cd
00 ( g − b ) m c , and DX ( p ) = − m ) sd ad cd g − b c . Denoting by λ ji the i -th eigenvalue of the equilibrium p j , j ∈ { , , m } , then we see that λ , = 0 ,λ = cd > ,λ = c > . λ m , = c ± (cid:112) ∆ m ,λ m , = c ± (cid:112) ∆ m d , λ , = c ± (cid:112) c + 4( g − b )2 ,λ , = c ± (cid:112) c + 4 d (1 − m ) s d . Since DX ( p ) has two zero and two positive eigenvalues, then p is an unstable non-hyperbolic equi-librium. The remaining statements are direct consequences of the Hartman-Grobman theorem [33]and the stable manifold theorem. As for p m , since 0 < m <
1, then ∆ m = c − d (1 − m ) m s < c .Furthermore, if b < g , then ∆ m = c + 4( g − b ) m > c . Thus we have that λ m = c + (cid:112) ∆ m > , λ m = c − (cid:112) ∆ m < , Re( λ m , ) > , which implies the desired result. On the other hand, if b > g , then ∆ m < c , which implies that λ m >
0. In the particular case that b > g , ∆ m > m >
0, the eigenvalues λ m , of DX ( p m )become real and positive and, hence, p m is a repelling node. Similarly, a sign analysis of the eigenvaluesof DX ( p ) reveals the stability of p and the dimensions of its invariant manifolds. (cid:3) (cid:3) Proposition 2. If b > g and ∆ > , then both equilibria p and q of (7) are in the domain Ω .Proof. It is immediate to see that if ∆ >
0, then p and q exist and are different. To see that p , q ∈ Ω,note that if b > g and ∆ >
0, then bsm > − a ( b − g ) ⇔ bsm > − a ( b − g ) − bsm ⇔ bs ( m + 2 m + 1) > − a ( b − g ) + bs ( m − m + 1) ⇔ b s (1 + m ) > bs ∆ ⇔ bs (1 + m ) > √ bs ∆ ⇔ q u > p u > . Finally, since q v = ( b − g ) q u /g and p v = ( b − g ) p v /g , the result follows. (cid:3) (cid:3) ultiple wave solutions in a diffusive predator-prey model 11 To perform a stability analysis of the equilibria p and q by standard methods is a challenging task.Indeed, the Jacobian matrix of X evaluated at p and q are given, respectively, by DX ( p ) = a a cd a a c and DX ( q ) = b b cd b b c , where a = 2 a (cid:0) bg + g − b (cid:1) + b (cid:16) b ( m − s − ( m + 1) √ bs ∆ (cid:17) bdg p u ,b = 2 a (cid:0) bg + g − b (cid:1) + b (cid:16) b ( m − s + ( m + 1) √ bs ∆ (cid:17) bdg q u , and a = agbd p u , a = − ( b − g ) p v , a = ( b − g ) p u .b = agbd q u , b = − ( b − g ) q v , b = ( b − g ) q u . Hence, the computation of analytic expressions for the eigenvalues of DX ( p ) and DX ( q ) turns outto be a cumbersome goal. However, direct inspection of equilibrium coordinates reveals evidence ofsome local bifurcations. If b = g , then from (11) we have q v = p v = 0, q u = 1, and p u = m , since bs (1 + m ) ± (cid:112) b s ( m − bs = bs (1 + m ) ± bs (1 − m )2 bs = 1 + m ± (1 − m )2 . This implies that q = p and p = p m . Since, according to Lemma 1, there is a stability change for p and p m when b = g , this is an indication of a transcritical bifurcation. The equilibrium points q and p collide and interchange their stability when b = g ; a similar statement follows for p and p m . On theother hand, if ∆ = 0, then q u = p u and q v = p v , so that q = p . Since these two equilibria exist onlyif ∆ >
0, this is evidence of a saddle-node (or fold) bifurcation of equilibrium points. Formal proofsfor these statements require the reduction of (7) into a one-parameter family of center manifolds oneach case, followed by verification of certain genericity conditions; we refer to [33, 43] for more details.However, we opt to omit these proofs in favor of a focus on the analysis of global bifurcations in (7)and the emergence of traveling waves in (5).4.
Bifurcation analysis
In this section we present a bifurcation analysis of (7) performed with the standard continuationpackage
Auto . As a starting point we consider parameters a = 24, b = 19, g = 1 and s = 100 fixedthroughout this section and let d and m to slowly vary. The fixed values of a, b, g and s correspond tothose in [4] after the transformation (4). Moreover, we take the initial value c = 1 for the wave speedin (7) for a first analysis. Later in Section 9, we let c to vary in order to capture its influence on theexistence and properties of the wave solutions.The resulting bifurcation scenario in the ( m, d )-plane is shown in figure 5. Of special importance forus are the curves h p and h q which represent homoclinic bifurcations to p and q , respectively. The curve h p is divided into two segments (labeled as h cp and h sp , respectively) by a codimension two Belyakovhomoclinic point, labeled as B . The right-hand side endpoints of both h cp and h q (marked with × )correspond to the last points where we could obtain convergence of the computed solutions with Auto .We will address the details and consequences of these homoclinic bifurcations in section 5.
IIIIIIVVII IV ❄ VI ✻ VIII ZH •• FN S ❈❈❈❲ R h sp h cp ✻ ♦ ✷ ×× h q ♦ ✷ P DH − p H + p B Figure 5.
Bifurcation diagram of (7) in the ( m.d )-plane space. Parameter values are the same as infigure 2.Let us continue with an inspection of the other elements in figure 5. The first of these remainingingredients is a curve of Hopf bifurcation at the equilibrium p . This bifurcation curve is divided intotwo segments. The first one is a segment of supercritical Hopf bifurcation, labeled as H − p ; the other oneis a segment of subcritical Hopf bifurcation, labeled as H + p . The curve H − p meets a Fold bifurcationcurve F in a quadratic tangency at a codimension two Zero-Hopf bifurcation point, labeled as ZH . Theseparation between the curves H − p and H + p occurs at a codimension two generalized Hopf bifurcationpoint located very close to the ZH point; we choose not to label this generalized Hopf point in figure 5to avoid an overlapping with the ZH point. The curve labeled as P D corresponds to a period doublingbifurcation while
N S is a Neimark-Sacker (or torus) bifurcation curve. These two curves meet at acodimension two strong resonant point R . The horizontal dotted line in figure 5 corresponds to d = 1(or equivalently, D = D in (3)). While this line does not represent any bifurcation, it is useful todistinguish the phenomena encountered above it from that which occurs below it. Indeed, one mustremember that if d > d < ultiple wave solutions in a diffusive predator-prey model 13 The bifurcation curves in figure 5 divide the ( m, d )-plane into the open regions I-VIII. Region VIII isbounded to the left by the curve F , while I is delimited to its right by F , and below by H − p . Region IIis surrounded by the curves H − p , N S , and
P D ; while region III is enclosed by the curves
P D , h + p , and h q . Furthermore, region IV is bounded by the segments h q , h + p , and P D , while region V is surroundedby the curves
P D , h − p , and h + p . Finally, region VI is enclosed by the curves h − p , h + p , P D , N S , and H + p , while region VII is delimited to the right by the curve F and above by H + p .It is relevant to note that the fold curve F corresponds to ∆ = 0 in section 3. Hence, equilibriumpoints p and q exist on the left-hand side of the F curve (regions I-VII). In particular, if ( m, d ) ∈ I,both equilibria are hyperbolic. If ( m, d ) passes through the supercritical Hopf bifurcation curve H − p from region I towards region II, a limit cycle branches out from p . While the Hopf bifurcation issupercritical, this criticality is restricted only to a suitable two-dimensional center manifold where thebifurcation takes place [33, 43]; the resulting periodic orbit in R is, in fact, of saddle type in region II.The stability properties of this cycle remain unchanged until this orbit undergoes a period doublingbifurcation when ( m, d ) ∈ P D . This periodic orbit faces a number of further bifurcations (not shownin figure 5) as the point ( m, d ) moves towards the curve h q where it gives rise to a homoclinic orbit.We will address this transitions again in section 6. On the other hand, if the point ( m, d ) crosses the N S curve from region II into region VI, an invariant torus bifurcates as a periodic orbit undergoes aNeimark-Sacker bifurcation.Our bifurcation diagram in figure 5 is just partially complete. Other codimension two strong reso-nances can be found along the
N S bifurcation curve. Although the complete bifurcation diagram nearthese bifurcation points is yet to be known in its full complexity, one should expect the appearanceof chaotic behavior when the point ( m, d ) is in a neighborhood of the
N S curve; for further details,see [43]. Furthermore, bifurcation theory tells us that there is an infinite number of bifurcation curvesin neighborhoods of both points B and R . However, the full bifurcation picture near each of thesepoints is not fully known from a theoretical point of view [43]. (We will address the complex dynamicsthat emerges due to the point B in section 5 below).5. Homoclinic bifurcations, wave pulses, and chaos
Figure 6 shows three different homoclinic orbits to the equilibrium p in the left-hand side column,and their respective time series in the right-hand side column. The values of parameters ( m, d ) for eachcase correspond to those marked as ♦ , (cid:3) and × along the curve h p in figure 5. In the left-hand sidecolumn of figure 6, the homoclinic trajectories develop a rotational movement near p before makinga large excursion and returning to p ; see the sequence of panels (a1)-(b1)-(c1). The amplitude ofthese oscillations increases as the point ( m, d ) moves to the right along the curve h p . As a result, thecorresponding wave pulses in panels (a2)-(b2)-(c2) feature an initial transient with increasingly largeroscillations —as ( m, d ) moves to the right along h p — around the equilibrium values; in each case, thisculminates in a large pulse before decaying back to the rest state. Indeed, the initial pattern of smalleramplitude oscillations of each wave takes most of a long interval of values of z ∈ ]0 ,
1[ (i.e, it is a “slow”build-up in terms of z ); while the large amplitude pulse occurs in a smaller interval (of order 10 − )of parameter z (i.e, a “fast” discharge). Further, notice that both populations tend to increase anddecrease their densities simultaneously along any given traveling pulse.The existence of the homoclinic orbit to p implies the presence of chaotic dynamics. Let us now statethe main reasons for this claim. For any ( m, d ) in a neighbourhood of the curve h p , the linearization of(7) at the equilibrium p has one (stable) eigenvalue λ s < λ u , ∈ C ,and λ u >
0. In particular, λ u , are complex conjugate with positive real part Re( λ u , ) >
0. Theequilibrium p is called a saddle-focus . Figure 7 shows all the possible values of the eigenvalues of p (in the complex plane) along the computed segment of the homoclinic bifurcation curve h p . Namely,as parameters ( m, d ) are allowed to vary along the computed segment of the curve h p in figure 5, each p • U VW Γ p (a1) ♦ p u p v U Vz (a2) • U VW Γ p (b1) (cid:3) p u p v U Vz (b2) • U VW Γ p (c1) × p u p v U Vz (c2)
Figure 6.
Homoclinic orbits to the equilibrium p along the bifurcation curve h p . Parameter valuesare ( m, d ) = (0 . , . m, d ) = (0 . , . m, d ) = (0 . , . p traces out a curve segment whose plots are shown in figure 7. Among the unstableeigenvalues, the pair λ u , are the closest to the imaginary axis Re( λ ) = 0; hence, we say that λ u , are the leading unstable eigenvalues. In this setting, if we define the so-called saddle quantity as σ = λ s + Re( λ u , ), Shilnikov’s theorems [33, 43] state that if σ >
0, the homoclinic bifurcationis simple or mild . In this simple Shilnikov homoclinic bifurcation, a single (repelling) periodic orbitbifurcates from the homoclinic orbit on one side of the curve h p . On the other hand, if σ <
0, thehomoclinic bifurcation is chaotic and gives rise to a wide range of complicated behaviour in phasespace. More specifically, one can find horseshoe dynamics in return maps defined in a neighbourhoodof the homoclinic orbit. The suspension of the Smale horsehoes forms a hyperbolic invariant chaoticset which contains countably many periodic orbits of saddle-type. The horseshoe dynamics is robustunder small parameter perturbations, i.e., the chaotic dynamics persist if the homoclinic connection isbroken; see [33, 43]. The segments labeled as h sp and h cp in figure 5 correspond to simple and chaoticregimes, respectively, and are separated by the Belyakov point B where σ = 0 [13]. (Actually, at( m, d ) = B , we have Re( λ u , ) = | λ s | ≈ . h sp and h cp alongthe curves for λ u , correspond to σ > σ < B where σ = 0. This same Belyakov transition is shown for the corresponding λ s value as well (and is also labeled as B ).The bifurcation picture near the curve h p in figure 5 is just a partial representation of the fullcomplexity one may encounter in this region of parameter space. Indeed, the saddle periodic orbitsassociated with the invariant chaotic set may also undergo further bifurcations such as period-doublingand torus bifurcations [33, 43]. Moreover, the presence of the chaotic h cp bifurcation and that of theBelyakov point B imply a very complicated structure (not shown) of infinitely many saddle-node andperiod-doubling bifurcations of periodic orbits as well as of subsidiary n -homoclinic orbits. Figure 8 ultiple wave solutions in a diffusive predator-prey model 15 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Re( ) -4-3-2-101234 I m () λ s ( m, d ) B ✻ λ u ( m, d ) λ u ( m, d ) Bσ < σ > λ u ( m, d ) Bσ < σ > Figure 7.
Eigenvalues of the Jacobian matrix DX ( p ) continued for every ( m, d ) along the homoclinicbifurcation curve h p in figure 5.shows a 2-homoclinic orbit to p (in panel (a1)) and a 4-homoclinic orbit to p (in panel (b1)), as wellas their corresponding time series in panels (a2) and (b2), respectively. In general, n -homoclinic orbitsare characterized by making n − n pulses before setting down to the steady state values; see the 2-pulse and 4-pulse waves in panels (a1)and (b1), respectively. Moreover, for each of these subsidiary n -homoclinic orbits, the system exhibitshorseshoe dynamics and chaos as in the original homoclinic scenario.As for the homoclinic bifurcation h q at the equilibrium q , the associated Jacobian matrix of (7)at q has a pair of complex-conjugate stable eigenvalues µ s , ∈ C and a pair of complex-conjugateunstable eigenvalues µ u , ∈ C , with Re( µ s , ) < µ u , ) >
0. The structure of the eigenvaluesof q as a function of ( m, d ) ∈ h q is shown in figure 9. We say that q is a focus-focus or bi-focus . Theresulting homoclinic orbit Γ q features a spiral-type convergence to q as z → ±∞ . Figure 10 showsthree different examples of such homoclinic orbit in the left column, and their respective time seriesin the right one. The values of parameters ( m, d ) for each case correspond to those marked as ♦ , (cid:3) and × along the curve h q in figure 5. In the left-hand side column of figure 10, the amplitude ofthe oscillations increases as the point ( m, d ) moves to the right along the curve h q . As a result, thecorresponding wave pulses in panels (a2)-(b2)-(c2) develop more oscillations —as ( m, d ) moves to theright along h q — before converging to the equilibrium values as z → ∞ . The spirals and oscillationsthat are visible in panels (a1)-(b1)-(c1) and in panels (a2)-(b2)-(c2), respectively, are associated withthe stable eigenvalues µ s , of q . There is another set of oscillations as z → −∞ which are associatedwith the unstable eigenvalues µ u , ; however, since Im( µ u , ) < Im( µ s , ) (see figure 9 again), these spiralsare relatively less pronounced and hard to see in figure 10. Nevertheless, like the case of the homoclinicorbit to p , here both populations tend to increase and decrease their densities simultaneously alongany given traveling pulse. p ❏❏❏❫ • U VW Γ p (a1) p u p v U Vz (a2) • U VW Γ p (b1) p u p v U Vz (b2)
Figure 8.
Panel (a1) shows Γ p , the 2-homoclinic orbit to p , while panel (a2) shows its time se-ries of U and V associated with Γ p . Similarly, panels (b1)-(b2) show a 2-homoclinic orbit andits associated 4-pulse wave, respectively. Parameter values are the same as in figure 2 except for( m, d ) = (0 . , . m, d ) = (0 . , . q induces chaotic dynamics for every( m, d ) ∈ h q . Indeed, the presence of the homoclinic orbit Γ q to a focus-focus equilibrium is accompaniedby horseshoe dynamics in cross sections near Γ q and, hence, an infinite number of saddle periodicorbits in a neighbourhood of Γ q [43, 55]. Furthermore, in this setting, the saddle quantity is definedas σ = Re( µ s , ) + Re( µ u , ). Since σ > m, d ) ∈ h q , it follows that there are no stableperiodic orbits near Γ q [27, 38].In sum, any solution in a neighborhood of either Γ p (in the chaotic case) or Γ q tends to behaveerratically and presents sensitive dependence to initial conditions. The corresponding orbit in thefour-dimensional phase space of (7) spends a long transient visiting a strange hyperbolic invariant setbefore converging to an attractor. Hence, any bounded solution of (7) passing near either Γ p (in the h cp side of the bifurcation) or Γ q is associated with a chaotic traveling wave, making spatial populationdynamics to behave in an unpredictable way [55].Figure 11(a) shows both homoclinic orbits Γ p and Γ q coexisting in phase space, while figures 11(b1)and 11(b2) show the time series of U and V associated with either trajectory. This special configura-tion occurs when the bifurcation curves h cp and h q cross each other at ( m, d ) ≈ (0 . , . q and Γ p , coexist in phase space. Moreover,one obtains the coexistence of both chaotic invariant sets (each associated with one of the homo-clinic trajectories) and, hence, the corresponding erratic behaviour and sensitive dependence on initialconditions of nearby solutions. ultiple wave solutions in a diffusive predator-prey model 17 -1 -0.5 0 0.5 1 1.5 2 Re( ) -1.5-1-0.500.511.5 I m () µ s ( m, d ) µ s ( m, d ) µ u ( m, d ) µ u ( m, d ) Figure 9.
Eigenvalues µ s ( m, d ), µ s ( m, d ), µ u ( m, d ) and µ u ( m, d ) of DX ( q ) continued for every ( m, d )along the homoclinic bifurcation curve h q in figure 5.6. Periodic orbits and wave trains
In this section we study the limit cycles existing in (7), their bifurcations, and their consequencesfor the nature of wave trains.6.1.
Period doubling phenomena.
Let us consider the periodic orbit Γ which originates at thesupercritical Hopf bifurcation H − p and track its successive bifurcations as parameter d is decreasedand m = 0 . m, d ) crosses the P D curve from region II to III, thecycle Γ undergoes a period doubling bifurcation. As ( m, d ) enters region III, Γ changes its stabilityand a secondary limit cycle Γ appears with approximately twice the period of Γ. As parameter d isfurther decreased, additional period doubling events occur (not shown in figure 5). This is illustratedin figure 12. Periodic orbits Γ , Γ and Γ of periods 2, 4, and 8 times that of Γ, respectively,are shown in panels (a1)-(b1)-(c1). Panels (a2)-(b2)-(c2) show one period of the corresponding timeseries of U and V . Here the actual periods of the solutions are rescaled to T = 1 for visualizationand computational purposes (We explain more about the numerical method to obtain these periodicorbits in the Appendix A). As a consequence, as parameter d decreases and system (7) undergoes thissequence of period doubling bifurcations, the associated wave trains in panels (a2)-(b2)-(c2) displayperiodic spatio-temporal patterns with doubling periods.6.2. Transition from wave trains to wave pulses.
It is essential to highlight that when the point( m, d ) crosses the PD curve from region II to region III, the cycle Γ does not disappear but just changesits stability. Figure 13 shows the graph of the period T of this cycle as a function of d . The bifurcationcurve oscillates around the critical value d ∗ = 1 . h q to q occurs. The amplitude of the oscillations decreases rapidly as the homoclinic limit is approached when d tends to d ∗ ; see [43] and the references therein. Indeed, the “snaking” behavior of the bifurcation curve q • UVW Γ q (a1) ♦ q u q v U Vz (a2) • UVW Γ q (b1) (cid:3) q u q v U Vz (b2) • UVW Γ q (c1) × q u q v U Vz (c2)
Figure 10.
Homoclinic orbits to the equilibrium q along the bifurcation curve h q . Parameter valuesare ( m, d ) = (0 . , . m, d ) = (0 . , . m, d ) = (0 . , . A , Γ B , and Γ C , respectively, corresponding to the points A , B and C in figure 13 As d approaches d ∗ , the cycles pass increasingly closer to q (see the sequence ofpanels (a1)-(b1)-(c1) in figure 14). As a result, one obtains wave trains which spend longer transientsclose to the equilibrium values (see the sequence of panels (a2)-(b2)-(c2) in which the period T ofeach cycle is rescaled to 1). Hence, one can think of the homoclinic orbit Γ q (and its correspondingwave pulse) as the limit of this sequence of periodic orbits (resp. wave trains) of increasing period as d → d ∗ . Furthermore, each of the periodic orbits bifurcated from the period doubling phenomena insubsection 6.1 may also increase their periods and undergo a convergence to n -homoclinic orbits ina similar fashion. Some of these secondary homoclinic bifurcations are mentioned before in section 5and shown in figure 8. 7. Heteroclinic connections and wave fronts
Figure 15(a) shows a heteroclinic orbit, labeled as Γ q,p , obtained for ( m, d ) = (0 . , .
4) inregion I. The heteroclinic connection is oriented from q to p as parameter z —that which parameterizesthe curve— is increased. The resulting time series are shown in figure 15(b) and they correspondto a wave front traveling from the steady state q that decays exponentially to p in synchronizedoscillatory fashion. This non-monotonic behavior is explained by the presence of a pair of stable ultiple wave solutions in a diffusive predator-prey model 19 p • U VW Γ p Γ q (a) • q q u q v U Vz (b) p u p v U Vz (c)
Figure 11.
Panel (a) shows a projection of Γ s and Γ p onto the U V W space, when ( m, d ) ≈ (0 . , . ∈ h cp ∩ h q . Meanwhile, panel (b1) (resp. (b2)) shows the time series of U and V rendered in different color tones, associated with Γ s (resp. Γ p ). The other parameter values arethe same as in figure 2.complex-conjugate eigenvalues (with negative real part) of p when parameters ( m, d ) are in region I.The connecting orbit Γ q,p ⊂ W s ( p ) ∩ W u ( q ) lies in the intersection of the global invariant manifolds W s ( p ) and W u ( q ). Namely, Γ q,p is contained in the two-dimensional unstable manifold W u ( q ) —computed with the method in the Appendix A, and represented in figure 15(a) as a transparent redsurface— and approaches p along its three-dimensional stable manifold W s ( p ). Moreover, since W s ( p )and W u ( q ) are, respectively, three and two-dimensional immersed smooth manifolds in R , it followsthat the intersection W s ( p ) ∩ W u ( q ) is transverse [35]. As a consequence, the heteroclinic orbit Γ q,p and, hence, its associated wave front persist under small parameter variations. Under small changesto the parameter values in region I, the resulting wave front may vary the amplitude of its oscillationsand the actual asymptotic values, but the qualitative behavior of the traveling front remains unalteredthroughout.As parameter ( m, d ) crosses the H − p curve from region I to region II, a pair of stable eigenvaluesof p cross the imaginary axis and become unstable; in the process, a limit cycle branches from p in a supercritical Hopf bifurcation. Hence, in region II, W s ( p ) is a one-dimensional manifold andthe connection Γ q,p does not exist. Rather, it is replaced by a heteroclinic orbit that joins q tothe bifurcated cycle. Figure 16(a) shows the bifurcated periodic orbit, labeled as γ after ( m, d ) hasentered region II from region I. The unstable manifold W u ( q ) rolls up around γ and intersects thethree-dimensional stable manifold W s ( γ ) transversally along a heteroclinic orbit, labeled as Γ q,γ . Thisheteroclinic connection is associated with the traveling wave shown in figure 16(b); this is a fronttransitioning from the steady state at q into a periodic pattern around p . The heteroclinic orbit Γ q,γ (and its traveling front) is preserved in an open subset of region II, and it disappears when ( m, d )crosses the PD curve towards region III as γ loses its stability in a period doubling bifurcation. UVW Γ (a1) U Vz (a2)
UVW Γ (b1) U Vz (b2)
UVW Γ (c1) U Vz (c2)
Figure 12.
The different periodic orbits (panels (a1), (b1) and (c1)) and associated wave trains(panels (a2), (b2) and (c2)) emerging from successive period doubling bifurcations. While the periodsin the time series are uniformly rescaled to 1 for computational purposes, the actual periods of thecycles are T = 12 . , T = 24 . , and T = 48 . . Parameter values are d = 1 . d = 1 . d = 1 . m = 0 . m, d ) is in an open subset of regions III, IV, and V, one can find a wave front traveling from p to q . This front travels in the opposite direction to that in figure 15 and, hence, it corresponds to athird kind of wave. The wave front traveling from p to q is shown in figure 17(b). The wave begins atthe steady state p with oscillations of increasing amplitude until it settles at q . This front correspondsto an intersection of the manifolds W s ( q ) and W u ( p ) forming a heteroclinic orbit in the phase spaceof (7); figure 17(a) shows the heteroclinic orbit (labelled as Γ p,q ) and the two-dimensional manifold W s ( q ) of q as a transparent blue surface. The connection Γ p,q is an orbit in the three dimensionalunstable manifold W u ( p ) which lies on W s ( q ) to converge to q .8. Multiple wave fronts at the focus-focus homoclinic bifurcation
In section 5 we described the complicated dynamics that can be found near the focus-focus ho-moclinic bifurcation Γ q when ( m, d ) ∈ h q . One of the consequences of this fact is the appearance ofmultiple wave fronts of type Γ p,q (described in section 7) which coexist with the main wave pulse Γ q .We explain this finding here by direct, close inspection of the invariant manifolds involved.Figure 18 shows the projection of W s ( q ) onto the U V W space, when ( m, d ) ≈ (0 . , . ∈ h q . Also shown is the homoclinic orbit Γ q (red curve). The two-dimensional manifold W s ( q ) is renderedas a transparent blue surface. Some orbits in W s ( q ) lie in the three-dimensional unstable manifold W u ( p ) forming heteroclinic connections. In our computations, we detected seven heteroclinic orbitscontained in W u ( p ) ∩ W s ( q ); the method to detect them is explained in the Appendix A. Most of ultiple wave solutions in a diffusive predator-prey model 21 B • C • A • T dd ∗ Figure 13.
Bifurcation curve of the period T of periodic orbits with respect to d , near the homoclinicbifurcation at p . Here m = 0 . p,q ) in figure 18. Each of these seven heteroclinicorbits correspond to a different wave front traveling from p to q which are present in the system atthe same time. Further, there are also 2- and 3-homoclinic orbits to q in W u ( q ) ∩ W s ( q ) coexistingwith the primary homoclinic orbit and the heteroclinic connections in phase space; we opted to notshow these subsidiary homoclinic orbits in figure 18 for visualization purposes. Figure 19(a) showsthe spatial profiles of the prey U associated with the primary wave pulse Γ q and the secondary 2- and3-pulse waves (labeled as Γ q and Γ q , respectively) associated with the secondary homoclinic orbits. Inturn, figure 19(b) shows three representative traveling fronts associated with the family of heteroclinicorbits in phase space. The existence of these families of wave fronts when ( m, d ) ∈ h q indicates thatthere must be a sequence of associated global bifurcations as ( m, d ) approaches the h q curve. Ateach of these bifurcation events, the manifolds W u ( p ) and W s ( q ) intersect tangentially in R along a(newly created) heteroclinic orbit. As ( m, d ) moves closer to the h q curve, the intersection becomestransversal and the heteroclinic orbit persists under small parameter variations. Similar events happenin the case of the secondary homoclinic orbits with W u ( q ) and W s ( q ) that explain the emergence of2- and 3-pulse waves as parameters ( m, d ) approach the h q curve.9. The influence of propagation speed and diffusion ratio
The study and results reported so far in sections 4–8 were produced with fixed wave speed c = 1.Here we ask ourselves if there is a minimum wave speed needed for the existence of some of the travelingwaves we have found. To this end, we consider the homoclinic orbits Γ p and Γ q existing at the curves h p -509.43 W -4
10 9.42 V U z u U v V -59.430 W -4
10 9.42 V U z u U v V -59.430 W -4
10 9.42 V U z u U v V q • UVW Γ A (a1) q u q v U Vz (a2) q • UVW Γ B (b1) q u q v U Vz (b2) q • UVW Γ C (c1) q u q v U Vz (c2)
Figure 14.
The different periodic orbits Γ A , Γ B and Γ C (in panels (a1), (b1) and (c1)) and associatedwave trains (panels (a2), (b2) and (c2)). While the periods in the time series are uniformly rescaled to1 for computational purposes, the actual periods of the cycles are T = 8 . A , T = 11 . B , and T = 27 . C . Parameter values are d = 1 . d = 1 . d = 1 . m = 0 . h q , respectively, when m ≈ . c and d . Figure 20 showsthe homoclinic bifurcation curves h p and h q in the ( c, d )-plane. The existence of both homoclinic orbitsis determined by a positive correlation between the wave speed c and the diffusion ratio d ; namely, as c decreases the wave pulses exist provided d becomes sufficiently small, or equivalently, as the predatordiffusion coefficient grows larger compared to the prey diffusion coefficient.In the case of h p , the relation between c and d in figure 20 is almost linear for c ≥ .
1. Indeedthe curve h p can be approximated as d ≈ . c + 0 . , for 0 . ≤ c <
1, with a root meansquare error e = 0 . h p is located in thehalfspace d <
1. Hence this kind of pulse wave with small wave speed c < c decreases below0 .
01, the diffusion ratio d drops abruptly in a non linear way in the form d = O ( c / ); see figure 20(c).As c is further decreased, the continuation procedure loses precision and the last point where we getconvergence of the numerical scheme is at c min = 0 . p when ( c, d ) = (0 . , . p performs many low-amplitude turns in W u loc ( p ) before developing the longexcursion. The corresponding wave in panel (b) shows a slow pattern (in terms of z ) of small amplitudeoscillations followed by a fast large amplitude pulse in a small interval (of order 10 − ) of parameter z , similar to typical dynamic behaviours with different time scales. Indeed, if the relation d = O ( c / )still holds for c →
0, then both (6) and (7) become singular as c → d →
0; while these systemsin the singular limit may be studied with tools from geometric singular perturbation theory [26], this ultiple wave solutions in a diffusive predator-prey model 23 p ❄ •• q Γ q,p W u ( q ) U VW (a) (b)
U q u p u Vp v q v z Figure 15.
The connecting orbit Γ q,p ⊂ W s ( p ) ∩ W u ( q ) lies in the intersection of the global invariantmanifolds W s ( p ) and W u ( q ) in panel (a). In panel (b) the associated wave front travels from the steadystate q and its amplitude decays exponentially fast to p showing oscillations. Parameter values are( m, d ) = (0 . , .
4) and the other parameters remain fixed as in figure 2. p ❄ •• q Γ q,γ W u ( q ) ✛ γU VW (a) (b) U q u p u Vp v q v z Figure 16.
The connecting orbit Γ q,γ ⊂ W s ( γ ) ∩ W u ( q ) lies in the intersection of the global invariantmanifolds W s ( γ ) and W u ( q ) in panel (a). The associated wave front travels from the steady state q and adopts a periodic behaviour oscillating around the equilibrium values of p in panel (b). Parametervalues are ( m, d ) = (0 . , .
7) and the other parameters remain fixed as in figure 2. q •• ✻ p Γ p,q W s ( q ) U VW (a) (b)
U q u p u Vp v q v z Figure 17.
The connecting orbit Γ p,q ⊂ W u ( p ) ∩ W s ( q ) lies in the intersection of the global invariantmanifolds W u ( p ) and W s ( q ) in panel (a). In panel (b) the associated wave front travels from the steadystate p and its amplitude increases exponentially fast before settling down at q . Parameter values are( m, d ) = (0 . , .
4) and the other parameters remain fixed as in figure 2.is beyond the scope of this work. Nevertheless, solutions u ( x, t ) = U ( x + ct ), v ( x, t ) = V ( x + ct ) of (5)as c → d = D /D → x .As for the bifurcation curve h q in figure 20(a)-(b), the dependence between c and d is approximatelyquadratic. The homoclinic orbit to the focus-focus q exists whenever c and d satisfy d ≈ . c − . c + 1 . ≤ c ≤
1; the root mean square error of this approximation is e =6 . × − . In particular, the computed segment of the curve h q is located in the half-space d > c < c > c . Indeed, the bifurcation curve h q can be continued down to c min = 0 with d min ≈ . d min > h p ). Inthe limit as c →
0, the resulting wave pulse corresponds to a stationary solution of (5) where bothspecies propagate with diffusion ratio D /D = d min .10. Traveling waves in the absence of the predator
System (7) has two invariant planes given byΠ U = { ( U, V, W, R ) ∈ R : V = R = 0 } , and Π V = { ( U, V, W, R ) ∈ R : U = W = 0 } . The two-dimensional planes Π U and Π V can be understood as the spaces in which prey or predators,respectively, are extinct. In particular, the origin p ∈ Π U ∩ Π V . The restriction of (7) to Π U is given ultiple wave solutions in a diffusive predator-prey model 25 q • p • ✛ Γ q Γ p,q W s ( q )(a) p • ✛ Γ p,q (b) p • ✟✟✟✙ Γ p,q Γ q • q (c) Figure 18.
The stable manifold W s ( q ) of q projected onto the U V W space in panel (a) when( m, d ) ≈ (0 . , . ∈ h q . The manifold W s ( q ) contains a family of coexisting heteroclinicorbits which join the equilibrium p to q as well as the primary focus-focus homoclinic connection Γ q to q . One of such heteroclinic orbits is represented by the trajectory Γ p,q . Panels (b) and (c) showenlargements near the equilibria p and q , respectively. The other parameter values are the same as infigure 2.by X U : dUdz = W,dWdz = 1 d (cid:0) cW − sU ( U − m )(1 − U ) (cid:1) . (12)System (12) has three equilibria: = (0 , m,
0) and (1 , p , p m and p , respectively, to Π U . The equilibrium ( m,
0) is a hyperbolic repeller and (1 ,
0) is ahyperbolic saddle of (12). This result is a direct consequence of Hartman-Grobman theorem. Indeed, p u q u . . . . . Γ q Γ q Γ q . . . . . . . . (a) p u q u . . . . . . . . . . . . . (b) Figure 19.
Spatial profiles of the prey U associated with the wave pulses Γ q , Γ q and Γ q in panel(a); and three representative traveling fronts Γ p,q , Γ p,q , and Γ p,q in panel (b), when ( m, d ) ≈ (0 . , . ∈ h q . The other parameter values are the same as in figure 2.the linear part of (12) is given by DX U ( U, W ) = (cid:32) sd (cid:0) U − m + 1) U + 2 mU (cid:1) cd (cid:33) . Therefore, evaluation of DX U at ( m,
0) and (1 ,
0) is given, respectively, by DX U ( m,
0) = (cid:32) sd m ( m − cd (cid:33) , DX U (1 ,
0) = (cid:32) sd (1 − m ) cd (cid:33) . ultiple wave solutions in a diffusive predator-prey model 27 (a) h p h q cd cd h q (b) h p d c (c) Figure 20.
The homoclinic bifurcation curves h p and h q in the ( c, d )-plane in panel (a). Panels (b)and (c) show enlargements near the curves h q and h p , respectively (the apparently different shapeof each curve is due to the different scales in each plot). Parameter m = 0 . • p Γ p U VW (a) (b)
U Vzp u p v Figure 21.
The homoclinic orbit Γ p to p , when ( c, d ) = (0 . , . Furthermore, the eigenvalues of DX U ( m,
0) and DX U (1 ,
0) are given, respectively, by λ m ± = c ± (cid:112) c − sdm (1 − m )2 d , λ ± = c ± (cid:112) c + 4 sd (1 − m )2 d . Since 0 < m <
1, then we have c − sdm (1 − m ) < c . Therefore, Re( λ m ± ) > m,
0) is a repeller. On the other hand, let us note that c + 4 sd (1 − m ) > c , so λ − < < λ .Thus, (1 ,
0) is a saddle.In the case of the origin (0 ,
0) of (12), we have DX U ( ) = (cid:32) cd (cid:33) , with associated eigenvalues λ = 0 , and λ = c/d > . Hence, is a non-hyperbolic equilibrium. Theeigenvectors of DX U ( ) associated with λ and λ are v = (1 , T and v = ( c, d ) T , respectively.According to the Centre Manifold theorem [33], the origin of (12) has a one-dimensional local centremanifold W c loc ( ) which is tangent to v at . This implies that W c loc ( ) can be represented locally asthe graph of a function W = W ( U ) that satisfies W (0) = W (cid:48) (0) = 0 (for further details, see [43]). Bytaking the Taylor series expansion of this function around U = 0, we have W ( U ) = r (cid:88) k =2 a k U k + O ( U r +1 ) , where the coefficients a k ∈ R , and O ( U r +1 ) are terms of order r + 1 and higher of the Taylor series of W ( U ). Here, the coefficients a k are determined by substitution of W = W ( U ) into (12). Thus, aftersome calculations, we obtain W ( U ) = − msc U + O ( U ) . If we restrict (12) to W c loc ( ), we obtain the one-dimensional differential equation U (cid:48) = W ( U ) = − msc U + O ( U ) . Then, for
U > U (cid:48) <
0. Therefore, is a local attractor in W c loc ( ). Since λ >
0, it follows that is a non-hyperbolic saddle point of (12).Figure 22 shows the phase portrait of (12). The global centre manifold W c ( ) extends itself (for z <
0) along the flow of (12) and forms a heteroclinic connection to ( m, W u ( m,
0) is a two-dimensional unstable manifold,the intersection along W c ( ) ∩ W u ( m,
0) is transversal. It follows that the associated wave front persistsunder small parameter variations. Furthermore, no other saddle connection occurs in (12). Only undersuitable parameter perturbations, the manifolds W u ( ) and W s (1 ,
0) may come to intersect along asecond heteroclinic orbit from (0 ,
0) to (1 , V reveals that there are eitherno homoclinic, heteroclinic orbits nor limit cycles with non-negative V -coordinates in Π V . Therefore,there are no plausible traveling waves of (5) occurring in the absense of prey population.11. Discussion
Upon taking into consideration a diffusion process in model (2), we highlight that its dynamicalrichness considerably raises. Indeed, it made appear the three main types of traveling waves we werelooking for: wave pulses, wave fronts, and wave trains. Most of these orbits allow us to conclude that ultiple wave solutions in a diffusive predator-prey model 29 ( m, ✂✂✂✂✂✍ •• • (1 , W u ( ) W c ( ) W s (1 , W u (1 , UW Figure 22.
The centre manifold W c ( ) extends itself (for z <
0) along the flow of (12) and forms aheteroclinic connection to ( m, c, d, m, s ) = (1 , . , . , (cid:40) u t = du zz − cu z + su ( u − m )(1 − u )( u + v ) − auv,v t = v zz − cv z + buv − gv ( u + v ) , (13)which coincides with system (5), once the change of variables ( x, t ) → ( x + ct, t ) has been made. Byusing this approach, we studied the stability of every shown solution, and we found that each of themis an unstable trajectory. This means that such trajectories will be left in the long term when anyperturbation is applied. For further details, see [58].On the other hand, the presence of chaos in system (7) indicates that setting initial conditionsfor populations that are well modeled by it in the real world may be unreliable as a consequence ofinitial condition sensitivity features. In other words, any change in that condition may give rise toan erratic population behavior as time goes by. As can be seen in sections above, key parametersare crucial as well as pivotal theoretical paradigms that should be met to regulate chaotic traits andhence avoid erratic behavior. Therefore, the multiple wave solutions coexisting with chaos emerge asa manifestation of regular spatiotemporal evolution in the midst of such complex behavior. The factthat we are able to pinpoint the specific associated homo/heteroclinic orbits and periodic solutionsamong chaotic dynamics becomes more relevant in light of these findings. Table 1.
Summary of the main traveling wave solutions.Wave type Region in ( m, d ) plane Orbit in phase spaceType A wave pulse Bifurcation curve h q Focus-focus homoclinic orbit Γ q Type B wave pulse Bifurcation curve h p Saddle-focus homoclinic orbit Γ p Type C wave pulse Neighborhoods of h q Focus-focus 2-homoclinic orbit Γ q Type D wave pulse Neighborhoods of h p Saddle-focus 2-homoclinic orbit Γ p Type E wave pulse Neighborhoods of h p Saddle-focus 4-homoclinic orbit Γ p Type F wave pulse Neighborhoods of h q Focus-focus 3-homoclinic orbit Γ q Type A wave front Region I Heteroclinic orbit from q to p Γ q,p Type B wave front h q and Regions III, IV, V Heteroclinic orbit from p to q Γ p,q Type A wave train Regions II, III, IV, V and VI Limit cycle, ΓType B wave train Regions III and V 2-turns limit cycle, Γ Type C wave train Region III 4-turns limit cycle, Γ Type D wave train Region III 8-turns limit cycle, Γ Regarding the diffusion process, let us notice that the condition d > d < m, d ) belongs to any of the h p , h q and NS curves, or some of their neighborhoods; see figure 5and discussion in section 5.Finally, table 1 shows a list of all the traveling waves exposed throughout this paper, together withthe regions of the bifurcation diagram where they can be found. We highlight that we made this tableaccording to the obtained information, and the mentioned regions may not be the only ones where thegiven traveling waves can be found.It is relevant to remember that there are regions where two or more types of traveling waves existsimultaneously. For example, if ( m, d ) is in the intersection between h q and h p , then there are twohomoclinic orbits, Γ q and Γ p . Furthermore, let us remember that Γ is branched from Γ when ( m, d )crosses the PD curve from region II to III. However, Γ does not disappear after this event. In fact, afterthe successive period doubling bifurcations, a large number of wave trains may be found together in thephase portrait of system (7). Finally, we should recall that when ( m, d ) ≈ (0 . , . ∈ h q ,there is a family of heteroclinic orbits, which go from p to q . These arise from a transverse intersectionbetween W u ( p ) and W s ( q ), which indicates that these solutions are robust under small changes ofparameter values. This result can be interpreted as follows: for suitable initial conditions and Alleethreshold and sensitive diffusion ratio values, traveling fronts mimicking invasion-like traits of bothpopulations are likely able to be observed, and they persist under any sufficiently small perturbation.A relevant fact regarding the homoclinic orbit Γ q is that it presents Shilnikov focus-focus homoclinicchaos in a concrete model vector field. Indeed, most of the studies about this bifurcation have beencarried out (so far) only theoretically [27, 63]. Moreover, as far as we know, this work representsthe first example of the actual computation of a global invariant manifold involved in a focus-focushomoclinic bifurcation in R . ultiple wave solutions in a diffusive predator-prey model 31 On the other hand, according to the summary given in table 1, to obtain an oscillatory behavior ofpopulations in space as time goes by, the condition is to be in any of the regions II-VI of the bifurcationdiagram. In these, many wave trains with different spatiotemporal periods can be found. Finally, toavoid chaos and ensure that both species establish in a stationary state in which they survive in thelong term, then the main requirement for the model parameters is to remain away from the homoclinicbifurcation curves and torus bifurcations.The amplitude of the oscillations varies in markedly different scales for the rescaled variables, andvariations in V are much more pronounced than that of U for pulses, fronts and trains. This differencein amplitudes can be explained by paying attention to the following observations: (i) the prey andthe predator populations spatially spread at distinct rates. This rate difference provides that the low-diffusivity population promotes a heterogeneous aggregation of the two populations along the domain;in consequence, both traveling-wave profiles propagate in a spatially non-homogenous shape, regardlessof d < d >
1; (ii) the prey-predator ratio in the original variables
N/P prevails as both amplitudesare of the same order in any escenarios here considered; (iii) upon integrating system (5) along thewhole finite spatial domain and assuming uniform convergence in time, we obtain that solutions satisfyformulae: L (cid:90) − L [ su ( u − m )(1 − u )( u + v ) − auv ] d x = d B dt , (14a) L (cid:90) − L [ buv − gv ( u + v )] d x = d B dt , (14b)where the total biomasses are given by B ( t ) = (cid:82) L − L u d x and B ( t ) = (cid:82) L − L v d x , regardless of whetherthe boundary conditions for u x and v x vanish or cancel out each other at x = ± L . Moreover, once wedefine the weighed total biomass by B ( t ) := b B ( t ) + a B ( t ), from (14), we obtain formula L (cid:90) − L [ bsu ( u − m )(1 − u ) − agv ] ( u + v ) d x = d B dt . (15)Upon taking into account L (cid:29)
1, we get traveling-wave solutions of (5), which are characterizedby having profiles that keep their shape over time. Now, we consider a moving interval-frame J =[ − L − ct, L − ct ], which “runs alongside” the traveling profiles with the same speed c ≥
0. In so doing,integral in the left-hand side term of (15) over J is constant for all t ≥
0. That is, since traveling-waveprofiles displacement only changes their position in time, the weighed total biomass traveling speed istherefore conserved, i.e. d B /dt ≡ C , where C is constant. Thus, as u + v ≥ | x | ≤ L , thetotal biomasses of the two populations follow a conservation-like property for traveling-wave solutions.Namely, as the wave variable z = x + ct corresponds to a spatial translation for t ≥
0, populationsamplitudes balance each other to satisfy identity (15) for a constant weighed total biomass rate ofchange C .It is interesting to note that identities in (14) are also satisfied for stationary solutions; that is, uponsetting u t = v t = 0, we have relation (15) with d B /dt ≡
0, when homogeneous Neumann boundaryconditions are in place. Similar identities are obtained for Dirichlet, periodic and mixed-type boundaryconditions, as well. Details about stationary solutions of (5) are omitted here, but will appear in anew manuscript based on current research by these authors.
Appendix A. Two-dimensional invariant manifolds as a families of orbit segments
The two-dimensional global manifolds in this paper are computed as families of orbit segments,which can be obtained as solutions of a suitable boundary value problem (BVP), irrespective of thevector field undergoing a homoclinic or heteroclinic bifurcation. This allows us to make use of
Auto toimplement and solve the BVP; then, the manifold is “built up” by continuation of the respective orbitsegments [24, 40, 41]. All images of W s ( p ) and W ss ( q ) are then obtained by rendering the respectiveportion of the manifold as a surface from the computed orbit segments with dedicated Matlab routines.In this appendix we focus on the case of a stable manifold (a similar approach can be adopted tocompute unstable manifolds by considering unstable eigenvalues and reversing time); see also [3, 1,7, 54] for further details and applications. We consider a function u : [0 , → R , that solves thedifferential equation ddτ u ( τ ) = T F ( u ( τ )) , (16a)where F : R → R stands for the vector field (7). Let us notice that (16a) is a scaled version of (7)in which the integration time T of an orbit segment appears now as an explicit parameter. In thisway, the actual integration time over an orbit segment of (16a) is always 1. The function u representsa unique orbit segment { u ( τ ) ∈ R : 0 ≤ τ ≤ } of (16a) provided suitable boundary conditions areimposed at one or both end points u (0) and u (1). (This same approach can be applied to calculateperiodic orbits by setting u (0) = u (1) and imposing additional suitable integral conditions [24]). In thecase of orbits in a stable manifold, these constraints arise naturally from the stable manifold theorem,which states that the stable manifold W s ( q ) of an equilibrium q is tangent to its stable eigenspace E s ( q ) at q [43, 33]. Hence, we can obtain an approximation of W s ( q ) as a one-parameter familyof orbit segments ending at E s ( q ), sufficiently near q . Formally, this family is parameterised usinga suitable sub-manifold of E s ( q ). If the eigenvalues of q are complex conjugate —as is the case of W s ( q ) in (7)—, then we consider the boundary conditions given by u (1) = w s + δ ( w s − w s ) ∈ R , δ ∈ [0 , , (16b)where w s and w s are chosen as follows. Let w s = q + ε v s , where ε > v s ∈ R is any non-zero vector in E s ( q ). Then, we take w s as the first return (backwards in time)of the orbit through w s to the local section generated by v s and by any other linearly independentvector that does not belong to E s ( q ) [7]. Equation (16b) defines an approximate fundamental domain ,as each orbit in W s ( q ) intersects this segment exactly once. Furthermore, this implies that δ ∈ [0 , (cid:99) W sδ,T ( q ) := { u δ ( t ) : t ∈ [0 , } that gives an approximationof W s ( q ).The set (cid:99) W sδ,T ( q ) is a ( δ, T )-dependent family of orbit segments. Indeed, for each T , we have aone-parameter family of orbit segments that approximates W s ( q ). Then, to obtain this approximationin Auto , we need to specify an initial orbit that satisfies (16a) and (16b), with δ = δ ∈ [0 ,
1) fixed.To get this trajectory, we consider the trivial orbit u ≡ w s + δ ( w s − w s ) with T = 0, and performa continuation in T until a desired integration time T = − T ∗ <
0, with u (1) satisfying (16b). Then,a continuation in δ ∈ [0 ,
1) allows us to obtain the desired approximation of W s ( q ) as a family oforbit segments for T = − T ∗ fixed. Moreover, by construction, the family (cid:99) W sδ,T ( q ) lies in an O ( (cid:15) )neighbourhood of W s ( q ).Instead of keeping a particular integration time T fixed, it is possible to restrict the point u (0).This allows us, for instance, to get the intersection W s ( q ) ∩ M with a codimension one submanifold M = { x ∈ R : G ( x ) = 0 } , by imposing another boundary condition given by G ( u (0)) = 0 . (16c) ultiple wave solutions in a diffusive predator-prey model 33 Solutions of (16a) with boundary conditions (16b) and (16c) are orbits that begin at M and end atthe approximate fundamental domain of W s ( q ) near q .To obtain a first trajectory in this case, we continue an orbit that satisfies (16b) for δ = δ ∈ [0 , G ( u (0)). We stop this continuation when condition (16c) issatisfied. This starting orbit segment can then be continued as a one-parameter family of solutions ofthe boundary value problem (16) in the free parameters δ and T , so that the point u (0) traces out aone-dimensional curve on M . The two-dimensional manifolds in figures 15, 16, 17 and 18 were obtainedby setting M to be the sphere S r = { x ∈ R : | x − q | = r } with radius r = 0 .
03, centered at q , where | · | denotes the Euclidean norm in R . We highlight that not every orbit ending at the fundamentaldomain may intersect M . Also, the intersection set W s ( q ) ∩ M could consist of many different curves.In that case, it is necessary to repeat the same procedure for different values of δ ∈ [0 ,
1) to get thewhole set. See [7, 40] for further details.In particular, if W s ( q ) contains a heteroclinic orbit flowing from p to q , such trajectory is approx-imated by a bounded orbit segment in the family (cid:99) W sδ,T ( q ) passing sufficiently close to p . Indeed, if aheteroclinic orbit exists, the two-parameter continuation procedure effectively stops as the integrationtime diverges. In practice, an approximation of such connecting orbit is obtained at some specific δ = δ ∗ ∈ [0 ,
1) with a large integration time T = − T ∗ . A similar criterium can be used to detectsecondary n -homoclinic orbits to q as orbit segments ending near q as the integration time T diverges.For instance, in figures 18 and 19, the fundamental domain δ ∈ [0 ,
1) is divided into 13 sub-segmentsby the values 0 < δ < δ < . . . < δ <
1. The heteroclinic connections correspond to δ ≈ . δ ≈ . δ ≈ . δ ≈ . δ ≈ . δ ≈ . δ ≈ . δ ≈ . δ ≈ . δ ≈ . δ ≈ . δ ≈ . . Conflict of interest
The authors declare that they have no conflict of interest.
References
1. P. Aguirre, E. Doedel, B. Krauskopf, and H. M. Osinga,
Investigating the consequences of global bifurcations fortwo-dimensional invariant manifolds of vector fields , Discr. Cont. Dynam. Syst. (2010), no. 4, 1309–1344.2. P. Aguirre, B. Krauskopf, and H. M Osinga, Global invariant manifolds near homoclinic orbits to a real saddle:(non)orientability and flip bifurcation , SIAM J. Appl. Dynam. Syst. (2013), no. 4, 1803–1846.3. Pablo Aguirre, Bifurcations of two-dimensional global invariant manifolds near a noncentral saddle-node homoclinicorbit , SIAM Journal on Applied Dynamical Systems (2015), no. 3, 1600–1643.4. Pablo Aguirre, Jos´e D Flores, and Eduardo Gonz´alez-Olivares, Bifurcations and global dynamics in a predator–preymodel with a strong allee effect on the prey, and a ratio-dependent functional response , Nonlinear Analysis: RealWorld Applications (2014), 235–249.5. Pablo Aguirre, Eduardo Gonz´alez-Olivares, and Eduardo S´aez, Three limit cycles in a leslie–gower predator-preymodel with additive allee effect , SIAM Journal on Applied Mathematics (2009), no. 5, 1244–1262.6. , Two limit cycles in a leslie–gower predator–prey model with additive allee effect , Nonlinear Analysis: RealWorld Applications (2009), no. 3, 1401–1416.7. Pablo Aguirre, Bernd Krauskopf, and Hinke M Osinga, Global invariant manifolds near a shilnikov homoclinicbifurcation , Journal of Computational Dynamics (2014), no. 1, 1–38.8. H Resit Akcakaya, Roger Arditi, and Lev R Ginzburg, Ratio-dependent predation: an abstraction that works , Ecology (1995), no. 3, 995–1004.9. S. Aly, I. Kim, and D.Sheen, Turing instability for a ratio-dependent predator-prey model with diffusion , AppliedMathematics and Computation (2011), no. 17, 7265–7281.10. Roger Arditi and Lev R Ginzburg,
Coupling in predator-prey dynamics: ratio-dependence. , Journal of theoreticalbiology (1989), no. 3, 311–326.11. Roger Arditi and Henni Saiah,
Empirical evidence of the role of heterogeneity in ratio-dependent consumption ,Ecology (1992), no. 5, 1544–1551.12. Nicolas Baca¨er, A short history of mathematical population dynamics , Springer Science & Business Media, 2011.
13. L. A. Belyakov,
Bifurcation of systems with homoclinic curve of a saddle-focus with saddle quantity zero , Mat. Zam. (1984), 681–689.14. A. Berryman, A. Gutierrez, and R. Arditi, Credible, parsimonious and useful predator-prey models: A reply toAbrams, Gleeson, and Sarnelle , Ecology (1995), 1980–1985.15. V´ıctor Bre˜na-Medina and Alan Champneys, Subcritical turing bifurcation and the morphogenesis of localized pat-terns , Physical Review E (2014), no. 3, 032923.16. Alan R Champneys, Yu A Kuznetsov, and Bj¨orn Sandstede, A numerical toolbox for homoclinic bifurcation analysis ,International Journal of Bifurcation and Chaos (1996), no. 05, 867–887.17. E.A. Coding, M.J. PLank, and SS. Benhamou, Random walk models in biology , J. R. Soc. Interface (2008), no. 25,813–834.18. Dana Contreras Julio and Pablo Aguirre, Allee thresholds and basins of attraction in a predation model with doubleallee effect , Mathematical Methods in the Applied Sciences (2018), no. 7, 2699–2714.19. Dana Contreras-Julio, Pablo Aguirre, Jos´e Mujica, and Olga Vasilieva, Finding strategies to regulate propagation andcontainment of dengue via invariant manifold analysis , SIAM Journal on Applied Dynamical Systems (2020),no. 2, 1392–1437.20. Franck Courchamp, Ludek Berec, and Joanna Gascoigne, Allee effects in ecology and conservation , Oxford UniversityPress, 2008.21. Franck Courchamp, Tim Clutton-Brock, and Bryan Grenfell,
Inverse density dependence and the allee effect , Trendsin ecology & evolution (1999), no. 10, 405–410.22. JHP Dawes and MO Souza, A derivation of holling’s type i, ii and iii functional responses in predator–prey systems ,Journal of theoretical biology (2013), 11–22.23. A. Dhooge, W. Govaerts, Yu. A. Kuznetsov, H. G.E. Meijer, and B. Sautois,
New features of the software matcontfor bifurcation analysis of dynamical systems , Mathematical and Computer Modelling of Dynamical Systems (2008), no. 2, 147–175.24. E. J. Doedel, Lecture notes on numerical analysis of nonlinear equations , Numerical Continuation Methods forDynamical Systems, Springer, 2007, pp. 1–49.25. Eusebius J Doedel, Thomas F Fairgrieve, Bj¨orn Sandstede, Alan R Champneys, Yuri A Kuznetsov, and XianjunWang,
Auto-07p: Continuation and bifurcation software for ordinary differential equations , 2007.26. N. Fenichel,
Geometric singular perturbation theory for ordinary differential equations , J. Differ. Equations (1979), no. 1, 53–98.27. AC Fowler and CT Sparrow, Bifocal homoclinic orbits in four dimensions , Nonlinearity (1991), no. 4, 1159–1182.28. Joanna C Gascoigne and Romuald N Lipcius, Allee effects driven by predation , Journal of Applied Ecology (2004), no. 5, 801–810.29. Eduardo Gonz´alez-Olivares, Betsab´e Gonz´alez-Ya˜nez, Jaime Mena-Lorca, and Jos´e D Flores, Uniqueness of limitcycles and multiple attractors in a gause-type predator-prey model with nonmonotonic functional response and alleeeffect on prey , Mathematical Biosciences and Engineering (2013), 345–367.30. Stephen D Gregory and Franck Courchamp, Safety in numbers: extinction arising from predator-driven allee effects ,Journal of animal ecology (2010), no. 3, 511–514.31. Xiaona Guan, Weiming Wang, and Yongli Cai, Spatiotemporal dynamics of a leslie–gower predator–prey modelincorporating a prey refuge , Nonlinear Analysis: Real World Applications (2011), no. 4, 2385–2395.32. J. Guckenheimer, B. Krauskopf, H. M. Osinga, and B. Sandstede, Invariant manifolds and global bifurcations , Chaos (2015), no. 9, 097604.33. John Guckenheimer and Philip Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of vector fields ,vol. 42, Springer Science & Business Media, 2013.34. AP Gutierrez,
Physiological basis of ratio-dependent predator-prey theory: the metabolic pool model as a paradigm ,Ecology (1992), no. 5, 1552–1563.35. Morris W. Hirsch, Differential topology , Graduate Texts in Mathematics, vol. 33, Springer, 1976.36. Crawford S Holling,
The components of predation as revealed by a study of small-mammal predation of the europeanpine sawfly , The Canadian Entomologist (1959), no. 5, 293–320.37. , Some characteristics of simple types of predation and parasitism , The Canadian Entomologist (1959),no. 7, 385–398.38. Ale Jan Homburg and Bj¨orn Sandstede, Homoclinic and heteroclinic bifurcations in vector fields , Handbook ofDynamical Systems (2010), 379–524.39. B. Krauskopf and T. Rieß, A Lin’s method approach to finding and continuing heteroclinic connections involvingperiodic orbits , Nonlinearity (2008), no. 8, 1655.40. Bernd Krauskopf and Hinke M Osinga, Computing invariant manifolds via the continuation of orbit segments ,Numerical Continuation Methods for Dynamical Systems, Springer, 2007, pp. 117–154. ultiple wave solutions in a diffusive predator-prey model 35
41. Bernd Krauskopf, Hinke M Osinga, Eusebius J Doedel, Michael E Henderson, John Guckenheimer, AlexanderVladimirsky, Michael Dellnitz, and Oliver Junge,
A survey of methods for computing (un) stable manifolds of vectorfields , Modeling And Computations In Dynamical Systems: In Commemoration of the 100th Anniversary of theBirth of John von Neumann, World Scientific, 2006, pp. 67–95.42. Yang Kuang and Edoardo Beretta,
Global qualitative analysis of a ratio-dependent predator–prey system , Journalof mathematical biology (1998), no. 4, 389–406.43. Yuri A Kuznetsov, Elements of applied bifurcation theory , vol. 112, Springer Science & Business Media, 2013.44. M. Lewis and P.Kareiva,
Allee dynamics and the spread of invading organisms , Theo. Popul. Biol. (1993),141–158.45. M.A. Lewis and B. Li, Spreading Speed, Traveling Waves, and Minimal Domain Size in Impulsive Reaction–DiffusionModels , Bull. Math. Biol. (2012), 2383–2402.46. Huiru Li and Haibin Xiao, Traveling wave solutions for diffusive predator–prey type systems with nonlinear densitydependence , Computers & mathematics with applications (2017), no. 10, 2221–2230.47. Xiaobiao Lin, Peixuan Weng, and Chufen Wu, Traveling wave solutions for a predator–prey system with sigmoidalresponse function , Journal of Dynamics and Differential Equations (2011), no. 4, 903–921.48. Robert F Luck, Evaluation of natural enemies for biological control: a behavioral approach , Trends in Ecology &Evolution (1990), no. 6, 196–199.49. Dolnik M., Zhabotinsky A.M., A.B. Rovinsky, and Epstein I.R., Spatio-temporal patterns in a reaction-diffusionsystem with wave instability , Chem. Eng. Sci. (2000), 223–231.50. James D Murray, Mathematical biology ii: Spatial models and biomedical applications , vol. 18, Springer, 2001.51. ,
Mathematical biology: I. an introduction , vol. 17, Springer Science & Business Media, 2007.52. E. Nelson,
Dynamical Theories of Brownian Motion , Princeton University Press, 2976.53. A. Okubo and S.A. Levin,
Diffusion and Ecological Problems: Modern Perspectives , 2nd ed., Springer, New York,NY, 2000.54. Hinke M Osinga,
Two-dimensional invariant manifolds in four-dimensional dynamical systems , Computers &Graphics (2005), no. 2, 289–297.55. IM Ovsyannikov and LP Shil’nikov, Systems with a homoclinic curve of multidimensional saddle-focus type, andspiral chaos , Mathematics of the USSR-Sbornik (1992), no. 2, 415–443.56. Lijuan Qin, Feng Zhang, Wanxiong Wang, and Weixin Song, Interaction between allee effects caused by organism-environment feedback and by other ecological mechanisms , PloS one (2017), no. 3.57. Michael L Rosenzweig, Paradox of enrichment: destabilization of exploitation ecosystems in ecological time , Science (1971), no. 3969, 385–387.58. Bj¨orn Sandstede,
Stability of travelling waves , Handbook of dynamical systems, vol. 2, Elsevier, 2002, pp. 983–1055.59. Jonathan A Sherratt and Matthew J Smith,
Periodic travelling waves in cyclic populations: field studies andreaction–diffusion models , Journal of the Royal Society Interface (2008), no. 22, 483–505.60. L. T. Takahashi, N. A. Maidana, W. C. Ferreira, P. Pulino, and H. M. Yang, Mathematical models for the Aedesaegypti dispersal dynamics: travelling waves by wing and wind , Bull. Math. Biol. (2005), no. 3, 509–528.61. Peter Turchin, Complex population dynamics: a theoretical/empirical synthesis , Princeton university press, 2003.62. Vito Volterra,
Variazioni e fluttuazioni del numero d’individui in specie animali conviventi , C. Ferrari, 1927.63. Stephen Wiggins,
Global bifurcations and chaos: analytical methods , vol. 73, Springer Science & Business Media,2013.64. GSK Wolkowicz,
Bifurcation analysis of a predator-prey system involving group defence , SIAM Journal on AppliedMathematics (1988), no. 3, 592–606.65. Xiujuan Wu, Yong Luo, and Yizheng Hu, Traveling waves in a diffusive predator-prey model incorporating a preyrefuge , Abstract and Applied Analysis, vol. 2014, Hindawi, 2014.66. Dongmei Xiao and Shigui Ruan,
Codimension two bifurcations in a predator–prey system with group defense , Inter-national Journal of Bifurcation and Chaos (2001), no. 08, 2123–2131.67. W. M. S Yamashita, L. T. Takahashi, and G. Chapiro, Traveling wave solutions for the dispersive models describingpopulation dynamics of Aedes aegypti , Math. Comput. Simulat. (2018), 90–99. (cid:63)
Departamento de Matem´atica, Universidad T´ecnica Federico Santa Mar´ıa, Casilla 110-V, Valpara´ıso,Chile.
E-mail address : [email protected], [email protected] † Department of Mathematics, ITAM, R´ıo Hondo 1, Ciudad de M´exico 01080, M´exico.
E-mail address ::