Spatiotemporal Patterns in Mathematical Models for Predator Invasions by Sandra M. Merchant B.Sc., The University of Northern British Columbia, 2001 M.Sc., The University of British Columbia, 2003 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) December 2009 c Sandra M. Merchant 2009 Abstract Much attention has been given to oscillatory reaction-diffusion predator-prey systems recently because, in the wake of predator invasions, they can exhibit complex spatiotemporal patterns, notably wave trains and associated irregular spatiotemporal oscillations, thought to occur in natural systems. This thesis considers the generation and stability of spatiotemporal patterns behind invasion in these models and an extension that includes non-local intraspecific prey competition. In the first part, we study the mechanism by which a single member is selected from a continuous family of wave train solutions behind the invasion. This was first studied by Sherratt (1998), where the author develops a selection criterion that is valid near a supercritical Hopf bifurcation in the kinetics and when the predator and prey diffuse at equal rates. We formulate a “pacemaker” selection criterion that generalizes the criterion of Sherratt (1998), but does not depend on these assumptions. We test this pacemaker criterion on three sample systems and show that it provides a more accurate approximation and can apply to unequal diffusion coefficients. In the second part of the thesis, we study the effect of including nonlocal intraspecific prey competition in these systems. We first study the qualitative effect of non-local competition on the spatiotemporal patterns behind predator invasions in these models, and in a related caricature system. We find that non-local prey competition increases the parameter range for spatiotemporal pattern formation behind invasion, and that this effect is greater for lower kurtosis competition kernels. We also find that sufficiently non-local competition allows the formation of stationary spatially periodic patterns behind invasion. Second, we revisit the selection and stability of wave train solutions. We modify the selection criterion from the first part, ii also applying it to the non-local system, and study how the properties of selected wave trains vary with the standard deviation of the non-local prey competition kernel. We find that the wavelength of selected wave trains decreases with the standard deviation of the non-local kernel and also that unstable wave trains are selected for a larger parameter range, suggesting that spatiotemporal chaos may be more common in highly non-local systems. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . x Dedication xi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Statement of Co-Authorship . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Oscillatory predator-prey systems 1.1.2 Biological invasions 1.1.3 Predator invasions and oscillatory reaction-diffusion 1 . . . . . . . . . . . 2 . . . . . . . . . . . . . . . . . . . 4 systems . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 1 7 1.1.4 Wave train selection and stability . . . . . . . . . . . 11 1.1.5 Non-local interactions . . . . . . . . . . . . . . . . . . 17 Objectives of the thesis . . . . . . . . . . . . . . . . . . . . . 21 1.2.1 Wave train selection in the local system (Chapter 2) . 21 1.2.2 Spatiotemporal patterns behind invasion in the nonlocal system (Chapter 3) . . . . . . . . . . . . . . . . 22 iv 1.2.3 References Wave train selection and stability in the non-local system (Chapter 4) . . . . . . . . . . . . . . . . . . . . . 23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2 Wave Train Selection Behind Invasion Fronts in ReactionDiffusion Predator-Prey Models . . . . . . . . . . . . . . . . 30 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2 Mathematical background . . . . . . . . . . . . . . . . . . . . 33 2.3 Coherent structures and selection in the lambda-omega system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.3.1 Coherent structures . . . . . . . . . . . . . . . . . . . 38 2.3.2 Wave train selection . . . . . . . . . . . . . . . . . . . 40 2.4 Beyond the lambda-omega system: the pacemaker criterion . 43 2.5 The criterion in practice 46 . . . . . . . . . . . . . . . . . . . . 2.5.1 Equal diffusion coefficients . . . . . . . . . . . . . . . 47 2.5.2 Unequal diffusion coefficients . . . . . . . . . . . . . . 54 2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.8 Appendix A - Linear spreading speeds . . . . . . . . . . . . . 61 2.8.1 Linear spreading speeds for the equal diffusion cases . 62 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3 Instabilities and Spatiotemporal Patterns Behind Predator Invasions in Systems with Non-local Prey Competition . . 69 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.2 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.3 Instabilities of the spatially homogeneous steady states . . . 76 . . . . . . . . . . . . . . . . . 79 . . . . . . . . . . . . . . . . 80 3.3.1 Prey-only steady state 3.3.2 Coexistence steady state 3.3.3 A caricature system and the effect of kurtosis . . . . 82 3.4 Non-local competition and predator invasions . . . . . . . . . 90 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 v References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4 Selection and Stability of Wave Trains Behind Predator Invasions in a Model with Non-local Prey Competition . . . 108 4.1 Introduction 4.2 The model and wave train families . . . . . . . . . . . . . . . 114 4.3 Wave train stability . . . . . . . . . . . . . . . . . . . . . . . 118 4.4 . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.3.1 Computation of essential spectra for wave trains . . . 119 4.3.2 Computation of the Eckhaus instability boundary . . 121 Wave train selection . . . . . . . . . . . . . . . . . . . . . . 124 4.4.1 Selection in the non-local system . . . . . . . . . . . . 126 4.4.2 Methods 4.4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . 130 . . . . . . . . . . . . . . . . . . . . . . . . . 128 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 5 Conclusion 5.1 5.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 5.1.1 Wave train selection . . . . . . . . . . . . . . . . . . . 143 5.1.2 The influence of non-local prey competition 5.1.3 Wave train stability . . . . . . . . . . . . . . . . . . . 144 Discussion and future directions References . . . . . 143 . . . . . . . . . . . . . . . . 145 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 vi List of Tables 2.1 Parameter values for models (i)-(iv) . . . . . . . . . . . . . . 47 2.2 Mesh interval and timestep for numerical simulations . . . . . 50 vii List of Figures 1.1 Predator-prey oscillations in the Lotka-Volterra model . . . . 3 1.2 Invasion fronts in the Fisher-KPP equation . . . . . . . . . . 6 1.3 Wavefront formed by predator invasion . . . . . . . . . . . . . 9 1.4 Example of wave train in wake of predator invasion . . . . . . 10 1.5 Hopf bifurcation for system (1.12) and examples of wave trains 13 1.6 Example of irregular oscillations in wake of predator invasion 15 2.1 Wave trains behind invasion fronts . . . . . . . . . . . . . . . 32 2.2 Heteroclinic orbits of (2.11) compared with simulations of the system (2.5) . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.3 Example of a transiently observed selected wave train . . . . 51 2.4 Wave train speeds in numerical simulations compared with predictions for cases with equal diffusion coefficients . . . . . 2.5 52 Wave train speeds in numerical simulations compared with predictions for cases with unequal diffusion coefficients . . . . 55 2.6 Example of an unstable selected wave train . . . . . . . . . . 58 3.1 Spatiotemporal patterns observed behind predator invasions in the local system . . . . . . . . . . . . . . . . . . . . . . . . 75 3.2 Probability density functions for non-local kernels . . . . . . . 77 3.3 Growth rates of perturbations and spatial period of pattern for coexistence steady state . . . . . . . . . . . . . . . . . . . 82 3.4 Wavelengths of stationary prey patterns . . . . . . . . . . . . 86 3.5 Coexistence state stability boundaries for the caricature system 88 3.6 Wavenumbers at the Turing instability boundaries for the coexistence state of the caricature system . . . . . . . . . . . . 89 viii 3.7 Example of stationary spatially periodic pattern in the wake of predator invasion . . . . . . . . . . . . . . . . . . . . . . . 3.8 Example of “very irregular” oscillations in the wake of predator invasion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 91 92 Stability boundaries and observed patterns following predator invasion with the Laplace competition kernel . . . . . . . . . 94 3.10 Stability boundaries and observed patterns following predator invasion with the Gaussian competition kernel . . . . . . . . . 95 3.11 Stability boundaries and observed patterns following predator invasion with the uniform competition kernel . . . . . . . . . 96 3.12 Spatial period of stationary periodic patterns behind invasion 97 3.13 Example of a “patchy” predator invasion . . . . . . . . . . . . 101 4.1 Spatiotemporal patterns behind invasion fronts . . . . . . . . 110 4.2 Stability boundaries for the spatially homogeneous coexistence state . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 4.3 Wavelengths L at Hopf bifurcation points for various σ . . . . 117 4.4 Essential spectra for wave train solutions of (4.4) . . . . . . . 121 4.5 Eckhaus instability boundary for wave trains in B − L space 4.6 Predicted and observed wave trains in B − L space . . . . . . 131 4.7 Selection criterion for σ = 7.0 . . . . . . . . . . . . . . . . . . 133 4.8 Eckhaus instability boundary for selected wave trains in B −σ 123 space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 ix Acknowledgements Most importantly, I would like to express my gratitude to my supervisor, Wayne Nagata, for all of his ideas, guidance and wisdom. My M.Sc. thesis supervisor, Michael Doebeli, has also been a great support, especially by continuing to welcome me in his lab group all of these past years and for suggesting that I look at non-local interactions again. I wish to thank the Natural Sciences and Engineering Research Council of Canada (NSERC) and the International Graduate Training Centre (IGTC) for providing financial assistance that has allowed me to dedicate a large amount of my time to research. Thank you also to the Institute of Applied Mathematics (IAM) for computing resources. I am also indebted to the IAM, Mathbio group and the Doebeli lab (the D-dudes) for providing such a great atmosphere and many opportunities to discuss and present my ideas. I thank Jens Rademacher for introducing me to the idea of numerical continuation that opened so many doors, as well as lending me his codes for essential spectra and answering my many questions about them. I’d also like to thank Mark Lewis for inviting me to the University of Alberta for my IGTC exchange and suggesting the idea for the caricature system of chapter 3. I also greatly appreciate the support and many suggestions from my supervisory committee members: Eric Cytrynbaum, Michael Doebeli, Leah Edelstein-Keshet and Michael Ward. Thank you also to Brian Wetton for advice on numerical methods and to Yue-Xian Li for first piquing my interest in stability and spatiotemporal chaos in Math 560. Finally, I’d like to thank my husband Cameron for around the clock IT support and the innumerable other ways in which he has helped me. x Dedication To my sister, Sheri. xi Statement of Co-Authorship The main content of this thesis is collaborative work with my supervisor Wayne Nagata. The initial idea to study wave trains behind predator invasions was his. The idea to study wave train stability using continuation methods arose from discussions between Wayne Nagata and Jens Rademacher (and eventually myself). The ideas for the wave train selection criteria developed in chapters 2 and 4 were my own. The consideration of non-local interactions in these models was suggested to me by Michael Doebeli. For chapter 2 of the thesis, the design of the study and research were performed jointly by Wayne Nagata and myself. I performed all numerical work and subsequent data analysis as well as wrote the first draft of the manuscript. Wayne Nagata helped revise the manuscript. For chapters 3 and 4, the design of the study and subsequent reseach were primarily my own, although Wayne Nagata provided invaluable guidance throughout and Mark Lewis suggested the moment expansion that leads to the caricature system in chapter 3. I performed all numerical work and data analysis as well as prepared the manuscripts, which, along with chapters 1 and 5, were edited by Wayne Nagata. xii Chapter 1 Introduction Oscillatory reaction-diffusion predator-prey models have recently been a focus of study because, in the wake of predator invasions, they are known to exhibit complex spatiotemporal patterns, in particular wave trains and associated irregular spatiotemporal oscillations, thought to occur in natural populations. This work has focused on two main problems: understanding the mechanism by which spatiotemporal patterns are generated in the wake of predator invasion, and determining the stability of wave train solutions generated behind the invasion. This thesis aims to increase our knowledge in both of these areas. For the oscillatory reaction-diffusion framework, we study how in the invasion wake a single wave train solution is selected from a one-parameter family of possible wave trains. In addition, we consider an extension of the oscillatory reaction-diffusion models that includes non-local intraspecific competition of the prey species. We study both the generation of spatiotemporal patterns behind invasion and the stability of wave train solutions in these non-local models. Specific objectives for each of the main projects, Chapters 2-4 of the thesis, are presented in Section 1.2. First, however, I summarize the relevant background in Section 1.1. 1.1 Background This summary begins with oscillatory populations and corresponding ordinary differential equation models in Section 1.1.1. Section 1.1.2 then provides an introduction to biological invasions and partial differential equation models, followed by models for predator invasion in oscillatory predator-prey systems in Section 1.1.3. Much of this thesis considers two central problems in these systems, the selection and stability of wave trains following inva1 sion, and known results for these are summarized in Section 1.1.4. Finally, Section 1.1.5 introduces the notion of non-local interactions and presents the models considered in Chapters 3 and 4 of the thesis. 1.1.1 Oscillatory predator-prey systems Populations of many species in nature are known to oscillate in time. Classical examples include the Canadian lynx (Lynx canadensis), Fennoscandian voles and lemmings (Microtus spp., Clethrionomys spp. and Lemmus spp.), and red grouse (Lagopus lagopus scoticus) in Scotland, which are known to exhibit temporal cycles of 9-11, 3-5 and 6-11 years, respectively [56, 60]. The cause of temporal oscillations observed in natural populations is a central question in ecology. One hypothesis is that oscillations arise through the interaction of a predator and its prey. Predation of one species by another has been considered a mechanism for oscillatory population dynamics since 1926, when Volterra proposed his famous model to explain the oscillations in predaceous fish catches in the Adriatic [16, 39]. This model consists of the system of coupled ordinary differential equations dh dt dp dt = rh − αhp, = −δp + βhp, (1.1) where h(t) is the density of the prey population at time t, p(t) is the density of the predator, and r, α, δ, and β are positive parameters. In this model prey are assumed to have a constant intrinsic birth rate r and the predator have a constant death rate δ. Prey deaths due to predation are assumed to follow the law of mass action, so that the rate of prey death is proportional, with proportionality constant α, to the product of prey and predator densities. Predators convert captured prey to new predators with a conversion efficiency of β [16, 39]. The Lotka-Volterra predator-prey model does indeed produce oscillations in both the predator and prey populations, as shown in Figure 1.1. However, this system is not structurally stable [39]. That is, arbitrarily small perturbations of the system can have qualitatively different dynam2 1.5 4.5 4 3.5 1 p(t) h(t) 3 2.5 2 0.5 1.5 1 0 0.5 0 20 40 60 80 100 0 20 t 40 60 80 100 t Figure 1.1: Predator-prey oscillations in the Lotka-Volterra model (1.1). Parameter values are r = 1.0, α = 0.5, δ = 0.1 and β = 0.5 and initial condition h(0) = 1.0, p(0) = 1.0. ics, such as cycles with different amplitude or possibly no cycles at all [64]. This behaviour limits the applicability of the Lotka-Volterra model, and ecologists have since searched for structurally stable models that exhibit periodic behaviour. Two modifications to these predator-prey dynamics that are deemed more realistic and allow for structurally stable limit cycles are logistic prey growth and a saturating predator functional response [13, 39]. Three systems that incorporate these assumptions are ∂h ∂t ∂p ∂t = f (h, p), = g(h, p), (1.2) with hp (i) f (h, p) = h(1 − h) − h+C , hp g(h, p) = −Bp + A h+C hp (ii) f (h, p) = h(1 − h) − A h+C , p g(h, p) = Bp 1 − h (iii) f (h, p) = h(1 − h) − p 1 − e−Ch , g(h, p) = Bp A − 1 − Ae−Ch (1.3) where A, B and C are parameters. These models all have logistic prey growth in the absence of the predator, but differ in how the saturating func3 tional response and resulting predator growth are modelled [39]. Model (i) is known as the Rosenzweig-MacArthur model, introduced in [46]. Model (ii) is known as the May model, proposed in [38] and (iii) is known as the Ivlev model, given in [28]. Unlike (1.1), these systems have been nondimensionalized by rescaling time and prey and predator densities to minimize the number of parameters in the “kinetics” (i)-(iii), and so A, B and C in fact represent products and ratios of various biological parameters. The system with any of the kinetics (1.3)(i)-(iii) possesses a prey-only steady state (h, p) ≡ (1, 0) where the prey exists at carrying capacity and the predator is absent, and a coexistence steady state (h, p) ≡ (h∗ , p∗ ) where both species exist at some non-zero levels. For each model, as a kinetic parameter is varied there is a supercritical Hopf bifurcation of the coexistence steady state to stable limit cycle solutions, and so these models have proved useful for studying oscillatory predator-prey systems [54–56]. 1.1.2 Biological invasions Biological invasion is the term used to describe the introduction and subsequent spread of an alien species into a novel environment. While some biological invasions are naturally occurring, developments in our civilization have made it so that, today, most introductions, both accidental and intentional, are made by humans [65]. Correspondingly, the frequency of biological invasions has increased dramatically over the past several decades. In many cases, invasive species destroy native communities and damage populations of benefit to humans. One famous invader is the Gypsy moth (Lymantria dispar ). Introduced accidentally around 1870, the Gypsy moth has since spread across much of the North American continent, damaging agricultural crops and resulting in massive economic losses [42]. Other invasive species that have caused much trouble for humans include rabies in central Europe, the brown tree snake (Boiga irregularis) in Guam and Africanized bees (Apis mellifera) in the Americas [24, 65]. Due to their negative impacts on native species and associated human interests, invasive species have been the subject of study for most of the 4 past century. One of the first attempts at modelling the spread of invasive species is attributed to Skellam (1951), who modelled the spread of the muskrat (Ondatra zibethica) across Europe following its accidental introduction in 1905 [42, 57]. These first models were reaction-diffusion systems for a single species, which couple traditional population dynamics with diffusive dispersal, and take the form ∂n = D∆n + f (n), ∂t (1.4) where n(x, t) is the population density at spatial location x in time t, f (n) describes the local population dynamics and the diffusion coefficient D measures the rate of dispersal. Diffusive dispersal assumes that movement is via Brownian random motion and if individuals are released at a central point the resulting distribution in space after time t is Gaussian with mean squared displacement 2tD, 4tD and 6tD for 1, 2 and 3 spatial dimensions, respectively [3, 13, 26]. Since the initial work by Skellam (1951), an array of other mathematical frameworks have been applied to model the spread of invasive species, including integro-difference, discrete space, and stochastic models [6, 30, 33, 37, 42]. Each modelling framework involves different assumptions on the nature of the system being studied and provides valuable insight into the spread of invasive species. The work in this thesis is based on the original reaction-diffusion framework. On infinite domains, these models typically possess travelling wave solutions, which are translationally invariant solutions that move in a single direction at a constant speed. For modelling invasion, travelling front solutions are particularly important, since local introduction of a species on finite domains in these models generally evolves to an invasion front with a spreading speed that asymptotically approaches that of a travelling front. An example of an invasion front for the model (1.4) in one spatial dimension with logistic dynamics f (n) = n(1−n), known as the Fisher-KPP equation, is shown in Figure 1.2 [14]. The asymptotic speed at which fronts invading unstable states travel has been the subject of much study [1, 22]. A review of this topic is provided 5 time → n(x,t) 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0 100 200 300 400 x 500 600 700 800 Figure 1.2: Invasion fronts in the Fisher-KPP equation. The solution for n(x, t) is plotted at ten equally spaced times as a function of space x, with vertical separation proportional to the time interval between solutions. From numerical simulation of the system (1.4) in one spatial dimension with f (n) = n(1−n) and D = 1. Initial conditions were n(x, 0) = 1 on [0, 50] and zero elsewhere. Dirichlet boundary conditions were used, with n(0, t) = 1, n(800, t) = 0. by van Saarloos (2003), where the author defines a linear spreading speed ν ∗ [61]. This speed is obtained by linearization of the system about the spatially homogeneous steady state ahead of the front. More specifically, ν ∗ ∈ R is given by the largest magnitude “dynamically relevant” solution of the saddle point equations ∗ ∗ 0 = S(k , ω ) 0 = (∂k + ν ∗ ∂ω )S|(k∗ ,ω∗ ) 0 = =(ω ∗ − ν ∗ k ∗ ), (1.5) solved for (k ∗ , ω ∗ , ν ∗ ) where S(k, ω) = 0 is the characteristic equation obtained by assuming the ansatz n ∼ ei(kx−ωt) , k, ω ∈ C in the linearization about the unstable steady state ahead of the front (see [61] for details). For example, for fronts invading an unstable state n = 0 in systems of the form (1.4) with one spatial dimension, the linearized system is given by nt = Dnxx + f 0 (0)n. Substitution of the ansatz n(x, t) = n̄ei(kx−ωt) then gives the characteristic equation where S(k, ω) = −Dk 2 + f 0 (0) + iω. Solution of system (1.5) then 6 gives k ∗ = i p f 0 (0)/D, ω ∗ = 2if 0 (0) and the well-known linear spreading speed ν∗ = p 4f 0 (0)D, for rightward travelling fronts. Fronts propagating into unstable states are grouped into two classes: “pulled” and “pushed” fronts. Pulled fronts are those that travel at speed ν ∗ and for (1.4) it has been shown that various conditions on the form of f (n) (for instance f (0) = f (1) = 0, f (n) > 0 in (0, 1) and f 0 (n) assumes its maximum for n = 0) are sufficient to ensure that compact initial conditions converge asymptotically to fronts with spreading speed ν ∗ [22]. Pushed fronts are those that travel at a speed ν > ν ∗ . Pushed waves occur when the front is driven forward not by population growth at the leading edge, but rather by growth at intermediate densities, for instance when the focal species exhibits an Allee effect in its growth function f (n) [22, 36, 41]. 1.1.3 Predator invasions and oscillatory reaction-diffusion systems In the reaction-diffusion literature, a few studies have also considered biological invasions in systems where in addition to the focal species other interacting species are also included in the model [23]. Multi-species models for invasion when the interactions are competitive or cooperative have been studied by various authors (for instance [40, 63]), but for this thesis we focus on invasions in predator-prey systems. The general form of a reaction-diffusion predator-prey system in one spatial dimension is ∂h ∂t ∂p ∂t 2 = Dh ∂∂xh2 + f (h, p), 2 ∂ p = Dp ∂x 2 + g(h, p), (1.6) where h(x, t) and p(x, t) are the densities of the predator and prey, respectively at position x in time t and Dh and Dp are the respective diffusion coefficients. The functions f (h, p) and g(h, p), which we refer to as the “kinetics,” describe the growth and interaction of predator and prey as in the 7 non-spatial model (1.2) and depend on space and time only through the densities h(x, t) and p(x, t). To be relevant for modelling invasion, it is generally assumed that the system (1.6) supports a spatially homogeneous prey-only state (h, p) ≡ (1, 0), where the predator is absent and the prey is at carrying capacity everywhere on the domain, and a spatially homogeneous coexistence state (h, p) ≡ (h∗ , p∗ ), where the predator and prey coexist at some intermediate levels. For predator-prey systems, there have recently been two main invasion scenarios considered. The first is the study of predators as biological control, and the conditions under which the introduction of a (faster) predator can slow, stop or even reverse the spread of an invasive prey species [41]. In this case the prey species has partially invaded the spatial domain when the predator is introduced and eventually catches up to the prey invasion front. The second invasion scenario is the one studied in this thesis. In this scenario the predator invades a spatial domain that is already completely occupied by the prey species. In this case, initial conditions are often taken as the prey-only steady state everywhere on the spatial domain, except for a localized region where the predator is introduced. For instance, h(x, 0) = 1, p(x, 0) = 0 if −l < x < l (1.7) otherwise where here l is the spatial extent of an initial predator introduction at density . As for single-species systems, such predator invasions typically take the form of invasion fronts, travelling away from the site of introduction, an example of which is shown in Figure 1.3. Note that the prey species also forms a front as it is invaded by the predator. As for the single-species case the asymptotic speed in the case of a pulled front is determined by the linear spreading speed ν ∗ from (1.5). Of course, this linear spreading speed will depend on the predator growth rate as well as its diffusion [26]. Recently, a great deal of attention has been given to predator invasions in reaction-diffusion systems (1.6) with oscillatory kinetics such as those in (1.3)(i) − (iii), because these have been shown to exhibit spatiotemporal 8 1.6 time → h(x,t) 1.2 0.8 0.4 0 0 100 200 300 400 500 300 400 500 x 0.8 time → p(x,t) 0.6 0.4 0.2 0 0 100 200 x Figure 1.3: Wavefront formed by predator invasion. The solutions for h(x, t) and p(x, t) are plotted at times t = 550, 555, · · · , 600 as a function of space x, with vertical separation proportional to the time interval between solutions. Shown is the right side of the spatial domain for a numerical simulation of (1.6) with kinetics (1.3)(i), initial conditions (1.7) with l = 50 and = p∗ , no-flux boundary conditions and parameters A = 0.15, B = 0.11, C = 0.2, and Dh = Dp = 1. patterns called “wave trains” in the wake of the invasion. The study of these patterns behind invasion is in part motivated by recent empirical work suggesting similar spatiotemporal patterns exist in natural populations. Field data indicative of temporal cycles in natural populations have been collected and studied for many decades. The most extensively studied such data set is the annual fur return data collected by the Hudson’s Bay Company for the Canadian lynx (Lynx canadensis), which contains over a century of population data [60]. In contrast, the inherent difficulties of collecting and analyzing spatiotemporal datasets have delayed the study of spatiotemporal patterns until the past decade or so. In this recent period, however, a number of field studies have found that temporal cycles previously observed in some populations are not synchronized over space, meaning that the cycles are at different phases at different spatial locations. Populations that have been demonstrated to have such asynchronous cycles include Kielder for9 est field voles (Microtus agrestis), Fennoscandian voles (Microtus spp. and Clethrionomys spp.) and the larch budmoth (Zeiraphera diniana) [56]. An obvious question now is the cause of spatial asynchrony in such oscillatory populations. If temporally oscillatory populations also vary periodically in one spatial direction then they may take the appearance of a spatially periodic density profile that travels with constant speed through space. An example of such a solution is shown in the region behind the predator invasion front in Figure 1.4. For reaction-diffusion systems on the real line, spatially periodic travelling wave solutions are also known as wave trains, and have a spatial period L and at a fixed location x the population density oscillates temporally with a period T . The wave train speed is given by the ratio of these periods, cwtrain = L/T . Since they are finite in extent, and hence not strictly speaking travelling waves, for convenience we refer also to solutions of the form observed behind the invasion in Figure 1.4 as wave trains. Since such wave train solutions have been shown to evolve behind predator invasions in a range of oscillatory reaction-diffusion systems, they have been proposed as a mechanism for the generation of these patterns in nature. 0.25 time → p(x,t) 0.2 0.15 0.1 0.05 0 0 200 400 600 800 x 1000 1200 1400 1600 Figure 1.4: Example of wave train in wake of predator invasion. The solution for p(x, t) is plotted at times t = 1900, 1910, · · · , 2000 as a function of space x, with vertical separation proportional to the time interval between solutions. Only the right side of the spatial domain is shown. From numerical simulation of system (1.6) with initial conditions (1.7) with l = 50 and = p∗ , no-flux boundaries, kinetics (1.3)(ii) and parameters A = 6.92, B = 0.005, C = 0.64, Dh = Dp = 1. 10 1.1.4 Wave train selection and stability Oscillatory reaction-diffusion predator-prey models that couple oscillatory kinetics such as (1.3)(i)-(iii) with diffusive dispersal are known to possess wave train solutions. The existence of wave train solutions for reactiondiffusion systems was first shown by Kopell and Howard (1973) [32]. In this study, the authors established a number of analytical results for a special class of systems called λ−ω systems. In one spatial dimension these systems have the form ∂u ∂t ∂v ∂t = = ∂2u ∂x2 ∂2v ∂x2 + λ(r)u − ω(r)v, + ω(r)u − λ(r)v, (1.8) where r2 = u2 + v 2 and λ(r), ω(r) are real functions of r [32]. For general oscillatory reaction-diffusion systems, we note that near the Hopf bifurcation in the kinetics, the kinetics can be “reduced to normal form” by applying appropriate transformations and neglecting higher order terms [34, 53]. The essential idea behind such a reduction is that the behaviour of the system near the bifurcation point can be approximated by this corresponding simpler system. In the case of equal diffusion coefficients, this corresponding system is a λ − ω system, as studied by Kopell and Howard (1973), and has the general form ∂u ∂t ∂v ∂t = = ∂2u ∂x2 ∂2v ∂x2 + λ0 u − ω0 v − (λ1 u + ω1 v)(u2 + v 2 ), + ω0 u + λ0 v + (ω1 u − λ1 v)(u2 + v 2 ), (1.9) where λ0 is a linear function of the kinetic bifurcation parameter and ω0 , ω1 and λ1 are parameters related to the original kinetic parameters but independent of the bifurcation parameter. Kopell and Howard (1973) have shown that the system (1.9) has a oneparameter family, parameterized by amplitude r, of wave train solutions that can be explicitly written down [32]. These solutions are of the form u = r cos((ω0 + ω1 r2 )t ± v = r sin((ω0 + ω1 r2 )t ± √ √ λ0 − λ1 r2 x), λ0 − λ1 r2 x). (1.10) 11 The plane wave with amplitude r therefore moves at speed ω0 + ω1 r 2 c̄ = ± √ . λ0 − λ1 r2 (1.11) It is also known that families of wave train solutions exist for generic reaction-diffusion systems, for instance with unequal diffusion coefficients [12, 27, 45]. Wave trains with phase speed c̄ are stationary in a co-moving frame. For solutions stationary in a frame moving with speed c̄, system (1.6) may be reduced to a system of first order ordinary differential equations by making the coordinate transformation ξ = x − c̄t and introducing new variables y = h0 (ξ) and z = p0 (ξ) to get dh dξ dp dξ dy dξ dz dξ = y, = z, = = 1 Dh 1 Dp (−c̄y − f (h, p)) , (1.12) (−c̄z − g(h, p)) , where this system has equilibria (h, p, y, z) = (1, 0, 0, 0) and (h, p, y, z) = (h∗ , p∗ , 0, 0) that are clearly related to the equilibria of (1.6). Wave train solutions to (1.6) of speed c̄ are thus periodic solutions of system (1.12). In general these periodic solutions are non-sinusoidal and cannot be expressed analytically. For oscillatory kinetics, however, at fixed speed c̄ the system (1.12) undergoes a Hopf bifurcation as the kinetic parameters are varied and periodic solutions of system (1.12) can be numerically computed by continuation from this Hopf bifurcation point. These Hopf points in the B − c̄ plane and some representative periodic solutions (wave trains of speed c̄) are shown in Figure 1.5 for kinetics (1.3)(ii). Bifurcation at infinite speed corresponds to the Hopf bifurcation in the kinetics-only system. Wave train solutions for parameters noted in the left panel were computed with the numerical continuation software auto [11]. Note that for given kinetic parameter values beyond the Hopf bifurcation there exists a continuous (oneparameter) family of wave train solutions, parameterized by amplitude, as for the λ − ω system (1.9). Each member of this family has a temporal 12 period T , spatial period L and corresponding speed c̄. 0 (i) L=209.808, T=-69.94 p(ξ) 0.14 * (i) 0.06 -5 0 500 1000 1500 2000 ξ (ii) L=704.111, T=-70.41 0.14 * (ii) -10 p(ξ) wave train speed 0.1 0.1 0.06 0 500 1000 1500 2000 ξ -15 (iii) L=1268.01, T=-70.45 p(ξ) 0.14 * (iii) -20 0.1 0.06 0.005 0.015 B 0.025 0 500 1000 1500 2000 ξ Figure 1.5: Hopf bifurcation for system (1.12) and examples of wave trains. The curve in the left panel is the location in B − c̄ space of Hopf bifurcation points for system (1.12) with kinetics (1.3)(ii) and parameters A = 6.92, C = 0.64 and Dh = Dp = 1. Wave train solutions exist to the left of this curve. On the right are the predator density p(ξ, t) for numerically computed wave trains at parameter values shown on the left with B = 0.01 and speeds (i) c̄ = −3, (ii) c̄ = −10, (iii) c̄ = −18. Numerical simulations of the system (1.6) with any of the oscillatory kinetics (1.3)(i)−(iii) for suitable parameter values and compact initial data have been shown by various authors to exhibit wave train solutions in the wake of invasion [54, 55]. However, as we can observe in Figure 1.4 in a given simulation a single member is somehow selected from the one-parameter family of possible wave train solutions. Also, this selection appears to be robust to perturbations of initial and boundary conditions. The members of the wave train family differ in ecologically important characteristics, for 13 instance amplitude of cycles, spatial and temporal periods and the presence or absence of extended low phases. It is therefore an important problem to be able to predict which member of the wave train family will evolve for given kinetic parameters. This wave train selection problem was studied by Sherratt (1998) using the λ − ω system (1.9) as an approximation to general oscillatory reactiondiffusion systems near the Hopf bifurcation in the kinetics [52]. For the system (1.9), Sherratt (1998) predicts that the wave train selected is the solution (1.10) with amplitude [52] 2λ0 r= ω12 1 q 2 2 2 λ1 + ω1 − λ1 . (1.13) Wave trains selected behind predator invasions are not always stable. For instance, Figure 1.6 shows an unstable wave train behind the invasion front. In this case, the selected wave train is only observed transiently in a spatial region immediately behind the invasion front and breaks up into irregular oscillations further back. The instability of selected wave trains has been proposed as a route to spatiotemporal chaos in oscillatory reaction-diffusion systems [54, 55]. Irregular oscillations are observed in many natural populations, and have been studied extensively for oscillatory reaction-diffusion predator-prey systems. In these systems, various initial conditions have been shown to produce irregular oscillations, and several numerical studies give evidence that they may indeed be genuinely chaotic [43, 51, 55]. The stability of selected wave train solutions is clearly important for predictions of spatiotemporal patterns expected in natural populations. In addition, study of long term behaviour after the predator has invaded the entire spatial domain (with no-flux boundaries) finds a significant difference in the fate of regular versus irregular oscillations: regular oscillations (wave trains) damp to spatially homogeneous temporal oscillations, whereas irregular oscillations persist for very long times [29, 43]. The stability of selected wave trains may therefore have important implications for the long term dynamics of natural systems. 14 time → p(x,t) 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 500 1000 x 1500 2000 Figure 1.6: Example of irregular oscillations in wake of predator invasion. The solution for p(x, t) is plotted at times t = 4980, 4982, · · · , 5000 as a function of space x, with vertical separation proportional to the time interval between solutions. Only the right side of the spatial domain is shown. From numerical simulation of system (1.6) with kinetics (1.3)(i), initial conditions (1.7) with l = 50 and = p∗ , no-flux boundaries and parameters A = 0.15, B = 0.05, C = 0.2, Dh = Dp = 1. The linear stability of wave train solutions of reaction-diffusion systems on the real line can be studied by linearizing the system in a comoving frame about the wave train. For instance, consider the reaction-diffusion system in a comoving frame of speed c̄ ∂u ∂2u ∂u = D 2 + c̄ + f (u), ∂t ∂ξ ∂ξ ξ = x − c̄t ∈ R (1.14) where N is the number of species, u is an N -component function of time t and the moving spatial coordinate ξ, D is a positive diagonal matrix of diffusion coefficients and f is smooth. Let u∗ (ξ) be a spatially periodic stationary solution of (1.14), and hence a wave train solution of speed c̄ in the stationary frame. The linear operator L is obtained by linearizing system (1.14) about the stationary solution u∗ (ξ), and so takes the form L = D∂ξξ + c̄∂ξ + ∂u f (u∗ (ξ)). (1.15) The linear stability of the wave train u∗ (ξ) as an equilibrium of (1.14) is then determined by the spectrum of the operator (1.15). We note that since wave trains are translation invariant, the origin in the complex plane is always an element of the spectrum of L. Wave trains are linearly stable if all other 15 elements of the spectrum are contained in the left half-plane, and unstable otherwise [47]. One of the advantages of the λ − ω system (1.9) is that for the wave train solutions (1.10) linear stability can be easily determined. From Kopell and Howard (1973), the wave train solution (1.10) of amplitude r is linearly stable if and only if r2 ≥ 2λ0 (λ21 + ω12 ) . λ1 (3λ21 + 2ω12 ) (1.16) Note that this means that, for the λ − ω system, plane waves with small amplitude are always unstable, whereas large amplitude plane waves are linearly stable. For wave train solutions of general reaction-diffusion systems there is no such simple analytical criterion. Recently, however, methods have been developed for computing the spectrum of operators of the form (1.15) using numerical continuation [5, 44, 48]. A full account of the method, including the implementation of the algorithms with the software auto can be found in Rademacher et al. (2007) [44]. Using these new methods, Smith and Sherratt (2007) have computed the linear stability of wave train solutions for the system (1.6) with kinetics (1.3)(i) as well as for the λ − ω system [59]. This allows for the division of the wave train family into stable and unstable parts and helps us understand the patterns observed behind predator invasions, in particular when irregular oscillations may be expected. Instabilities of wave trains on the real line can be grouped into two classes: “convective” instabilities and “absolute” instabilities. An absolute instability corresponds to the growth of perturbations with time at every fixed point in the domain. In the case of convective instabilities, however, even though the overall norm of the perturbation grows in time, perturbations decay at every fixed point in the domain, because as they grow they are transported, or convected, towards infinity [49]. On finite domains with separated boundary conditions, therefore, convectively unstable solutions may in fact appear stable, because perturbations may convect away through the boundary, whereas absolute instabilities are expected to always manifest [49]. The essential spectrum, as described above does not distin16 guish between convective and absolute instabilities, and should therefore be viewed as a conservative predictor of wave train stability on finite domains, in that wave trains with stable essential spectra should remain stable upon domain truncation, but those with unstable essential spectra may appear either stable or unstable. In fact, in a recent study Smith et al. (2009) have found that, for λ − ω systems absolute stability is the appropriate criterion for predicting whether wave trains will exhibit irregular spatiotemporal dynamics [58]. In their study, Smith et al. (2009) distinguish between convective and absolute instabilities by computation of the so-called “absolute spectra” for wave train solutions [58]. While techniques for computation of absolute spectra do exist, it is substantially more complex for general reaction-diffusion systems, and so in this thesis I consider only the more conservative essential spectra for wave trains [44]. 1.1.5 Non-local interactions As shown above, reaction-diffusion systems are widely used to model spatial predator-prey dynamics, in part because they are amenable to analysis with standard mathematical methods. A standard assumption of these models, however, is that the only spatial process is the random Brownian motion of individuals. In particular, the kinetics functions f (h, p) and g(h, p) are “local”, meaning that interactions and resulting growth rates at a position x are determined only by the predator and prey densities at this same position. In some systems, for instance when individuals compete for a rapidly equilibrating common resource, this may not be a reasonable assumption [15, 17]. In this case it may be more reasonable for growth or death rates arising from the interaction to depend also on neighbouring population densities, or even on the full spatial distribution of the population. These “non-local” interactions have been considered in several contexts, including reaction-diffusion, reaction-convection and metapopulation models [7, 10, 35]. Perhaps the most basic example of non-local interactions in the reaction-diffusion framework is the “non-local Fisher equation,” which models mortality due to intraspecific competition in a non-local manner. The non-local Fisher equa- 17 tion is given by ∂h ∂2h + h(1 − hef f ), = ∂t ∂x2 x ∈ (−∞, ∞) (1.17) where h(x, t) is the population density at x in time t and hef f (x, t) is the effective population density at (x, t) and can depend on population densities at other spatial locations. A common practice is to use a weighted average of population densities over the spatial domain, and so take hef f to be a spatial convolution of the population density distribution, i.e. define Z ∞ G(x − y)h(y, t)dy, hef f (x, t) := (G ∗ h)(x, t) = (1.18) −∞ where G(z) is a weighting function that describes how the strength of the non-local competition falls off with distance. Standard probability density functions are natural choices for the non-local kernel G, and three of these are most commonly studied in the literature: √ 1. Laplace (Double exponential): G(z) = 2 − e 2σ √ 2|z| σ z2 1 e− 2σ2 2πσ √ if |z| ≤ 3σ 2. Gaussian (Normal): G(z) = √ 3. Uniform: G(z) = √1 2 3σ 0 otherwise where σ is the standard deviation of the distribution [62]. The non-local Fisher equation (1.17) with hef f given by (1.18) with any of kernels 1-3 can have very different dynamics from the “local” system (where hef f (x, t) := h(x, t)). This problem has been studied in the context of biological invasions by Sasaki (1997) and Gourley (2000) [17, 50]. In these studies, the authors found that sufficiently large standard deviation σ for any of these kernels can result in non-monotone or quasi-periodic travelling front solutions, which do not exist for the “local” Fisher equation. The uniform kernel 3, but not 1 or 2, can destabilize the h ≡ 1 state behind the 18 invasion and thereby allow the formation of stationary spatially periodic patterns in the wake of the invasion [17]. Non-local interactions have also been studied using an extension of the non-local Fisher equation that includes an explicit benefit for aggregation as well as competition between individuals for space. This model is given by ∂h ∂2h 2 = + h 1 + αh − βh − (1 + α − β)G ∗ h , ∂t ∂x2 (1.19) where α and β are parameters for the strength of aggregative and competitive forces, respectively [4, 7, 19]. Finally, there is also a substantial literature where “non-local” interactions are incorporated as spatiotemporal convolutions (for instance [2, 8, 20, 21]). The motivation for this form of non-locality is different from that above: the spatiotemporal convolution models a distributed temporal delay in the dynamics and the spatial component arises from the fact that individuals take time to move [8]. However, purely spatial convolutions can be considered in these models as a degenerate case. Since they are specific to different kernel choices and local kinetics, the specific findings of these studies are too numerous to list here. A common theme in these studies of scalar non-local systems, however, is that the non-locality of interactions can facilitate the formation of spatiotemporal patterns, including uniform temporal oscillations, stationary spatially periodic patterns, standing waves and wave trains. Only a few studies have considered the effect of non-local interactions in multi-species systems. For the classical Lotka-Volterra competition model, Britton (1989) extended this to model both intra- and inter-specific competition non-locally [7]: ∂u1 ∂t ∂u2 ∂t = = ∂ 2 u1 ∂x2 ∂ 2 u2 ∂x2 + u1 (1 + α1 u1 − (1 + α1 )G1 ∗ u1 − β1 G2 ∗ u2 ) , + u2 (1 + α2 u2 − (1 + α2 )G2 ∗ u2 − β2 G1 ∗ u1 ) , (1.20) where u1 (x, t) and u2 (x, t) are the population densities of the competing species, α1 , α2 , β1 , β2 are positive parameters and G1 , G2 are (possibly different) non-local competition kernels. In this study, the author finds that non-locality of the interactions promotes aggregation of the popula19 tions (stationary periodic patterns) and can even allow the two species to coexist when local kinetics would predict that they do not (the “principle of competitive exclusion”) [7]. For predator-prey systems, Gourley and Britton (1996) studied a LotkaVolterra predator-prey system where prey dynamics were modelled nonlocally as in (1.19) with β = 0 and a spatiotemporal convolution [18]: ht = Dh ∆h + h (1 + αh − (1 + α)G ∗ ∗h) − hp, pt = ∆p + ap(h − b), (1.21) where h(x, t) and p(x, t) are prey and predator densities, α, a, and b are parameters and G ∗ ∗h is the spatiotemporal convolution with kernel G Z Z t (G ∗ ∗h)(x, t) = G(x − y, t − s)h(y, s)dsdy. Rn −∞ Here the predator-prey interaction is for a type 1 functional response, i.e. mass action, and is purely local. In this study, Gourley and Britton (1996) find that the non-local system (1.21) possess bifurcations of the spatially homogeneous coexistence state that do not occur with purely local prey dynamics. In particular, for purely spatial convolutions (i.e. G(x, t) = R g(x)δ(t), where g ∈ L1 (Rn ) and Rn g(x)dx = 1) there exist bifurcations to stationary spatially periodic, standing wave, and wave train solutions [18]. As in [18], this thesis considers predator-prey systems with non-local intraspecific prey competition, but where the predator functional response is saturating as in the local kinetics (1.3)(i)−(iii). Specifically, I consider here system (1.6) with local kinetics (1.3)(i) and (ii) but use a spatial convolution for the prey competition. These systems are then of the form (i) ht = Dh hxx + h(1 − G ∗ h) − hp h+C , hp pt = Dp pxx − Bp + A h+C (1.22) (ii) ht = Dh hxx + h(1 − G ∗ h) − pt = Dp pxx + Bp 1 − hp hp A h+C , 20 where all variables and parameters are as in (1.3) and G ∗ h is a spatial convolution as in (1.18). Note that unlike (1.21) these systems do not provide an explicit benefit for aggregation. 1.2 Objectives of the thesis In this thesis, I study the generation and stability of spatiotemporal patterns formed in the wake of predator invasions in oscillatory reaction-diffusion predator-prey models and the extensions (1.22) of these models to include non-local prey competition. From the background presented in Section 1.1, it is evident that quite a bit is already known about this topic. Notably, oscillatory reaction-diffusion systems have been shown to support wave train solutions in the wake of invasion, and these wave trains may be unstable and lead to spatiotemporal chaos. There is evidence in natural systems of both of these dynamics, and so predator invasion has been considered as a possible mechanism for the generation of these patterns. There are, however, still a number of open problems on the theoretical study of these systems and this thesis attempts to address three of these problems, described in detail below. While the work presented here furthers our understanding of spatiotemporal pattern formation behind predator invasions, it by no means completes it. Further questions and open problems are discussed in the concluding chapter, Chapter 5, of the thesis. 1.2.1 Wave train selection in the local system (Chapter 2) As described in the background, wave trains generated in the wake of predator invasions in oscillatory reaction-diffusion models are somehow selected from a one-parameter family of possible wave trains. Sherratt (1998) studied this problem using the λ − ω system approximation and developed a criterion that applies when the system is close to the Hopf bifurcation in the kinetics and has nearly equal diffusion coefficients [52]. Away from the Hopf bifurcation, or with highly different diffusivities, this approximation and the resulting selection criterion fail to provide accurate predictions of 21 selected wave trains. It is desirable to have accurate predictions of selected wave trains because each member of this family has different properties (amplitude, speed, spatial and temporal periods, stability) that may have very different ecological and management implications. The development and validation of a more accurate selection criterion is the objective of Chapter 2. In this study, we develop a selection criterion that does not depend specifically on the λ − ω system, but rather on the full reaction-diffusion system and its wave train solutions. This criterion is then tested for reaction-diffusion systems (1.6) with kinetics (1.3)(i) − (iii) with equal diffusion coefficients Dh , Dp as well as some cases with unequal diffusion coefficients. 1.2.2 Spatiotemporal patterns behind invasion in the non-local system (Chapter 3) Non-local interactions have been shown in a wide range of models to destabilize spatially homogeneous states and lead to spatiotemporal patterns. This form of interaction may be more realistic than the local interactions of reaction-diffusion models, and it is therefore important to establish when and how the dynamics for non-local models differ significantly from their local counterparts. This problem has not been considered yet for predator-prey systems that have saturating predator functional responses and corresponding oscillatory local kinetics, or in terms of the spatiotemporal patterns formed behind predator invasions. In Chapter 3, we study the systems (1.22)(i) and (ii) that incorporate non-local intraspecific competition of the prey species, with a particular focus on the spatiotemporal patterns formed behind predator invasions. The objective of this study is to determine if and when non-local competition of the prey will qualitatively change the spatiotemporal patterns that evolve in the wake of predator invasions. In previous studies the emergence of spatiotemporal patterns was shown to depend on the range of the non-local interaction, measured by the standard deviation σ of the kernel G(z), as well as on the exact kernel form. The kernels 1-3 vary greatly in kurtosis, 22 a measure of the kernel “peakedness” and equal to the fourth standardized moment of the distribution. Kernels with the same peakedness as the Gaussian distribution are termed mesokurtic, whereas those with greater kurtosis, such as kernel 1, are termed leptokurtic and are more peaked and have fatter tails. Kernels with smaller kurtosis than the Gaussian kernel, such as kernel 3, are termed platykurtic and have broader humps and thinner tails [31]. Platykurtic kernels, in particular, have been shown to be destabilizing for a range of models (see, for instance [9, 25]). For this study we therefore consider all three of the kernels (Laplace, Gaussian and Uniform) mentioned above, with the aim of determining the sensitivity of the pattern behind invasion to both the standard deviation and kurtosis of the non-local kernel. 1.2.3 Wave train selection and stability in the non-local system (Chapter 4) In the study of Chapter 3, we find that wave trains continue to be a fundamental solution form observed behind predator invasions in systems with non-local prey competition. The theoretical results shown in the introduction, regarding the existence, selection and stability of wave train solutions are only for local kinetics and no such work has been done for wave train solutions of systems (1.22). In Chapter 4 we study the system (1.22)(ii) with non-local competition kernel 1 further. The primary goal of this study is to determine how the non-local competition changes the wave train family and consequently the properties of wave trains observed behind predator invasions. As a prerequisite to this objective, this study also aims to extend the selection criterion developed in Chapter 2 to this non-local system, as well as the methods of Rademacher et al. (2007) for the numerical computation of wave train stability [44]. 23 References [1] D. G. Aronson and H. F. Weinberger. Nonlinear diffusion in population genetics, combustion, and nerve pulse propagation. In A. Dold and B. Eckmann, editors, Partial Differential Equations and Related Topics, volume 446 of Lecture Notes in Mathematics, pages 5–49. SpringerVerlag, 1975. [2] P. Ashwin, M. V. Bartuccelli, T. J. Bridges, and S. A. Gourley. Travelling fronts for the KPP equation with spatio-temporal delay. Z. angew. Math. Phys., 53:103–122, 2002. [3] H. C. Berg. Random Walks in Biology. Princeton University Press, second edition, 1993. [4] J. Billingham. Dynamics of a strongly nonlocal reaction-diffusion population model. Nonlinearity, 17:313–346, 2004. [5] G. Bordiougov and H. Engel. From trigger to phase waves and back again. Physica D, 215:25–37, 2006. [6] P. C. Bressloff and G. Rowlands. Exact travelling wave solutions of an “integrable” discrete reaction-diffusion equation. Physica D, 106:255– 269, 1997. [7] N. F. Britton. Aggregation and the competitive exclusion principle. J. Theor. Biol., 136:57–66, 1989. [8] N. F. Britton. Spatial structures and periodic travelling waves in an integro-differential reaction-diffusion population model. SIAM J. Appl. Math., 50:1663–1688, 1990. [9] M. Doebeli, H. J. Blok, O. Leimar, and U. Dieckmann. Multimodal pattern formation in phenotype distributions of sexual populations. Proc. R. Soc. Lond. B, 274:347–357, 2007. [10] M. Doebeli and T. Killingback. Metapopulation dynamics with quasilocal competition. Theor. Popul. Biol., 64:397–416, 2003. 24 [11] E. J. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, and X. Wang. AUTO97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont). Concordia University, 1997. [12] S. R. Dunbar. Traveling waves in diffusive predator-prey equations: periodic orbits and point-to-periodic heteroclinic orbits. SIAM J. Appl. Math., pages 1057–1078, 1986. [13] L. Edelstein-Keshet. Mathematical Models in Biology. Birkhauser Mathematics Series. McGraw-Hill, first edition, 1988. [14] R. A. Fisher. The wave of advance of an advantageous gene. Annals of Eugenics, 7:355–369, 1937. [15] J. Furter and M. Grinfeld. Local vs. non-local interactions in population dynamics. J. Math. Biol., 27:65–80, 1989. [16] N. J. Gotelli. A Primer of Ecology. Sinauer Associates Inc., second edition, 1998. [17] S. A. Gourley. Travelling front solutions of a nonlocal Fisher equation. J. Math. Biol., 41:272–284, 2000. [18] S. A. Gourley and N. F. Britton. A predator-prey reaction-diffusion system with nonlocal effects. J. Math. Biol., 34:297–333, 1996. [19] S. A. Gourley, M. A. J. Chaplain, and F. A. Davidson. Spatio-temporal pattern formation in a nonlocal reaction-diffusion equation. Dynamical Systems, 16:173–192, 2001. [20] S. A. Gourley and J. W.-H. So. Dynamics of a food-limited population model incorporating nonlocal delays on a finite domain. J. Math. Biol., 44:49–78, 2002. [21] S. A. Gourley, J. W.-H. So, and J. H. Wu. Nonlocality of reactiondiffusion equations induced by delay: Biological modeling and nonlinear dynamics. J. Math. Sci., 124:5119–5153, 2004. 25 [22] K. P. Hadeler and F. Rothe. Travelling fronts in nonlinear diffusion equations. J. Math. Biol., 2:251–263, 1975. [23] A. Hastings, K. Cuddington, K. F. Davies, C. J. Dugaw, S. Elmendorf, A. Freestone, S. Harrison, M. Holland, J. Lambrinos, U. Malvadkar, B. A. Melbourne, K. Moore, C. Taylor, and D. Thomson. The spatial spread of invasions: new developments in theory and evidence. Ecol. Lett., 8:91–101, 2005. [24] R. Hengeveld. Dynamics of biological invasions. Chapman and Hall, 1989. [25] E. Hernandez-Garcia, C. Lopez, S. Pigolotti, and K. H. Andersen. Species competition: coexistence, exclusion and clustering. Phil. Trans. R. Soc. A, 367:3183–3195, 2009. [26] E. E. Holmes, M. A. Lewis, J. E. Banks, and R. R. Veit. Partial differential equations in ecology: Spatial interactions and population dynamics. Ecology, 75:17–29, 1994. [27] J. Huang, G. Lu, and S. Ruan. Existence of traveling wave solutions in a diffusive predator-prey model. J. Math. Biol., 46:132–152, 2003. [28] V. S. Ivlev. Experimental Ecology of the Feeding of Fishes (English translation of Ivlev 1955). Yale University Press, 1961. [29] A. L. Kay and J. A. Sherratt. On the persistence of spatiotemporal oscillations generated by invasion. IMA J. Appl. Math., 63:199–216, 1999. [30] T. H. Keitt, M. A. Lewis, and R. D. Holt. Allee effects, invasion pinning, and species’ borders. Am. Nat., 157:203–216, 2001. [31] R. E. Kirk. Statistics: An Introduction. Holt, Rinehart and Winston, Inc., third edition, 1990. [32] N. Kopell and L. N. Howard. Plane wave solutions to reaction-diffusion equations. Stud. Appl. Math., 52:291–328, 1973. 26 [33] M. Kot, M. A. Lewis, and P. van den Driessche. Dispersal data and the spread of invading organisms. Ecology, 77:2027–2042, 1996. [34] Y. A. Kuznetsov. Elements of Applied Bifurcation Theory, volume 112 of Applied Mathematical Sciences. Springer-Verlag, 1995. [35] S. A. Levin and L. A. Segel. Pattern generation in space and aspect. SIAM Review, 27:45–67, 1985. [36] M. A. Lewis and P. Kareiva. Allee dynamics and the spread of invading organisms. Theor. Popul. Biol., 43:141–158, 1993. [37] M. A. Lewis and S. Pacala. Modeling and analysis of stochastic invasion processes. J. Math. Biol., 41:387–429, 2000. [38] R. M. May. Stability and Complexity in Model Ecosystems, volume 6 of Monographs in population biology. Princeton University Press, second edition, 1974. [39] J. D. Murray. Mathematical Biology I: An Introduction. SpringerVerlag, third edition, 2002. [40] A. Okubo, P. K. Maini, M. H. Williamson, and J. D. Murray. On the spatial spread of the gray squirrel in Britain. Proc. R. Soc. Lond. B, 238:113–125, 1989. [41] M. R. Owen and M. A. Lewis. How predation can slow, stop or reverse a prey invasion. Bull. Math. Biol., 63:655–684, 2001. [42] S. V. Petrovskii and B. L. Li. Exactly solvable models of biological invasion. Mathematical biology and medicine series. Chapman & Hall, 2006. [43] S. V. Petrovskii and H. Malchow. Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics. Theor. Popul. Biol., 59:157–174, 2001. [44] J. D. M. Rademacher, B. Sandstede, and A. Scheel. Computing absolute and essential spectra using continuation. Physica D, 229:166–183, 2007. 27 [45] J. D. M. Rademacher and A. Scheel. Instabilities of wave trains and Turing patterns. Int. J. Bif. Chaos, 17:2679–2691, 2007. [46] M. L. Rosenzweig and R. H. MacArthur. Graphical representation and stability conditions of predator-prey interactions. Am. Nat., 97:209– 223, 1963. [47] B. Sandstede. Stability of travelling waves. In B. Fiedler, editor, Handbook of Dynamical Systems, volume 2, pages 984–1055. Elsevier Science B. V., 2002. [48] B. Sandstede and A. Scheel. Absolute versus convective instability of spiral waves. Phys. Rev. E, 62:7708–7714, 2000. [49] B. Sandstede and A. Sheel. Absolute and convective instabilities of waves on unbounded and large bounded domains. Physica D, 145:233– 277, 2000. [50] A. Sasaki. Clumped distribution by neighborhood competition. J. Theor. Biol., 186:415–430, 1997. [51] J. A. Sherratt. Unstable wavetrains and chaotic wakes in reactiondiffusion systems of λ − ω type. Physica D, 82:165–179, 1995. [52] J. A. Sherratt. Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave. Physica D, 117:145–166, 1998. [53] J. A. Sherratt. Periodic travelling waves in cyclic predator-prey systems. Ecol. Lett., 4:30–37, 2001. [54] J. A. Sherratt, B. T. Eagan, and M. A. Lewis. Oscillations and chaos behind predator-prey invasion: mathematical artifact or ecological reality? Phil. Trans. Roy. Soc. Lond. B, 352:21–38, 1997. [55] J. A. Sherratt, M. A. Lewis, and A. C. Fowler. Ecological chaos in the wake of invasion. Proc. Natl. Acad. Sci., 92:2524–2528, 1995. 28 [56] J. A. Sherratt and M. J. Smith. Periodic travelling waves in cyclic populations: field studies and reaction-diffusion models. J. R. Soc. Interface, 5:483–505, 2008. [57] J. G. Skellam. Random dispersal in theoretical populations. Biometrika, 38:196–218, 1951. [58] M. J. Smith, J. D. M. Rademacher, and J. A. Sherratt. Absolute stability of wavetrains can explain spatiotemporal dynamics in reactiondiffusion systems of lambda-omega type. SIAM J. App. Dyn. Sys., 8:1136–1159, 2009. [59] M. J. Smith and J. A. Sherratt. The effects of unequal diffusion coefficients on periodic travelling waves in oscillatory reaction-diffusion systems. Physica D, 236:90–103, 2007. [60] P. Turchin. Complex Population Dynamics: A Theoretical/Empirical Synthesis. Monographs in Population Biology. Princeton University Press, 2003. [61] W. van Saarloos. Front propagation into unstable states. Phys. Rep., 386:29–222, 2003. [62] D. D. Wackerly, W. Mendenhall, and R. L. Scheaffer. Mathematical Statistics with Applications. Duxbury Press, fifth edition, 1996. [63] H. F. Weinberger, M. A. Lewis, and B. Li. Analysis of linear determinacy for spread in cooperative models. J. Math. Biol., 45:183–218, 2002. [64] S. Wiggins. Global Bifurcations and Chaos: Analytical Methods, volume 73 of Applied Mathematical Sciences. Springer-Verlag, 1988. [65] M. Williamson. Biological invasions, volume 15 of Population and community biology series. Chapman & Hall, 1996. 29 Chapter 2 Wave Train Selection Behind Invasion Fronts in Reaction-Diffusion Predator-Prey Models 1 2.1 Introduction The cause of temporal cycles in natural populations has been a focus of study by ecologists for many decades. A classical hypothesis is that this oscillatory behaviour arises from the interaction between a predator population and its prey, and many models have been constructed and studied to support this hypothesis (see, for example [17]). Such models have often taken the form of kinetics systems: ordinary differential equation models that describe the time evolution of predator and prey densities that are assumed to be spatially constant. More recently, however, field studies have shown that in some natural populations oscillations are not synchronized in space, and when viewed in one spatial dimension take the form of a wave train [2, 14–16, 21, 22]. Wave trains, or periodic travelling waves, are spatiotemporal patterns that are periodic in both time and space and have the appearance of a spatially periodic solution that maintains its shape and moves at a constant speed. Consequently, there has been a great deal of study recently on oscillatory reaction-diffusion systems because these par1 A version of this chapter has been submitted. S. M. Merchant and W. Nagata, Wave train selection behind invasion fronts in reaction-diffusion predator-prey models. 30 tial differential equation models possess wave train solutions (see [28] and references therein). One way that wave trains can arise in oscillatory reaction-diffusion systems is following a predator invasion [26, 27, 29]. The initial condition for such a scenario consists of the prey at carrying capacity everywhere in the spatial domain, except a localized region in which a predator is introduced. Typically, a travelling front evolves that maintains its shape and moves at a constant speed. In some cases, behind this primary invasion front a secondary transition occurs and the solution takes the form of a wave train. Two numerical simulations where wave trains evolve following a predator invasion are illustrated in Figure 2.1. We can see from these examples that the wave train behind the front does not necessarily move at the same speed, or even in the same direction, as the invasion front itself. For oscillatory reaction-diffusion systems near a Hopf bifurcation in the corresponding kinetics there exists a one-parameter family of wave train solutions and a range of corresponding speeds [11]. In a particular numerical simulation of an invasion we typically observe only a single member of this family and this seems robust to changes in initial or boundary conditions. Therefore, it appears that a particular wave train is somehow selected out of the family. We would like find some means of predicting the selected wave train. Sherratt (1998) has, in fact, already produced an explanation of the selection mechanism and a prediction for the wave train selected behind invasion fronts in reaction-diffusion systems with oscillatory kinetics [24]. The basis of his prediction is an approximating lambda-omega (λ-ω) system. The behaviour of an oscillatory reaction-diffusion system near a nondegenerate supercritical Hopf bifurcation can be described by the simpler λ-ω system whose coefficients are obtained from the normal form of the Hopf bifurcation in the kinetics system. Predictions derived in this way are applicable near the Hopf bifurcation and when the predator and prey have nearly equal diffusion coefficients. For more widely applicable predictions, such as in cases where there are larger amplitude oscillations or unequal diffusion coefficients, it would be beneficial to develop a criterion to predict the selected 31 (a) 0.3 time → 0.25 p(x,t) 0.2 0.15 0.1 0.05 0 0 500 1000 1500 2000 x (b) 0.25 time → p(x,t) 0.2 0.15 0.1 0.05 0 0 200 400 600 800 x 1000 1200 1400 1600 Figure 2.1: Wave trains behind invasion fronts. The horizontal axis is the spatial coordinate x. The solution for p(x, t) is plotted as a function of space x at ten fixed times, with vertical separation proportional to the time interval between solutions. In (a) the invasion front first decays to the spatially homogeneous coexistence state, followed by a transition to a wave train further back (Case I). In (b) the invasion front transitions directly to a wave train (Case II). wave train that does not directly depend on the λ-ω system. In the remainder of this chapter we derive and test such a criterion. We first introduce in Section 2.2 the class of two-component reaction-diffusion systems we consider. These systems describe the evolution of population density distributions of two species, one a prey and the other a predator, in one space dimension. Two spatially homogeneous steady states are relevant: an unstable prey-only state that is invaded by a travelling front, and a coexistence state unstable to oscillatory modes that interacts with the invasion. In some cases, such as illustrated in Figure 2.1(a), there is a secondary front that invades the coexistence state. In many cases, the speed of a front invading an unstable steady state can be predicted by the 32 linear spreading speed (see the review [32] and references therein) which depends only on linearization about the unstable state. In Section 2.3 we consider coherent structures in the complex Ginzburg-Landau (CGL) equation [1, 30–33], of which the λ-ω system is a special case. The unstable state in this case is the origin, which corresponds to the coexistence state in predator-prey systems, and coherent structures represent travelling fronts that connect the steady state to wave trains. The linear spreading speed selects a particular coherent structure and wave train, and this retrieves the prediction developed in [24]. Coherent structures have been generalized as defects in general reaction-diffusion systems by Sandstede and Scheel (2004) in [23]. In Section 2.4 we extend the prediction for the λ-ω system to a new “pacemaker” criterion for defects in predator-prey reaction-diffusion systems that connect the unstable prey-only state with wave trains associated with oscillatory instabilities of the coexistence state. For the speed of the selected defect we take the minimum of the linear spreading speeds for the prey-only and coexistence states, and for the frequency of the selected wave train measured in the frame comoving with the defect we take the frequency of the linear Hopf instability of the coexistence state. The performance of the pacemaker criterion is then numerically tested in Section 2.5 on sample oscillatory reaction-diffusion systems. We find that the pacemaker criterion gives accurate predictions for a wider range of parameter values than the λ-ω criterion does, but still falls off in accuracy farther away from the Hopf bifurcation. Finally, Sections 2.6 and 2.7 discuss and summarize the key results. 2.2 Mathematical background We consider predator-prey reaction-diffusion systems in one space dimension, of the form ∂h ∂t ∂p ∂t 2 = Dh ∂∂xh2 + f (h, p), 2 ∂ p = Dp ∂x 2 + g(h, p), (2.1) 33 where h(x, t) is the density of prey at position x and time t and p(x, t) is the density of predator at (x, t). Both h and p are real-valued functions. The positive parameters Dh and Dp are the diffusion coefficients of the prey and predator respectively, while the functions f (h, p) and g(h, p) depend on parameters not explicitly shown here, and describe the local population dynamics. For the invasion scenario of interest, we require (2.1) to have two spatially homogeneous steady states: a prey-only steady state h(x, t) ≡ 1, p(x, t) ≡ 0 and a coexistence steady state h(x, t) ≡ h∗ , p(x, t) ≡ p∗ where both species persist at some non-zero levels. We assume that both the prey-only state (1, 0) and the coexistence state (h∗ , p∗ ) are unstable as fixed points for the corresponding kinetics system dh dt dp dt = f (h, p), = g(h, p). (2.2) In particular we assume that the linearization of (2.2) about the prey-only state has real eigenvalues of opposite sign, while the linearization about the coexistence state has complex conjugate eigenvalues with positive real part and nonzero imaginary parts, and for some nearby parameter values the coexistence state (h∗ , p∗ ) undergoes a supercritical Hopf bifurcation for (2.2). Figure 2.1(a) illustrates an invasion that appears to involve two travelling fronts, a primary front invading the unstable prey-only state, and a secondary front invading the unstable coexistence state. The two fronts do not necessarily travel at the same speed. Figure 2.1(b) shows there may be a single front invading the prey-only state, but we consider only cases where the dynamics are still influenced by the coexistence state. The speed at which fronts invade an unstable steady state has been the subject of much study. A comprehensive review of this topic is provided by [32]. In this, van Saarloos (2003) defines a linear spreading speed ν ∗ given by solving the 34 saddle point equations ∗ ∗ 0 = S(k , ω ) 0 = (∂k + ν ∗ ∂ω )S|(k∗ ,ω∗ ) 0 = =(ω ∗ − ν ∗ k ∗ ), (2.3) for (k ∗ , ω ∗ , ν ∗ ) where S(k, ω) is the characteristic equation for the linearization about the unstable steady state ahead of the front. While there may be multiple solutions to (2.3), only those for which <(D) > 0, where D= −i(∂k + ν ∗ ∂ω )2 S 2∂ω S (k∗ ,ω ∗ ) are relevant. The equivalence of this approach with the historical pinch point analysis is discussed in [7]. When there are several dynamically relevant saddle points we take the one with the largest corresponding ν ∗ to give the linear spreading speed. Details of how to compute linear spreading speeds using (2.3) for the system (2.1) are given in Appendix A. Fronts propagating into unstable states are grouped into two classes: pulled fronts that travel at speed ν ∗ and are in some sense generic, and pushed fronts that travel at a speed ν > ν ∗ . If initial conditions decay sufficiently rapidly in space, faster than eλ ∗x as x → ∞, where λ∗ := =(k ∗ ), it is argued that a front invading the unstable state forms whose asymptotic speed is ν ∗ in the case of a pulled front [32]. In this chapter, we assume that all fronts we encounter are of the pulled variety, and so travel at the speed given by (2.3). This assumption is a crucial component of the criterion developed, and we note that for other systems, for example the Lotka-Volterra competition system, it has been shown to be false [9]. However, whenever possible in our numerical simulations we compute rates of spread and note here that in all cases they do indeed appear to asymptotically approach a pulled front speed. Behind the front invading the prey-only state, simulations of predatorprey reaction-diffusion models exhibit a number of different spatio-temporal patterns, including the spatially homogeneous coexistence state, wave trains, and irregular oscillations that have been identified by some as spatio-temporal 35 chaos [18, 26, 27]. In this chapter we restrict our attention to cases where wave trains appear behind the invasion front. The existence of families of wave trains for reaction-diffusion systems has been established by a number of authors [10, 11, 20]. For solutions that are stationary in a frame moving with speed c̄, (2.1) may be reduced to the system of ordinary differential equations dh dξ dp dξ dy dξ dz dξ = y, = z, = = 1 Dh 1 Dp (−c̄y − f (h, p)) , (2.4) (−c̄z − g(h, p)) , where ξ = x − c̄t. It was shown by [11] that if diffusion coefficients Dh and Dp are sufficiently close and the kinetics system (2.2) has a stable limit cycle near a supercritical Hopf bifurcation, then there exists a one-parameter family of limit cycle solutions for the system (2.4). This family of limit cycles corresponds to a family of wave train solutions for the reaction-diffusion system (2.1) near a Hopf bifurcation. In fact, it is known that families of wave trains exist for generic reaction-diffusion systems, for example with unequal diffusion coefficients [20]. The numerical continuation methods used in this chapter involve continuation along a one-parameter wave train family, and so require the existence of such a family of wave train solutions to (2.1). Therefore, we assume that such a family exists. In many cases, for parameters near the Hopf bifurcation, the rear of the front invading the prey-only state in the reaction-diffusion system (2.1) remains near the coexistence state for a significant distance before the solution transitions by a secondary front to a wave train, as in Figure 2.1(a). In these cases we assume that the speeds of the two fronts are independently determined. Then we focus on the secondary front invading the coexistence state, which connects directly to a wave train. With this motivation, Sherratt (1998) has explained the selection of a wave train from a one-parameter family in the λ-ω system ∂u ∂t ∂v ∂t = = ∂2u ∂x2 ∂2v ∂x2 + λ0 u − ω0 v − (λ1 u + ω1 v)(u2 + v 2 ), + ω0 u + λ0 v + (ω1 u − λ1 v)(u2 + v 2 ), (2.5) 36 which can be used to approximate any two-component reaction-diffusion system near a nondegenerate supercritical Hopf bifurcation, provided the diffusion coefficients are equal, or nearly so. Here the real-valued functions u(x, t), v(x, t) represent perturbations from the coexistence state, so the zero state in (2.5) corresponds to the coexistence state in (2.1). The constants λ0 , ω0 , λ1 and ω1 are the coefficients of the normal form of the Hopf bifurcation in the kinetics system (2.2), and therefore depend on the parameters of the kinetic terms f and g. In particular, the linearization of the kinetic terms about the zero state has eigenvalues λ0 ± iω0 . See [13, 25] for more details. We assume in the following that λ0 > 0, ω0 > 0, λ1 > 0 and ω1 6= 0, thus the zero state is unstable and the Hopf bifurcation in the kinetics system is nondegenerate and supercritical. The advantage of the λ-ω system is that its one-parameter family of wave train solutions can be explicitly written down [11]. These solutions are u(x, t) = r cos(kx − Ωt), v(x, t) = r sin(kx − Ωt), p where Ω = −ω0 − ω1 r2 and k = ± λ0 − λ1 r2 are parameterized by the amplitude r. The wave train with amplitude r therefore has phase speed c̄ = Ω ω0 + ω1 r2 = ±√ . k λ0 − λ1 r 2 (2.6) For the λ-ω system (2.5), Sherratt (1998) predicts that the wave train selected behind a front invading the zero state is the one with amplitude [24] 2λ0 R= ω12 2.3 1 q 2 2 2 λ1 + ω1 − λ1 . (2.7) Coherent structures and selection in the lambda-omega system In this section we use the linear spreading speed of a coherent structure in the λ-ω system to find the wave train selected behind an invading front, 37 and thus recover the prediction (2.7) of [24]. This derivation is implicit in the physics literature cited, but we reproduce it here to motivate the more general selection criterion we use in the case when the λ-ω system is not a good approximation to (2.1). This derivation also points out that the selected wave train only requires that initial conditions decay sufficiently rapidly in space. 2.3.1 Coherent structures The λ-ω system (2.5) is a special case of the CGL equation, well-studied in the physics literature (see, for example, [1]). We are interested in coherent structures that correspond to travelling fronts connecting a steady state to a wave train. The one-dimensional (cubic) CGL equation is At = A + (1 + ic1 )Axx − (1 − ic3 )|A|2 A, (2.8) where A(x, t) is a complex-valued function, and , c1 and c3 are real parameters. To write (2.5) in this form, we put A(x, t) iω0 t u(x, t) + iv(x, t) = √ e λ1 (2.9) in (2.5) and obtain (2.8) with = λ0 , c1 = 0, c3 = ω1 . λ1 A coherent structure of the CGL equation (2.8) is a solution of the form A(x, t) = e−iωt a(ξ)eiφ(ξ) , ξ = x − νt, (2.10) where a and φ are real-valued functions and ω and ν are real parameters. Substituting the ansatz (2.10) into the CGL equation (2.8) and defining q = φ0 , κ= a0 , a 38 we obtain the three-dimensional system of ordinary differential equations a0 = κa, κ0 = −λ0 − νκ + a2 + q 2 − κ2 , q0 = −ω − νq − ω1 2 λ1 a (2.11) − 2κq. The existence and linearized stability of fixed points of the three-dimensional reduced system (2.11) is conveniently summarized in [31, Appendix A]. There are L (“linear”) fixed points (aL , κL , qL ) given by aL = 0, κL + iqL = −ν ± p ν 2 − 4(λ0 + iω) , 2 (2.12) and N (“nonlinear”) fixed points (aN , κN , qN ) with aN = p λ0 − (qN )2 , κN = 0, qN νλ1 1 ± = 2ω1 2 s ν 2 λ21 λ1 ω + 4λ0 . +4 2 ω1 ω1 (2.13) The L fixed points of the reduced system (2.11) correspond to the homogeneous zero steady state solution for both the CGL equation and the λ-ω system, and the N fixed points correspond to wave train solutions of the CGL equations A(x, t) = aN ei[qN x−(ω+qN ν)t] , with constant amplitude aN and phase speed ν + ω/qN . For the λ-ω system p the corresponding wave trains, using (2.9), have amplitude r = aN / λ1 and p phase speed c̄ = ν + (ω − ω0 )/qN = ±(ω0 + ω1 r2 )/ λ0 − λ1 r2 . We are interested in heteroclinic orbits of the reduced system (2.11), connecting an L fixed point and an N fixed point. These heteroclinic orbits correspond to coherent structures or travelling fronts of the CGL equation (2.8), connecting the zero solution and a wave train. Due to symmetry, we consider only the case ν > 0, for a front travelling to the right. Correspondingly, in the reduced system we seek heteroclinic orbits which approach an N fixed point as ξ → −∞, and approach an L fixed point as ξ → +∞. If such orbits exist, they must leave along the unstable manifold of the N fixed 39 point and enter along the stable manifold of the L fixed point. Hence, the dimensions of the local stable and unstable manifolds of the L and N fixed points indicate whether the heteroclinic orbits that we seek are at least possible, and whether they can be expected to be robust. Inspecting the eigenvalues of the fixed points reveals that there is an N fixed point with a p one-dimensional local unstable manifold, and if ν > |ω|/ λ0 > 0 then there is an L fixed point with a three-dimensional local stable manifold. (If ω = 0 then (2.11) restricted to the invariant plane a = 0 is orbitally equivalent to a Hamiltonian system, and in this case there is an L fixed point with a p three-dimensional local stable manifold provided ν ≥ 2 λ0 .) Hence, there is the possibility of a two-parameter family of heteroclinic orbits, parameterized by ν and ω. We have not proved that such a two-parameter family of heteroclinic orbits exists, but for various values of ν and ω we have numerically computed heteroclinic orbits of (2.11) using the continuation software package auto and the process described in [5, 8, 12]. The computations suggest that a two-parameter family of heteroclinic orbits of the reduced system indeed exists. 2.3.2 Wave train selection As we have just described, for given coefficients λ0 , ω0 , λ1 and ω1 there can be a continuous family of heteroclinic orbits of the reduced system (2.11), parameterized by ν and ω. In a particular simulation of the λ-ω system we see only one of the corresponding front solutions. That is, particular values of ν and ω appear to be selected. Moreover, the selection seems robust to changes in initial conditions and boundary conditions. We predict the selected front and corresponding wave train in the λ-ω system, by using the linear spreading speed of a pulled front at the zero state A(x, t) ≡ 0 of the CGL equation. Application of (2.3) to the reduced system (2.11) at the relevant L fixed point gives p κ∗L = − λ0 qL∗ = 0, 40 and the linear spreading speed and corresponding frequency are p ν ∗ = 2 λ0 , ω ∗ = 0. (2.14) For the parameter values (2.14), a heteroclinic orbit for (2.11) connects the L ∗ fixed point (0, κ∗L , qL∗ ) with the N fixed point (a∗N , 0, qN ), where from (2.13) and (2.14) we have √ a∗N = q λ0 − ∗ )2 , (qN ∗ qN λ0 =− ω1 q 2 2 λ1 + ω1 − λ1 (2.15) (only one of the two possible qN is valid, since (aN )2 > 0). Many authors have observed that fronts in the CGL equation indeed appear to move at the linear spreading speed, but we have no proof that for all initial conditions that decay sufficiently rapidly in space, the selected front must be as predicted, a pulled front moving at the linear spreading speed given by (2.14). For a few cases, we have compared numerically computed heteroclinic orbits of the reduced system to the front solutions produced by direct simulation of the λ-ω system. One example is given in Figure 2.2. In p this figure, we have shown the amplitude u2 + v 2 of the numerical solution p of the λ-ω system (2.5), along with the amplitude a/ λ1 of the fronts corresponding to the heteroclinic orbits of the reduced system (2.11) computed using auto. In all the cases we studied, we found close agreement between the predicted and observed coherent front solutions. From other simulations we observe that this agreement also appears to be robust with respect to perturbations of the initial condition. With the parameters given by (2.14), the front solution for the λ-ω system (2.5) is a coherent structure a(ξ) u(x, t) + iv(x, t) = √ eiω0 t eiφ(ξ) , λ1 ξ = x − ν ∗ t, with frequency ω0 and moving to the right with speed ν ∗ = 2 p λ0 . As ξ → +∞ the front solution approaches the homogeneous zero steady state, 41 0.07 0.06 amplitude 0.05 0.04 0.03 0.02 0.01 0 1000 1020 1040 1060 1080 1100 1120 1140 x Figure 2.2: Heteroclinic orbits of (2.11) compared p with simulations of the system (2.5). The solid curve is the amplitude u2 + v 2 from numerical simulation of (2.5). Dashed curves are the (spatially translated) amplitude p a/ λ1 of front solutions corresponding to heteroclinic orbits of (2.11) computed in auto for ω = ω ∗ and ν = ν ∗ + ∆ν with (from right to left) ∆ν = 0.15, 0.1, 0.05, 0.005. Other parameter values are λ0 = 0.05, ω0 = 0.14142, λ1 = 9.25926 and ω1 = −14.32219. and as ξ → −∞ it approaches the wave train solution a∗ ∗ u(x, t) + iv(x, t) = √N eiω0 t eiqN ξ . λ1 (2.16) To see the relation with the prediction (2.7) from [24], we note that a wave train solution of the λ-ω system (2.5) of the form u(x, t) + iv(x, t) = rei(kx−Ωt) under the change of variables ξ = x − ν ∗ t becomes rei(ν ∗ k−Ω)t eikξ . Comparing with (2.16), we see that in the comoving frame the selected wave 42 train has temporal frequency ν ∗ k − Ω = ω0 , (2.17) and spatial wavenumber ∗ k = qN . Substitution of this wavenumber using (2.15) to find the amplitude of the p selected wave train r = a∗N / λ1 retrieves the prediction (2.7) derived by [24]. We have verified this selection with numerical simulations of the λ-ω system using randomly-generated initial conditions. See Figure 2.4(i). We note in (2.17) that the frequency of the selected wave train in the frame comoving with the selected coherent structure coincides with the imaginary part of the complex eigenvalue λ0 + iω0 of the kinetics system for (2.5) linearized around the unstable zero state. We can therefore think of the linear Hopf frequency ω0 at the zero state as a “pacemaker”, in the sense p that the selected wave train in the frame moving with speed ν ∗ = 2 λ0 must have the same temporal frequency ω0 . 2.4 Beyond the lambda-omega system: the pacemaker criterion In this section we consider wave train selection behind invading fronts in the predator-prey reaction-diffusion systems described in Section 2.2. In general for the full oscillatory reaction-diffusion systems (2.1), wave trains are not sinusoidal and we do not have exact solutions for them, so a prediction of the form (2.7) is not possible. However, as shown in Figure 2.1, numerical simulations of predator invasions in the full system even well away from the Hopf bifurcation in the kinetics still appear to have travelling fronts connecting steady states to wave trains. For the full system we assume there exist two wave train selection regimes, which we refer to here as Case I and Case II. In Case I, suggested by Figure 2.1(a), we assume there is a primary front 43 invading the prey-only state with speed cinv that connects to the coexistence state that is stationary in the frame comoving with speed cinv , and a secondary front invading the coexistence state with speed csec that connects to a wave train and is temporally periodic in the frame comoving with speed csec . We assume the two fronts can be treated independently, and the relevant front for selection is the secondary one connecting the coexistence state to the wave train. In this case selection is directly analogous to selection in the λ-ω system. In Case II we assume there is a front with speed cinv , temporally periodic in the frame comoving with speed cinv , that connects the prey-only state directly with the selected wave train, as suggested by Figure 2.1(b). In addition, we assume that there exists, at least nearby, a front connecting the coexistence state to the wave train. In both cases we assume the fronts are pulled and so the speeds cinv and csec are the linear spreading speeds obtained by linearization about the spatially homogeneous prey-only and coexistence steady states respectively. We also assume that there do not exist wave trains of the same speed cinv as the primary invasion front, for in this case the invasion may instead take place through a pointto-periodic heteroclinic connection in a single frame moving with speed cinv , as studied by [6]. In both Cases I and II, we propose that the speed ccoh of the front that selects the wave train is the minimum of the two speeds defined above, ccoh := min{csec , cinv }, and we assume the front moves to the right: ccoh > 0 (for fronts moving to the left, replace x with −x). Our justification for this value of the selecting front speed is as follows. First, if csec < cinv then we predict that the tail of the primary invasion front decays to the coexistence state (h∗ , p∗ ) and behind the primary front a secondary transition occurs from the coexistence state via a front with speed ccoh = csec . This corresponds to Case I of the criterion. If, on the other hand csec > cinv , such a secondary front could not exist for long since it would catch up to the primary front. In this case, we therefore predict that eventually the primary invasion front moving with 44 speed ccoh = cinv will not decay to the coexistence steady state and instead connect directly from the prey-only state to the selected wave train, which is Case II. In analogy with the selection criterion (2.17) for the CGL equation we conjecture that in the comoving frame of speed ccoh the selected wave train has the frequency ωcoh of the linear unstable oscillatory mode of the coexistence state. To express this selection criterion in terms of the wave train parameters, let L > 0 be the spatial period, cwtrain be the phase speed and T = L/cwtrain be the (signed) temporal period of the selected wave train. Then the frequency of the wave train in the comoving frame of speed ccoh is (2π/L)ccoh − (2π/T ). Our criterion for the selected wave train is that it satisfies 2π 2π ccoh − = ωcoh , L T (2.18) where the “pacemaker” frequency ωcoh is the imaginary part of the eigenvalue of the linearization of the kinetics system about the coexistence state: ωcoh 1 = −sgn(ω0 ω1 ) 2 s ∂f ∂g 2 ∂f ∂g − − −4 ∂h ∂p ∂p ∂h (2.19) (h∗ ,p∗ ). where ω0 and ω1 are two of the normal form coefficients of the Hopf bifurcation in the corresponding kinetic system. We note that since here the selection criterion depends on these kinetic normal form coefficients we are assuming that the system is sufficiently close to a supercritical Hopf bifurcation in the kinetics. If f (h, p) and g(h, p) are given by the kinetics of the λ-ω system, the selection criterion (2.18) is the same as the selection criterion (2.17) using the p linear spreading speed for the λ-ω system, with ωcoh = ±ω0 , ccoh = 2 λ0 , L = 2π/|k| and T = ±2π/Ω, and so we have in some sense extended the prediction for the λ-ω system to one that may be used on the full system. We acknowledge that in Case II, as in Figure 2.1(b), the coexistence state (h∗ , p∗ ) is not approached and so there is no good reason to expect that the pacemaker frequency ωcoh plays a role as the frequency in the selection crite- 45 rion. Nevertheless, it performs surprisingly well as a predictor in numerical simulations and so we continue to use it for the selection criterion. While the motivation for the criterion has a fairly sound mathematical basis, the criterion itself is a somewhat naive extension of those ideas. Essentially, we are assuming that a generalized temporally periodic “coherent front” solution exists that connects a steady state (either prey-only or coexistence) to the selected wave train solution. Furthermore, we are assuming that this coherent front is of the pulled variety and, critically, we are assuming that we know the period 2π/ωcoh of this coherent front. The criterion (2.18) therefore has the drawbacks that these assumptions may be false. However, the criterion has a number of advantages. Most obviously, it may be applied to systems other than the λ-ω system. In addition, although the motivation for the criterion depends on the diffusion coefficients of the two species being equal, the criterion itself does not require this. 2.5 The criterion in practice To study the validity of the criterion (2.18) we have considered three particular forms for the kinetic equations in the oscillatory reaction-diffusion system (2.1) as well as a λ-ω system. That is, we consider the four sets of kinetics (i) f (h, p) = λ0 h − ω0 p − (λ1 h + ω1 p)(h2 + p2 ), g(h, p) = ω0 h + λ0 p + (ω1 h − λ1 (iii) p)(h2 f (h, p) = h(1 − h) − p 1 − e−Ch , g(h, p) = Bp A − 1 − Ae−Ch + (ii) p2 ) g(h, p) = −Bp + (iv) hp , h+C hp A h+C f (h, p) = h(1 − h) − hp f (h, p) = h(1 − h) − A h+C , p g(h, p) = Bp 1 − h where ω0 , ω1 , λ0 , λ1 , A, B and C are parameters. Models (ii)-(iv) all assume a logistic growth for the prey species in absence of the predator and assume a saturating functional response for the predation term, but differ in how this saturation and resulting predator growth are modelled [17]. In addition, (ii)-(iv) have all been utilized in the past to model predator invasions and the production of wave trains following the primary invasion front [26, 27]. Our 46 goal here is to both illustrate how one would apply the selection criterion (2.18), as well as to study the accuracy of the criterion for these well-known kinetics and compare with the prediction that would be arrived at using the λ-ω system. 2.5.1 Equal diffusion coefficients For each model we set Dh = Dp = 1 and chose a bifurcation parameter (B for models (ii) and (iv) and C for model (iii)) and computed the normal form parameters for the kinetics. The bifurcation parameter is then linearly related to the parameter λ0 as shown in Table 2.1, with the supercritical Hopf bifurcation occurring at λ0 = 0. For computational convenience, for all models the remaining parameters were fixed at values that resulted in the same values of the corresponding normal form parameters ω0 and c3 = ω1 /λ1 . Hence, for a given value of the parameter λ0 , all models have the same predictions for selected speeds of wave trains using the λ-ω system. The values of these various parameters are provided in Table 2.1. Table 2.1: Parameter values for models (i)-(iv) Model A B C ω0 λ1 ω1 λ0 (i) 0.14142 9.25926 -14.32219 - (ii) 0.15 0.2 0.14142 9.25926 -14.32219 8(0.1 − B) (iii) 3.02173 0.02923 0.14142 0.21101 -0.32638 0.03684(C − 2.54504) (iv) 6.92368 0.63542 0.14142 18.35223 -28.38715 0.02268 − B For each of the models (i)-(iv) we computed the predicted selected wave train using (2.18) by numerical continuation with the software package auto. We continued the system (2.4) in the kinetic bifurcation parameter λ0 as a boundary value problem with periodic boundary conditions subject to the constraint (2.18). In tandem with this, we also continue the linear spreading speed equations as outlined in Appendix A for the speeds csec and cinv that correspond to Cases I and II of the criterion respectively. We note 47 that since the secondary front speed goes to zero at the Hopf bifurcation, Case I of the criterion applies near the Hopf bifurcation at λ0 = 0 and Case II may apply farther away from the bifurcation point, with the exchange between the two cases occurring if these two curves cross. We briefly mention here how we determined initial solutions for the continuation of the selection criterion. Two different means were used to find an initial wave train that satisfies (2.18). Case I of the criterion always applies close to the Hopf bifurcation in the kinetics and so in this case we start at λ0 = 0.0001 and use the analytical criterion (2.7) from the λ-ω system to find an approximate wave train speed to start. We then continue in the speed parameter the system (2.4) from the Hopf bifurcation to a wave train of the desired speed. For Case II of the criterion we started the continuation at λ0 = 0.005 and so the λ-ω system may not be a good approximation. In this case, we first generated a series of wave trains of different speeds by continuation in c̄ of system (2.4). For each of these wave trains we computed the temporal frequency (2π/L)cinv − (2π/T ) in the cinv frame and the wave train with temporal frequency nearest to ωcoh was used as the starting point for the continuation. The speeds of these predicted wave trains found by continuation are shown as curves in Figure 2.4. In this figure solid curves correspond to criterion (2.18) with ccoh = csec , dashed curves use the criterion with ccoh = cinv . In each case shown here the criterion with the proposed speed ccoh = min{csec , cinv } corresponds to the uppermost of these two curves. To compare with these predictions, we performed numerical simulations of each system over a range of the bifurcation parameter. Simulations were performed using a finite-difference discretization of the system (forward differences in time and centered differences in space) in the comoving frame of speed −ccoh (computed as above) on the interval [−l, 0] with boundary conditions h(−l, t) = 1, p(−l, t) = 0 ∂h ∂p (0, t) = 0, (0, t) = 0 ∂ξ ∂ξ 48 for models (ii)-(iv), and for model (i) ∂h ∂p (0, t) = 0, (0, t) = 0. ∂ξ ∂ξ h(−l, t) = 0, p(−l, t) = 0 where ξ = x + ccoh t. We note that our numerical simulations are of leftward moving fronts for computational convenience and so we perform the coordinate transformation ξ 7→ −ξ for later comparison with predictions made from the criterion (2.18). The domain length l was chosen so that the primary invasion front remained far from the left boundary for the duration of the simulation; we set 2000 + T (c − c ), c < c sec sec inv end inv l= 2000, csec ≥ cinv where Tend is the final time for the numerical simulation. The initial condition used for simulations of models (ii)-(iv) was h∗ , h(ξ, 0) = 1, ξ ≥ −1000 ξ < −1000 , p∗ , p(ξ, 0) = 0, ξ ≥ −1000 ξ < −1000 and for model (i) we initialized the system at the (h, p) = (0, 0) state everywhere with the addition of white noise of magnitude = 0.001 on the subinterval [−1000, 0]. Most simulations were run up to a final time Tend = 10000, with a few cases run to Tend = 20000 for the model (iv) because the evolution of the selected wave train in these cases was so slow that no wave train was observed at time t = 10000. A snapshot of the spatial distribution of the predator and prey populations was saved at times t = 981, 982, 983, · · · , 1000, t = 1981, 1982, · · · , 2000 etc. for the computation of the “observed” speed of the wave train selected by the primary front. Similarly to the reaction-diffusion systems, linear spreading speeds may also be computed for finite difference systems such as our discretization here and in general differ from the continuum system (for example, see [7]). Since 49 the linear spreading speed ccoh is a key component of our selection criterion, we chose the spatial mesh size ∆ξ and the timestep ∆t so that the linear spreading speed for the discretized system was within 0.5% of the continuum value. This required different choices for ∆ξ and ∆t in different parameter regions for each of the models, which we summarize Table 2.2. Table 2.2: Mesh interval and timestep for numerical simulations Dh = 1, Dp ≤ 1 Dh = 1, Dp = 2 (i) 1 4 1 ∆t = 50 [0.02,0.295] Range (ii) [0.023,0.199] in λ0 (iii) [0.03,0.195] [0.005,0.025] - [0.03,0.19] [0.005,0.01] (iv) [0.014,0.022] [0.007,0.0135] [0.002,0.0065] - - ∆ξ = 1 8 1 ∆t = 200 [0.005,0.015] 1 16 1 ∆t = 800 - 1 4 1 ∆t = 50 - 1 8 1 ∆t = 200 - [0.007,0.021] [0.003,0.005] [0.025,0.195] [0.005,0.02] ξ= ∆ξ = ∆ξ = ∆ξ = For each set of parameter values, the spatial profiles were plotted at times t = 1000, 2000, 3000, · · · and visually inspected. For the majority of cases, a wave train did appear to evolve behind the primary invasion front and maintained its form over time. For some parameter sets of models (ii) and (iii) in the Case I parameter region, however, a selected wave train was only observed transiently. In these cases, there does initially appear to be a secondary “front” of speed csec that selects a wave train, but after some finite time the secondary transition moves faster than the predicted speed csec , as illustrated in Figure 2.3, and the region behind the invasion became irregular for later times. For the speeds of wave trains in the simulations, the final set of 20 profiles was used to compute this speed if the wave train remained regular throughout the time-series. If, however, the wave train in the final data set appeared irregular as in Figure 2.3, then we instead used the last set of profiles prior to the change in speed of the secondary transition. In all cases the speed of the wave train was computed as cwtrain = d/19 − ccoh where d is the distance travelled over the course of the 19 time units by the second 50 2.5 time → p(ξ,t) 2 1.5 1 0.5 0 -4500 -4000 -3500 -3000 -2500 -2000 -1500 -1000 -500 0 ξ Figure 2.3: Example of a transiently observed selected wave train. The solution for p(ξ, t) is plotted at times t = 1000, 2000, · · · , 10000 as a function of ξ, with vertical separation proportional to the time interval between solutions. From simulation of model (iii) with parameters as in Table 2.1 and λ0 = 0.015, Dh = 1, Dp = 1. peak of the wave train behind the primary invasion front. The distance d was numerically computed by a simple C routine that tracked locations of local maxima chronologically through the set of profiles in the time-series. These “observed” speeds for the various models are shown in comparison with the predicted speeds in Figure 2.4. We can see from the results shown in Figure 2.4 that the performance of the criterion (2.18) is dependent on the particular model considered, as well as how far the parameter values are from the Hopf bifurcation in the kinetics. We first note that for the λ-ω model (i), the selection criterion developed in [24] and also here using the CGL equation performs extremely well, providing evidence that this selection mechanism is indeed valid at least for the λ-ω system. We can also observe from the data shown in Figure 2.4 that for all three models (ii)-(iv) the criterion appears to provide a higher order prediction than the λ-ω system prediction near the Hopf bifurcation and so retains accuracy farther from it. For models (ii) and (iii), we eventually see the loss of accuracy for the pacemaker criterion as the parameter values are increased still farther from the bifurcation point. In the case of model (ii) this occurs after there is an apparent switch from Case 51 Kinetics (i) with Dh = 1, Dp = 1 1 0 cwtrain -1 -2 -3 -4 -5 0 0.05 0.1 0.15 λ0 0.2 0.25 0.3 Kinetics (ii) with Dh = 1, Dp = 1 0 0 -1 -0.5 -2 cwtrain -3 -4 -1 -5 0 0.01 0.02 0.03 0.04 0.05 -1.5 -2 0 0.05 0.1 λ0 0.15 0.2 Figure 2.4: Wave train speeds in numerical simulations compared with predictions for cases with equal diffusion coefficients. Horizontal axis is the bifurcation parameter λ0 . Vertical axis are the corresponding wave train speeds cwtrain . Predictions using the λ-ω criterion (2.7) are dot-dashed lines. Predictions using the new pacemaker criterion (2.18) are solid lines for Case I and dashed lines for Case II. Simulation observations are point symbols. Solid circles are speeds computed with final time t = 10000. Squares are speeds computed with final time t = 20000. Open circles are cases where the selected wave train was only observed transiently and speeds were then computed as described in the text. Inset of (ii) shows the region near the kinetic Hopf bifurcation. 52 Kinetics (iii) with Dh = 1, Dp = 1 0 -0.5 -1 cwtrain -1.5 -2 -2.5 -3 -3.5 -4 0 0.05 0.1 λ0 0.15 0.2 Kinetics (iv) with Dh = 1, Dp = 1 0 -1 -2 cwtrain -3 -4 -5 -6 -7 -8 -9 0 0.005 0.01 λ0 0.015 0.02 Figure 2.4: continued. 53 I to Case II of the criterion whereas for model (iii) the observed speeds fall away from the prediction while still in the Case I region. In contrast, the pacemaker criterion performs extremely well for the entire parameter range chosen for model (iv). 2.5.2 Unequal diffusion coefficients We also conducted a limited study for models (i) and (iv) with Dh = 1, Dp = 0.5 and models (ii) and (iii) with Dh = 1, Dp = 2. All other parameters were taken as in Table 2.1 and the procedure used is as described for equal diffusion coefficients above. In this case, initial solutions for the continuation of the selection criterion are obtained from the equal diffusion cases by first continuing to the appropriate value of the diffusion coefficient parameter Dp . With unequal diffusion coefficients the analytical prediction (2.7) breaks down, and so we illustrate that in this case the criterion (2.18) provides an improved prediction even for the kinetics (i). We also note that for model (ii) we only show a small kinetic parameter range since beyond the range shown wave trains of speed cinv exist (and are observed in the simulations), in contradiction to our assumption for the application of the criterion. The results for the numerical simulations and the predicted speeds using (2.18) are shown in Figure 2.5. The results of this study for the case of unequal diffusion coefficients show that, at least for these examples, this criterion performs quite well even in the case of unequal diffusion coefficients. This suggests some promise for this criterion to apply to the more general and biologically relevant case of unequal predator and prey diffusion coefficients. 2.6 Discussion The pacemaker criterion (2.18) developed here provides a relatively simple method for predicting the selection of wave trains following predator invasions in oscillatory reaction-diffusion models. Our study of the performance of the pacemaker criterion for models (ii)-(iv) does suggest that this crite- 54 Kinetics (i) with Dh = 1, Dp = 0.5 0 -0.5 cwtrain -1 -1.5 -2 -2.5 -3 0 0.02 0.04 0.06 0.08 0.1 λ0 Kinetics (ii) with Dh = 1, Dp = 2 -1 -1.5 -2 cwtrain -2.5 -3 -3.5 -4 -4.5 -5 0 0.02 0.04 0.06 λ0 0.08 0.1 0.12 Figure 2.5: Wave train speeds in numerical simulations compared with predictions for cases with unequal diffusion coefficients. Horizontal axis is the bifurcation parameter λ0 . Vertical axis are the corresponding wave train speeds cwtrain . Curves and point symbols are as described for Figure 2.4. For model (i), the analytical prediction (2.7) assuming Dh = 1 = Dp is shown (dot-dashed curve) for comparison. 55 Kinetics (iii) with Dh = 1, Dp = 2 0 -0.5 -1 cwtrain -1.5 -2 -2.5 -3 -3.5 -4 0 0.05 0.1 λ0 0.15 0.2 Kinetics (iv) with Dh = 1, Dp = 0.5 0 -1 cwtrain -2 -3 -4 -5 -6 -7 0 0.005 0.01 λ0 0.015 0.02 Figure 2.5: continued. 56 rion may in general be valid for some range of parameters near the Hopf bifurcation in the kinetics. However, the performance of the pacemaker criterion is clearly model-dependent, and indeed since for each of the models chosen we studied only a small slice of parameter space it is possible that it would also be parameter-dependent. In fact, we have no a priori means of predicting when the pacemaker criterion will provide an accurate prediction. We can, however, identify a number of factors that may be responsible for deviation from the predictions of the pacemaker conjecture and suggest areas for future work. Our development of the criterion assumes that there exists a front between a spatially homogeneous steady state (whether prey-only or coexistence) and the selected wave train. While visual inspection of numerical simulations suggests that such fronts do exist, we have no mathematical proof of this. Fronts of the type we assume to exist are called defects in [23]; they are solutions that are temporally periodic in some comoving frame and asymptotic in space to wave trains (or steady states). Defects are heteroclinic connections in an infinite-dimensional space of temporally periodic functions, and what is required is a proof that these heteroclinic connections indeed exist in the reaction-diffusion systems we consider. We also assume that this heteroclinic connection has the particular frequency ωcoh , given by the imaginary part of the eigenvalues of the linearization of the kinetics at the coexistence steady state. We chose this frequency because it extends the λ-ω system prediction and seems to perform well for the criterion. For Case I this is natural, and it is somewhat surprising that it can perform well even for Case II where the primary front does not approach the coexistence state. Clearly, if our assumptions about the existence and frequency of the heteroclinic connection are false then we cannot expect the pacemaker criterion to be accurate. It is therefore important to establish the existence of such heteroclinic connections for the full model under study, rather than just the λ-ω system. The coherent front solutions for the CGL equation are a special case of the more general “coherent pattern forming fronts” discussed in [32] and there may be a method to compute such coherent pattern forming fronts numerically by continuation in auto using the techniques described 57 2.5 time → p(x,t) 2 1.5 1 0.5 0 -4000 -3500 -3000 -2500 -2000 x -1500 -1000 -500 0 Figure 2.6: Example of an unstable selected wave train. The solution for p(x, t) is plotted at times t = 1000, 2000, · · · , 10000 as a function of space x, with vertical separation proportional to the time interval between solutions. From simulation of model (ii) with parameters as in Table 2.1 and λ0 = 0.04, Dh = 1, Dp = 1. in [4] and references therein. While the numerical computation of these heteroclinic connections would not rigorously prove their existence, it would offer additional evidence for this proposed wave train selection mechanism and could open a path to further numerical work on the multiplicity and stability of such structures. In this study we have made no reference to the stability of the wave trains and heteroclinic connections considered. In fact, in our simulations we observe two potential instabilities. The first is as shown in Figure 2.3 where a front of (approximately) the predicted speed is only observed transiently. One possible reason for this transient behaviour is that the heteroclinic connection is unstable. Since our numerical simulations are for finite times, it is in fact possible that even the simulations that remained regular (closed circles in Figs. 2.4 and 2.5) are transients. The study of the stability of the heteroclinic connections assumed by the selection criterion would therefore be an important area of future study. A second form of instability that we see in our simulations is shown in Figure 2.6. This figure is shown in the stationary frame for illustration. In this case, a selected wave train is observed behind the initial predator 58 invasion front, but becomes irregular further back, a solution form that has been noted for predator invasions in a number of studies and attributed to the instability of the wave train [28]. The physics literature has long noted that wave trains may be either convectively or absolutely unstable [3]. When a solution is convectively unstable in a given frame of reference, localized perturbations will grow, but will be convected away quickly enough that at a fixed position in the frame of reference the perturbation in fact decays. If the solution is absolutely unstable in the frame of reference, however, perturbations grow at any fixed position. In our simulations in the comoving frame ccoh , irregularities such as those shown in Figure 2.6 were convected away through the boundary and so these selected wave trains are presumably only convectively unstable. It would be interesting to determine whether this is generally the case or whether absolutely unstable wave trains may also be selected. We expect that methods developed in [19] would be useful to numerically study this problem. We have used these methods to compute the essential spectra of wave trains in the predator-prey models for a few cases and found the selected wave trains to be unstable. However, we have been unable to compute the absolute spectra of non-sinusoidal wave trains. Finally, as we saw in our study of models with unequal predator and prey diffusion coefficients the developed criterion appears to be useful for prediction of wave train selection when the predator and prey species do not diffuse with the same strength. Since in natural systems species likely have different movement rates, this would be a substantial advance. However, we have only shown here a few cases and so the application of the criterion to unequal diffusion coefficients in general requires further study. Here again rigourous mathematical work on the existence of the assumed heteroclinic connections should be useful. 2.7 Conclusion For oscillatory reaction-diffusion equations to be of practical use in applied ecology, they must provide testable predictions for the behaviour of the 59 natural system they are attempting to model. Numerical simulation of such equations reveals that for some parameter regions wave train solutions evolve behind predator invasion fronts. Moreover, although mathematical theory states that a continuous family of wave trains is possible, for each set of kinetic parameters typically a single wave train is selected and the selection appears to be independent of initial conditions and boundary conditions. Using the λ-ω approximation, Sherratt (1998) provided an understanding of the mechanism and a prediction for the selected wave train that is valid near a supercritical Hopf bifurcation in the kinetics [24]. However, as we have demonstrated, this prediction rapidly loses accuracy as the system parameters are moved away from the Hopf bifurcation point. We have therefore reconsidered the selection criterion from [24] from the viewpoint of coherent structures in the complex Ginzburg-Landau equation, which has led us to formulate a “pacemaker” selection criterion which generalizes that of [24] but using the parameters and wave train solutions of the full system, rather than the corresponding λ-ω system approximation. Testing this criterion on three sample oscillatory predator-prey reaction-diffusion systems has shown that, for these systems at least, the pacemaker criterion provides an improved approximation for the selected wave train near the Hopf bifurcation in the kinetics. The pacemaker criterion remains accurate for a much larger region of parameter space than the λ-ω system criterion, but does eventually lose accuracy as well. In addition, the new criterion provided an improved prediction in examples for which the diffusivities of the predator and prey differed. Mathematically rigourous studies of the existence of the fronts involved in selection would likely show in what sense the pacemaker criterion gives a higher order approximation than one using solely the λ-ω system. The new criterion moves us a step closer to understanding the process by which wave trains are formed following predator invasions, but further mathematical analysis is still required to fully understand the selection mechanism. 60 Acknowledgements This work was partially supported by funding from the Natural Sciences and Engineering Research Council of Canada (NSERC). As well, S.M. would like to acknowledge the Pacific Institute for the Mathematical Sciences (PIMS) for funding through the IGTC fellowship program. We would also like to thank Michael Doebeli and his lab group for many helpful comments and suggestions. Finally, we are grateful for the comments and suggestions of two anonymous referees as we think they helped greatly improve the manuscript. 2.8 Appendix A - Linear spreading speeds In this appendix we describe how to compute the linear spreading speeds using (2.3) for the oscillatory reaction-diffusion system (2.1). We assume the exponential ansatz h, p ∼ ei(kx−ωt) for perturbations of the linearized system. This gives the characteristic equation S(k, ω) = −ω 2 + iω −k 2 (Dh + Dp ) + a1,1 + a2,2 + Dh Dp k 4 −Dp a1,1 k 2 − Dh a2,2 k 2 + a1,1 a2,2 − a1,2 a2,1 , where " a= ∂f ∂h ∂g ∂h ∂f ∂p ∂g ∂p (2.20) # (h0 ,p0 ) is the linearization matrix about the steady state (h0 , p0 ). Therefore, for the system (2.1) the linear spreading speed equations (2.3) are (dropping 61 the stars) 0 = −ω 2 + iω −(Dh + Dp )k 2 + a1,1 + a2,2 + Dh Dp k 4 − Dp a1,1 k 2 (2.21) 2 −Dh a2,2 k + a1,1 a2,2 − a1,2 a2,1 , 0 = −2iω(Dh + Dp )k + 4Dh Dp k 3 − 2Dh a2,2 k − 2Dp a1,1 k − 2νω (2.22) −iν(Dh + Dp )k 2 + iν(a1,1 + a2,2 ), 0 = =(ω − νk). (2.23) Typically, for given values of a, Dh and Dp the system of equations (2.21)(2.23) and hence the linear spreading speed ν ∗ must be solved for numerically. Indeed, for our use in this study we numerically continued solutions of (2.21)-(2.23) as 5 real equations in tandem with the wave train selection criterion. In addition, we track the quantity D= −ω(Dh + Dp ) − 6iDh Dp k 2 + iDh a2,2 + iDp a1,1 − 2ν(Dh + Dp )k + iν 2 −2ω − i(Dh + Dp )k 2 + ia1,1 + ia2,2 to ensure that we are at a dynamically relevant saddle point. For the special case Dh = Dp = 1 considered for models (i)-(iv) in this study, we can solve for the linear spreading speeds cinv and csec analytically, as we show in the following section. We use these analytical results as the initial solution for the continuation of (2.21)-(2.23). For the unequal diffusion cases in this study we begin from the analytical prediction for Dh = Dp = 1 and first continue in the diffusion coefficient Dp . 2.8.1 Linear spreading speeds for the equal diffusion cases When Dh = Dp = 1 we have from (2.22) that either ν = −2ik or ω = −ik 2 + i(a1,1 + a2,2 )/2. The second case is in general inconsistent with (2.21) and so we assume the first case, ie. that ν = −2ik. In this case, k must be purely imaginary and taking k = iγ, γ ∈ R in (2.23) gives =(ω) = 2γ 2 . Hence, we take ω = α + 2iγ 2 , α ∈ R in (2.21) and get the following two real 62 equations for α and γ. 0 = γ 4 + γ 2 (−a1,1 − a2,2 ) − α2 + a1,1 a2,2 − a1,2 a2,1 , 0 = α(2γ 2 − a1,1 − a2,2 ). (2.24) Computing csec for Dh = Dp = 1: This is the speed obtained by linearizing about the spatially homogeneous coexistence steady state (h0 , p0 ) = (h∗ , p∗ ). The second equation in (2.24) has the possible solution α = 0. In this case, solution of the first equation for γ using α = 0 yields γ2 = a1,1 + a2,2 1 ± 2 2 q (a1,1 − a2,2 )2 + 4a1,2 a2,1 . We assume for the selection criterion that ωcoh as defined in (2.19) is real and so the expression on the right has non-zero imaginary part. We therefore assume α 6= 0 and so from the second equation of (2.24) we have γ 2 = (a1,1 + a2,2 )/2. Then using ν ∗ = 2γ we have the linear spreading speed csec = q 2(a1,1 + a2,2 ) for rightward moving fronts emanating from the coexistence steady state. Computing cinv for Dh = Dp = 1: This is the linear spreading speed for fronts emanating from the spatially homogeneous prey-only steady state (h0 , p0 ) = (1, 0). In this case, for models (ii)-(iv) we have a1,1 < 0, a1,2 < 0, a2,1 = 0 and a2,2 > 0. Therefore, if α 6= 0 the first equation of (2.24) gives 1 α2 = − (a1,1 + a2,2 )2 + a1,1 a2,2 , 4 which contradicts α being real. We therefore take α = 0 and solve to get γ 2 = a2,2 . We thus have the linear spreading speed √ cinv = 2 a2,2 63 for rightward moving fronts emanating from the prey-only steady state. 64 References [1] I. S. Aranson and L. Kramer. The world of the complex GinzburgLandau equation. Rev. Mod. Phys., 74:99–143, 2002. [2] O. N. Bjørnstad, M. Peltonen, A. M. Liebhold, and W. Baltensweiler. Waves of larch budmoth outbreaks in the European Alps. Science, 298:1020–1023, 2002. [3] R. J. Briggs. Electron-stream interaction with plasmas. MIT Press, 1964. [4] A. R. Champneys and B. Sandstede. Numerical computation of coherent structures. In B. Krauskopf, H. M. Osinga, and J. Galan-Vioque, editors, Numerical Continuation Methods for Dynamical Systems, pages 331–358. Springer, 2007. [5] E. J. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, and X. Wang. AUTO97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont). Concordia University, 1997. [6] S. R. Dunbar. Traveling waves in diffusive predator-prey equations: periodic orbits and point-to-periodic heteroclinic orbits. SIAM J. Appl. Math., pages 1057–1078, 1986. [7] U. Ebert and W. van Saarloos. Front propagation into unstable states: universal algebraic convergence towards uniformly translating pulled fronts. Physica D, pages 1–99, 2000. [8] M. J. Friedman and E. J. Doedel. Numerical computation and continuation of invariant manifolds connecting fixed points. SIAM J. Numer. Anal., 28:789–808, 1991. [9] Y. Hosono. The minimal speed of traveling fronts for a diffusive Lotka Volterra competition model. Bull. Math. Biol., 60:435–448, 1998. 65 [10] J. Huang, G. Lu, and S. Ruan. Existence of traveling wave solutions in a diffusive predator-prey model. J. Math. Biol., 46:132–152, 2003. [11] N. Kopell and L. N. Howard. Plane wave solutions to reaction-diffusion equations. Stud. Appl. Math., 52:291–328, 1973. [12] B. Krauskopf, H. M. Osinga, E. J. Doedel, M. E. Henderson, J. Guckenheimer, A. Vladimirsky, M. Dellnitz, and O. Junge. A survey of methods for computing (un)stable manifolds of vector fields. Int. J. Bifurcat. Chaos Appl. Sci. Eng., 15:763–791, 2005. [13] Y. A. Kuznetsov. Elements of Applied Bifurcation Theory, volume 112 of Applied Mathematical Sciences. Springer-Verlag, 1995. [14] X. Lambin, D. A. Elston, S. J. Petty, and J. L. MacKinnon. Spatial asynchrony and periodic travelling waves in cyclic populations of field voles. Proc. R. Soc. Lond., Ser. B: Biol. Sci., 265:1491–1496, 1998. [15] J. L. MacKinnon, S. J. Petty, D. A. Elston, C. J. Thomas, T. N. Sherratt, and X. Lambin. Scale invariant spatio-temporal patterns of field vole density. Ecology, 70:101–111, 2001. [16] R. Moss, D. A. Elston, and A. Watson. Spatial asynchrony and demographic traveling waves during red grouse population cycles. Ecology, 81:981–989, 2000. [17] J. D. Murray. Mathematical Biology I: An Introduction. SpringerVerlag, third edition, 2002. [18] S. V. Petrovskii and H. Malchow. Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics. Theor. Popul. Biol., 59:157–174, 2001. [19] J. D. M. Rademacher, B. Sandstede, and A. Scheel. Computing absolute and essential spectra using continuation. Physica D, 229:166–183, 2007. [20] J. D. M. Rademacher and A. Scheel. Instabilities of wave trains and turing patterns. Int. J. Bif. Chaos, 17:2679–2691, 2007. 66 [21] E. Ranta and V. Kaitala. Travelling waves in vole population dynamics. Nature, 390:456, 1997. [22] E. Ranta, V. Kaitala, and J. Lindstrom. Dynamics of Canadian lynx populations in space and time. Ecography, 20:454–460, 1997. [23] B. Sandstede and A. Scheel. Defects in oscillatory media: toward a classification. SIAM J. Appl. Dyn. Syst., 3:1–68, 2004. [24] J. A. Sherratt. Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave. Physica D, 117:145–166, 1998. [25] J. A. Sherratt. Periodic travelling waves in cyclic predator-prey systems. Ecol. Lett., 4:30–37, 2001. [26] J. A. Sherratt, B. T. Eagan, and M. A. Lewis. Oscillations and chaos behind predator-prey invasion: mathematical artifact or ecological reality? Phil. Trans. Roy. Soc. Lond. B, 352:21–38, 1997. [27] J. A. Sherratt, M. A. Lewis, and A. C. Fowler. Ecological chaos in the wake of invasion. Proc. Natl. Acad. Sci., 92:2524–2528, 1995. [28] J. A. Sherratt and M. J. Smith. Periodic travelling waves in cyclic populations: field studies and reaction-diffusion models. J. R. Soc. Interface, 5:483–505, 2008. [29] M. J. Smith and J. A. Sherratt. The effects of unequal diffusion coefficients on periodic travelling waves in oscillatory reaction-diffusion systems. Physica D, 236:90–103, 2007. [30] M. van Hecke. Coherent and incoherent structures in systems described by the 1D CGLE: experiments and identification. Physica D, 174:134– 151, 2003. [31] M. van Hecke, C. Storm, and W. van Saarlos. Sources, sinks and wavenumber selection in coupled CGL equations and experimental im- 67 plications for counter-propagating wave systems. Physica D, 134:1–47, 1999. [32] W. van Saarloos. Front propagation into unstable states. Phys. Rep., 386:29–222, 2003. [33] W. van Saarloos and P. C. Hohenberg. Fronts, pulses, sources and sinks in generalized complex Ginzburg-Landau equations. Physica D, 56:303–367, 1992. 68 Chapter 3 Instabilities and Spatiotemporal Patterns Behind Predator Invasions in Systems with Non-local Prey Competition 2 3.1 Introduction Reaction-diffusion systems are a classical tool for the modelling of ecological systems. In one spatial dimension, such systems take the form ut = Duxx + f (u), x∈R (3.1) where N is the number of species, u is an N -component function of position x and time t that represents the population densities of the species, D is a positive diagonal matrix of diffusion coefficients and the components of f (u), referred to as the “kinetics,” describe the local dynamics of the interacting species. Models of this form have a long history in the study of biological invasion, where a novel species is introduced to a localized region of the spatial domain and spreads across the landscape over time [15, 17]. Perhaps the most 2 A version of this chapter will be submitted. S. M. Merchant and W. Nagata, Instabilities and spatiotemporal patterns behind predator invasions in systems with non-local prey competition. 69 famous example of this is the Fisher-KPP equation ht = Dhxx + h(1 − h), (3.2) for the population density h(x, t) of a single species, with diffusion coefficient D and logistic growth f (h) = h(1 − h) [7]. For this model, a localized non-zero initial density produces travelling waves of invaders, and has been successfully used to describe the range expansion of a number of natural species [17]. In reaction-diffusion systems (3.1), however, the kinetics f (u) are “local”, in that they do not depend directly on space or time, but only on the population densities u(x, t). Growth and death rates at a position x are therefore determined only by the population densities at this same spatial location. This may be unreasonable for some systems, for instance when individuals compete for a common (rapidly equilibrated) resource [8]. Consequently, more recently a series of studies have considered models with so-called “non-local” interactions. Perhaps the most basic example of such non-local interactions is the “non-local Fisher equation” that incorporates non-local intraspecific competition into the standard Fisher-KPP equation (3.2). It has the form ht = hxx + h(1 − G ∗ h), x∈R (3.3) where h(x, t) ∈ R is the population density at x in time t and G ∗ h is the spatial convolution of this density with a weighting function G(z) that describes how the strength of the non-local competition falls off with distance, i.e. Z ∞ G(x − y)h(y, t)dy. (G ∗ h)(x, t) = −∞ Often such weighting functions G are taken to be standard probability density functions. The NLF equation has been studied by Sasaki (1997) and Gourley (2000) in the context of a single species invasion [9, 28]. In these studies it was found that both a Gaussian and an exponential (also known as Laplace) non-local kernel G could result in non-monotone or quasi-periodic 70 travelling front solutions. In addition, a uniform kernel has the potential to destabilize the h ≡ 1 homogeneous state and may allow stationary spatially periodic patterns to evolve in the wake of the invasion [9]. There have also been a number of studies on the equation ht = hxx + h 1 + αh − βh2 − (1 + α − β)G ∗ h (3.4) that can include a benefit of aggregation (the αh term) and/or competition for space (the −βh2 term) [2, 3, 12]. There is also a substantial literature studying models such as these where the non-local term is a spatiotemporal convolution, and so includes a distributed temporal delay for the interaction [1, 4, 13, 14]. While we will be considering here only purely spatial convolutions, we note that a common theme from these studies is that nonlocal competition in such single-species models facilitates the formation of spatiotemporal patterns, including uniform temporal oscillations, stationary spatially periodic patterns, standing waves and wave trains. Compared to the scalar case, to our knowledge relatively few studies have considered the effect of non-local interactions in multi-species systems. Britton (1989) considers the classical Lotka-Volterra competition model where the competition occurs in a non-local manner [3]. He finds that not only does the non-locality promote aggregation of the populations, but that this pattern formation can allow for coexistence of the two species. Gourley and Britton (1996) study the Lotka-Volterra predator-prey system with type I functional response and prey dynamics as in (3.4) with β = 0 and a spatiotemporal convolution [10]. Here they find that even a purely spatial convolution can lead to new spatiotemporal patterns including stationary spatially periodic, standing wave, and wave train solutions. In light of the significant effects of including non-local prey competition in the predator-prey system above, we consider here the influence of such non-local competition in two “oscillatory” predator-prey models commonly 71 used to model predator invasions. These have the general form ht = Dh hxx + f (h, p), pt = Dp pxx + g(h, p), (3.5) where here h(x, t) and p(x, t) are the prey and predator densities respectively, and Dh and Dp are their respective diffusion coefficients. Since a key interest for ecologists is to predict the pattern that will be left behind after a predator invasion has passed, our focus here is on the patterns that develop behind predator invasion fronts in these systems. Due to the oscillatory kinetics, the local versions of these systems are known to possess complex spatiotemporal solutions, including wave trains and spatiotemporal chaos in the wake of invasion [31, 32]. We show here that, even in this case, the introduction of non-local competition can generate new spatiotemporal patterns, including stationary spatially periodic states, and generally increases the propensity for the system to form spatiotemporal patterns behind the invasion. Since the bifurcations found in previous studies have depended on the form of the non-local kernel, we consider here three commonly used kernels (Laplace, Gaussian and uniform) as well as a caricature system that includes the standard deviation and kurtosis of the kernel as parameters. We study the dependence of patterns observed on the extent of the non-locality, as measured by the standard deviation (width) and kurtosis, a measure of the “peakedness” and fatness of tails, of the kernel. We show with this “survey” of kernels and the two models considered that non-local competition may be expected to play the most significant role in pattern generation behind predator invasions when the competition kernel has a large standard deviation and a low kurtosis (i.e. broad peak and thin tails). In contrast, relatively narrow kernels and highly kurtotic ones had very little effect on the spatiotemporal patterns observed behind the invasion. In the next section we present the model and kernels considered, as well as outline the invasion dynamics and spatiotemporal patterns observed in the local model. In Section 3.3 we perform a linear stability analysis of the spatially homogeneous prey-only and coexistence state. We focus on 72 the instabilities of the coexistence state as predictors of the patterns to evolve behind predator invasions and also consider these instabilities in a related caricature system that depends only on the standard deviation and kurtosis of the non-local kernel. Section 3.4 then compares the coexistence state instabilities with the type of spatiotemporal pattern observed behind predator invasions in numerical simulations. We outline the basic trends observed and their implications as well as propose future work needed in this area in Section 3.5. 3.2 The model In this study we consider the two predator-prey systems with non-local prey competition given by (i) ht = Dh hxx + h(1 − G ∗ h) − hp h+C , hp pt = Dp pxx − Bp + A h+C hp (ii) ht = Dh hxx + h(1 − G ∗ h) − A h+C , p pt = Dp pxx + Bp 1 − h where A, B, C > 0 are biological parameters. For both models we consider only the case with equal diffusion coefficients Dh = Dp = 1. We also consider the parameters A and C as fixed, and vary the parameter B as a bifurcation parameter. Each of systems (i) and (ii) is a clear extension of the corresponding “local” model, obtained by taking G(x) = δ(x), with reaction kinetics hp h+C , hp A h+C (i)0 f (h, p) = h(1 − h) − g(h, p) = −Bp + hp (ii)0 f (h, p) = h(1 − h) − A h+C , p g(h, p) = Bp 1 − h . Both models incorporate logistic prey growth and a saturating functional response of the predator [23]. As well, both have been used to model predator invasions and the spatiotemporal patterns that arise behind them (see for instance [31] and [32]). We note that for suitable parameters both kinetics 73 support a prey-only steady state (h, p) ≡ (1, 0) with the prey at carrying capacity and a coexistence state (h, p) ≡ (h∗ , p∗ ) where predator and prey coexist at intermediate levels. The local versions of models (i) and (ii) are so-called “oscillatory” predator-prey models because, for suitable choices of A and C, in the kinetics system the coexistence state undergoes a supercritical Hopf bifurcation as the parameter B is decreased through a critical value B ∗ . For the reaction-diffusion systems with local kinetics this Hopf bifurcation gives the existence of wave trains, also known as periodic travelling waves, beyond the bifurcation point. Many studies have observed such wave train solutions behind predator invasions in numerical simulations of these models. Indeed, we find in our numerical simulations of the local systems that initial conditions of the form h(x, 0) = 1, p∗ , |x| ≤ 50 p(x, 0) = 0, |x| > 50 (3.6) evolved to a predator invasion front behind which we observed the spatially homogeneous state if B > B ∗ , as shown in Figure 3.1(a). For B < B ∗ , a wave train solution was observed at least transiently behind the invasion. Examples of these wave train solutions are shown in Figures 3.1(b)-(d). These wave trains arise in the invasion wake either from the spatially homogeneous state, as in 3.1(b), or directly behind the primary front, as in 3.1(c) and (d). In addition, in some cases these wave trains are not stable and either transitioned to some other wave train solution, as in Figure 3.1(d), or to irregular oscillations as in 3.1(c). As mentioned above, we consider three different kernel forms (Laplace, Gaussian and uniform) for the non-local prey competition. The first two of these are members of the generalized Gaussian distribution, and the last is a limiting case of this distribution [24]. The probability density function for the generalized Gaussian distribution with mean 0, standard deviation σ and shape parameter α is given by x α G(x) = ae−| b | 74 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 time → p(x,t) (a) B = 0.1075 0 200 400 600 800 1000 x time → p(x,t) (b) B = 0.0975 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 200 400 600 800 1000 x time → p(x,t) (c) B = 0.0825 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 200 400 600 x 800 1000 1200 800 1000 1200 time → p(x,t) (d) B = 0.0075 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 200 400 600 x Figure 3.1: Spatiotemporal patterns observed behind predator invasions in the local system. The solution for p(x, t) is plotted at times t = 2000, 2005, · · · , 2045 as a function of space x, with vertical separation proportional to the time interval between solutions. Only the right side of the domain is shown. Simulations are of model (i)0 with Dh = Dp = 1, A = 0.15, C = 0.2 and B as shown. Patterns shown are (a) homogeneous coexistence state (b) decay to homogeneous coexistence, followed by transition to wave train (c) wave train immediately behind front and (d) wave train immediately behind front, followed by transition to a different wave 75 train where 1 a= , 2Γ(1 + 1/α)b s b=σ ∞ Z and Γ(z) is the Gamma function: Γ(z) = Γ(1/α) , Γ(3/α) tz−1 e−t dt. The shape param- 0 eter α determines the kurtosis of the distribution, with α = 2 corresponding to the mesokurtic Gaussian distribution. For values of α < 2 the distribution is leptokurtic (more peaked and fatter-tailed than a Gaussian distribution) and for α > 2 it is platykurtic (broader hump and thinner tails) [18]. We consider here kernels with α = 1, α = 2 and the limiting case α → ∞. In these cases we have the well-known kernels: √ 2 − 1. Laplace (Exponential): α = 1, G(x) = e 2σ 2. Gaussian (Normal): α = 2, G(x) = √ 3. Uniform: α → ∞, G(x) = √1 2 3σ 0 √ 2|x| σ x2 1 e− 2σ2 2πσ if |x| ≤ √ 3σ otherwise and their probability density functions are shown in Figure 3.2. Finally, we note that in all that follows, we have fixed the parameters A and C at A = 0.15, C = 0.20 for model (i), and A = 6.92, C = 0.64 for model (ii). In this case the kinetics systems (i)0 and (ii)0 undergo a nondegenerate supercritical Hopf bifurcation as the parameter B is decreased through the critical value B ∗ , with B ∗ = 0.1 for (i), and B ∗ = 0.02268 for (ii). 3.3 Instabilities of the spatially homogeneous steady states Here we look at the linear stability of the spatially homogeneous steady states of models (i) and (ii) to spatial perturbations of the form ∼ eikx , 76 0.9 1 2 3 0.8 0.7 G(x) 0.6 0.5 0.4 0.3 0.2 0.1 0 -3 -2 -1 0 x 1 2 3 Figure 3.2: Probability density functions G(x) for non-local competition kernels. Kernels shown here are for σ = 1.0 (1. Laplace, 2. Gaussian and 3. Uniform). k ∈ R. To begin we transform the steady state (h0 , p0 ) to the origin and linearize. That is, we take in systems (i) and (ii) with Dh = Dp = 1 h = h0 + h̄, p = p0 + p̄, and discard higher order terms in the perturbations (h̄, p̄) as well as use the fact that (h0 , p0 ) is a spatially homogeneous solution of the system. We then have the linear system for the perturbations h̄t = h̄xx + a1,1 h̄ + h0 h̄ − h0 G ∗ h̄ + a1,2 p̄, p̄t = p̄xx + a2,1 h̄ + a2,2 p̄, where a= a1,1 a1,2 a2,1 a2,2 ! = ∂f ∂h ∂g ∂h ∂f ∂p ∂g ∂p ! (h0 ,p0 ) is the linearization matrix for the local system about the steady state and so is computed from kinetics (i)0 and (ii)0 . To test stability we take the initial conditions to be plane wave perturbations of the steady state and 77 seek solutions of the form h̄(x, t) = Φ1 (t)eikx , p̄(x, t) = Φ2 (t)eikx . Noting that we then have Z ∞ G ∗ h̄ = G(x − y)h̄(y)dy = Φ1 (t)e ikx Z ∞ G(z)e−ikz dz −∞ −∞ ∂Φ1 ikx ∂Φ2 ikx e , p̄t = e , h̄xx = −k 2 Φ1 (t)eikx and p̄xx = ∂t ∂t and substituting this into the linearized system results in a and that h̄t = −k 2 Φ2 (t)eikx linear system for the time-dependent coefficients ∂Φ1 ∂t ∂Φ2 ∂t ! = Z −k 2 + a1,1 + h0 − h0 H(k) a1,2 a2,1 −k 2 + a2,2 ∞ where H(k) = ! ! Φ1 (t) Φ2 (t) , G(z)e−ikz dz is the Fourier transform of the non-local −∞ competition kernel G. In particular, for the three kernels considered here we have 1. Laplace: H(k) = 1 1+ 1 2 2 2σ k σ 2 k2 2. Gaussian: H(k) = e− 2 √ sin(√ 3kσ) 3kσ 3. Uniform: H(k) = 1 if k 6= 0 otherwise. The stability of the homogeneous solution (h0 , p0 ) to perturbations of wavenumber k is then determined by the eigenvalues of the linearization matrix above. For each spatially homogeneous steady state considered here we have then two curves of eigenvalues λ1 (k) and λ2 (k), given by −b(k) ± λ1,2 (k) = 2 p ∆(k) , 78 where b(k) = 2k 2 − a1,1 − a2,2 − h0 + h0 H(k), ∆(k) = (a1,1 + h0 − h0 H(k) − a2,2 )2 + 4a1,2 a2,1 . Together, these curves make up the essential spectrum of the linearization about the steady state. The steady state is linearly unstable if there are any modes k for which the essential spectrum crosses into the right half of the complex plane, that is if <(λ1 (k)) > 0 or <(λ2 (k)) > 0 for any k ∈ R. In the following, we refer to the rightmost (up to complex conjugacy) element of the essential spectrum λ∗ as the leading eigenvalue, and the correponding mode k ∗ as the leading mode. We are interested in determining parameter values at which the spatially homogeneous states lose stability via this rightmost element λ∗ crossing the imaginary axis. We therefore consider the following three types of instability: 1. Hopf: λ∗ = iω ω ∈ R, k ∗ = 0 2. Turing: λ∗ = 0, k ∗ 6= 0 3. Turing-Hopf: λ∗ = iω ω ∈ R, k ∗ 6= 0. 3.3.1 Prey-only steady state For the prey-only homogeneous steady state we have a1,1 = −1 and a2,1 = 0 for both sets of kinetics and so the essential spectrum of the linearized operator is given by λ1 (k) = −k 2 − H(k), λ2 (k) = −k 2 + a2,2 . In the parameter region that allows for predator invasion we necessarily have a2,2 > 0 and so we can see from the second branch of eigenvalues that this steady state is always unstable, as in the local system. We note however, that the first branch of eigenvalues λ1 (k) is in fact the essential spectrum for the h ≡ 1 state of the non-local Fisher equation (3.3). This state therefore loses 79 stability when the non-local kernel has a Fourier transform H(k) for which −H(k) > k 2 for some wavenumber k. Since H(k) > 0 for the Laplace and Gaussian kernels this is clearly not possible. However, for the uniform kernel this state can be destabilized for sufficiently large standard deviation (with onset at σ ≈ 5.297795). Since the leading wavenumber is non-zero and the corresponding eigenvalue λ∗ is real this instability is of Turing type, and may lead to pattern formation in the prey only system, as noted in [9]. In Section 3.3.3 we study this instability further for comparison with the caricature system developed there, and we do find that perturbation of the h ≡ 1 state in numerical simulations leads to stationary spatially periodic patterns. We note that our primary interest, however, is on predator invasions into spatially homogeneous prey distributions and so when necessary we restrict to σ < 5.29 for the uniform kernel. 3.3.2 Coexistence steady state The stability of the spatially homogeneous coexistence state is dependent on both the parameter B and the standard deviation σ of the non-local competition kernel, and all three types of instability (Hopf, Turing and Turing-Hopf) can occur for this state. While we have not been able to compute the stability boundaries for this state analytically, we have numerically computed them by continuation with the software auto by solving for local maxima of <(λ) at which <(λ) = 0 [6]. These stability boundaries in the B − σ plane are shown as curves in Figure 3.9 for the Laplace non-local kernel, in Figure 3.10 for the Gaussian kernel and in Figure 3.11 for the uniform kernel. Although the two models behave quite differently, we note some general trends in the location of the stability boundaries. For both models and all three kernels, small values of σ have the Hopf bifurcation in the kinetics as the primary instability. This is in agreement with the local system behaviour. As well, for all kernels as the value of σ is increased, the parameter range of B for which the coexistence state is linearly stable shrinks, and so longer range non-local competition has a greater destabilizing effect. Finally, this trend is also true for the non-local kernels ordered 1-3. That is, for both 80 models kernel 1 (Laplace) has a larger stable parameter space than kernel 2 (Gaussian) and kernel 2 has a larger stable region than 3 (uniform). This suggests that the shape of the kernel has a significant effect on the stability of the homogeneous coexistence state. Since these kernels differ significantly in kurtosis (with decreasing kurtosis for kernels 1-3) we examine the dependence of the stability boundaries on kurtosis in a caricature model in the next section. In terms of the nature of the primary instability, we note that the models (i) and (ii) differ somewhat. As σ is increased, the primary instability transitions from Hopf to Turing-Hopf for both models with the Laplace kernel. For the Gaussian kernel, however, both models transition from Hopf to Turing-Hopf instability, but model (i) has an additional transition to Turing instability for even larger σ. For the uniform kernel, model (ii) again transitions from Hopf to Turing-Hopf instability whereas model (i) has only a transition from Hopf to Turing instability as σ is increased. The nature of the primary instability is therefore dependent on both the form of the non-local kernel and on the predator-prey local kinetics. Model (i) with the Gaussian or uniform kernels (2 and 3) possesses a Turing instability which may lead to the bifurcation of stationary spatially periodic patterns. Such patterns do not exist for any of the parameter range considered here for the the local model. To test whether such stationary patterns bifurcate we chose a parameter range adjacent to the Turing instability boundary. We then performed numerical simulations, using the methods described in Section 3.4, of perturbations of the coexistence steady state starting from the initial conditions h(x, 0) = h∗ + 0.0001 cos(k ∗ x), p(x, 0) = p∗ + 0.0001 sin(k ∗ x) on [−250, 250] and at the coexistence steady state elsewhere, where here k ∗ is the fastest growing mode from the linear stability analysis. At each time step, we computed the L2 -norm of the difference, δ(t), between the perturbed and unperturbed states on [−250, 250]. The rate of growth of the perturbation was then taken as the slope of a linear fit to log(δ(t)) from t = 81 15 to t = 30 (chosen to avoid initial phase adjustment and saturation effects at later times). From the time-series, simulations appeared to approach a stationary pattern by time t = 100, and so at this time we computed the wavelength of the spatial pattern, measured as the distance between successive troughs in the prey density. The growth rate of perturbations and corresponding wavelengths of spatial patterns are shown in Figure 3.3(a) and (b). As we can see from this figure, for both kernels the Turing instability does indeed appear to lead to supercritical bifurcation of stationary spatial patterns of approximately the predicted wavelength. (a) (b) 0.5 24 22 spatial period rate of growth 0.4 0.3 0.2 20 18 16 14 0.1 12 0 3 4 5 6 7 σ 8 9 3 4 5 6 7 8 9 σ Figure 3.3: Growth rates of perturbations and spatial period of pattern for coexistence steady state. Growth rate vs. σ is shown in panel (a), wavelength vs. σ in panel (b) for the system (i) with B = 0.105 and uniform (blue open circles and dashed curves) and Gaussian (red filled circles and solid curves) kernels. Observed values from numerical simulations are point symbols, predictions using linearized stability are shown as curves. 3.3.3 A caricature system and the effect of kurtosis In this section we study the effect of the standard deviation and kurtosis of the competition kernel on the stability and spatial patterns formed from the spatially homogeneous steady states in a related “caricature” of the predator-prey system. First, we derive the caricature system to show the 82 relationship to the non-local system. ForZ the convolution term in systems ∞ G(z)h(x − z, t)dz and expand (i) and (ii), we note that (G ∗ h)(x, t) = −∞ h(x − z, t) in a Taylor series about x (as in [19] and [23], section 11.5) hxx (x, t) 2 hxxx (x, t) 3 hxxxx (x, t) 4 h(x−z, t) ≈ h(x, t)−hx (x, t)z+ z − z + z −· · · 2 3! 4! The convolution above can therefore be written as (G ∗ h)(x, t) = ∞ X (−1)j j=0 Z j! Gj ∂j h (x, t), ∂xj ∞ where Gj = z j G(z)dz is the j th central moment of the distribution with −∞ density function G(z). Gj → 0 as j → ∞ for most probability distributions, We expect that j! so we truncate this to the 4th moment. In addition, we note that G0 = 1 for all probability density functions and we assume symmetric probability distributions, so that the odd moments are zero. Finally, the second central moment of a distribution is the variance, i.e. G2 = σ 2 and the fourth central moment is proportional to the kurtosis in that G4 = Kσ 4 , where K is the kurtosis of the distribution. In summary, we take an approximation to the convolution of the form (G ∗ h)(x, t) ≈ h(x, t) + σ2 Kσ 4 hxx (x, t) + hxxxx (x, t). 2 4! We choose here to truncate at the first two moments because we are interested in exploring how the standard deviation and kurtosis of the distribution affect the potential for spatiotemporal pattern formation, and because statistical techniques exist to estimate these parameters (σ and K) in natural systems. This approximation of G ∗ h gives us a caricature system for each of models (i) and (ii) ht = hxx − σ2 2 hhxx − Kσ 4 4! hhxxxx + f (h, p), pt = pxx + g(h, p). We now study the spatially homogeneous steady states in this reduced 83 system with the aim of determining how their stability depends on the standard deviation and kurtosis of the non-local prey competition kernel. In addition, this caricature can allow us to make some predictions for how the wavelength of spatial patterns will depend on these properties of the kernels. Linearizing the caricature system about the spatially homogeneous state (h0 , p0 ) as for the full system gives us again two branches of eigenvalues for perturbations of the form ∼ eikx . In this case these are given by −b(k) ± λ1,2 (k) = 2 p ∆(k) , (3.7) where Kσ 4 h0 4 σ2 k + 2 − h0 k 2 − a1,1 − a2,2 , 4! 2 4 2 Kσ h0 4 σ 2 h0 2 ∆(k) = a1,1 − k − k − a2,2 + 4a1,2 a2,1 . 4! 2 b(k) = The prey-only state As for the full system, we focus here on the carrying capacity steady state h ≡ 1 in the prey-only system. In this case for each perturbation wavenumber k we have a single corresponding real eigenvalue, that we write here as a quadratic in q = k 2 , λ(q) = −1 − q + Kσ 4 2 σ2 q− q , 2 4! 6 2 σ − 2 . Since q = k 2 Kσ 4 we note that σ 2 > 2 is necessary for this to exist. This steady state is then giving us a parabola with maximum at q ∗ = linearly stable in the prey-only system provided λ(q ∗ ) < 0 and unstable if λ(q ∗ ) > 0. In terms of the standard deviation and kurtosis, we therefore 2 have 3 (σ − 2)2 that the steady state is linearly stable provided that K > 2 σ4 and unstable otherwise. Since for this caricature system we have an expression for the most unp stable mode k ∗ = ± q ∗ of the linearized system, we can see that this most 84 unstable mode is a decreasing function of K in the relevant parameter regions. Therefore, we expect that if stationary spatially periodic patterns arise from the prey-only state (in the absence of predator introduction), the wavelength of the spatial pattern will increase as the kurtosis of the non-local competition kernel is increased. In terms of the standard deviation, since σ 2 > 2 is assumed we see that for fixed kurtosis the wavenumber of the most unstable mode increases with σ up to a maximum at σ = 2, and a decreasing function of the standard deviation for σ > 2. We therefore expect, based on the caricature system, that spatial patterns formed from the prey-only state will have a wavelength that decreases with standard deviation down to a minimum, and increases thereafter. For the full system, the only kernel that had an instability of this prey only state is the uniform competition kernel, and in this case we compare the predictions of the caricature system with that of the full kernel. The kurtosis of the uniform kernel is K = 9/5. For the caricature system the spatially homogeneous state at this value of kurtosis is in fact linearly stable for all standard deviations σ. However, for now we ignore this and still consider the leading mode as a predictor of the wavelength of spatial patterns. In Figure 3.4 we plot the predicted and observed wavelengths of spatial patterns using the full uniform kernel, along with the predictions from the caricature system, against the standard deviation σ of the kernel. From this figure we can see that knowledge of only these first two moments of the non-local competition kernel is sufficient to provide a reasonable approximation of the spatial period of prey patterns that may form. In addition, the full kernel does indeed exhibit an increase in spatial period with increasing standard deviation σ, as predicted from the caricature system. The coexistence state The analysis for the linear stability of the coexistence state is not as analytically tractable as the prey-only state and for the most part we must proceed numerically, as for the system with the full kernels. For the coexistence state of the caricature system, the linear stability is determined by 85 35 30 spatial period 25 20 15 10 5 0 5.5 6 6.5 7 σ 7.5 8 8.5 9 9.5 Figure 3.4: Wavelengths of stationary prey patterns. The solid (blue) line is the predicted wavelength using linearized stability with the uniform kernel, the dotted (red) line is the prediction using the caricature system and K = 1.8. Point symbols are the spatial periods of patterns observed in numerical simulations with the uniform kernel, with initial condition h(x, 0) = 1 + 0.0001 sin(kx) where 2π/k is the predicted wavelength. 86 the two curves of eigenvalues (3.7). For the three dimensional parameter space here, we numerically computed the stability boundary in K for a grid of (B, σ) values using Maple [20]. Figure 3.5 shows how the Turing and Turing-Hopf instability boundaries for the homogeneous coexistence state of the caricature system vary with these parameters. We note that in this figure the blue (more horizontal) surface corresponds to Turing instabilities: local maxima of <(λ) at which λ(k) = 0 for some k 6= 0. The red (more vertical) surface corresponds to Turing-Hopf instabilities: local maxima of <(λ) at which λ(k) = iω for some k 6= 0, ω 6= 0. Since these correspond to only local maxima, for given parameter values the upper (i.e. higher K) of the two surfaces in the figure is the primary instability, and parameter values above both surfaces correspond to a linearly stable homogeneous state, whereas below either surface the state is unstable. Also, for both models the Hopf bifurcation in the kinetics occurs at the leftmost edge of the range shown for B. While the stability boundaries for the caricature system do depend on the local kinetics terms through the ai,j terms, we note that there are some features of these surfaces in common for both models, and so may be general for oscillatory predator-prey dynamics. First, the occurrence of Turing instability (and hence potential stationary patterns) is very strongly dependent on the kurtosis of the non-local kernel. For the caricature system of both models there is a threshold of kurtosis above which no Turing instability exists, even for very large-standard deviation kernels. This threshold is however much higher for model (i) than (ii), suggesting that pattern formation in such a system is not entirely determined by the prey non-local competition, but is truly due to a combination of processes. As for the prey-only state, it would be beneficial to study how the wavelength of stationary spatially periodic patterns in the full system compare with those predicted by the leading mode in the caricature system. Unlike the prey-only system, we have not been able to show analytically how this wavelength varies with σ. However, in Figure 3.6 we show how the leading mode wavenumber varies with the standard deviation σ of the kernel and the kinetic parameter B at the critical value of the kurtosis K. 87 M odel (i) 6 5 4 kurtosis 3 2 1 00.1 0.105 0.11 B 0.115 0.12 1 0.125 2 3 4 5 6 7 8 9 standard deviation M odel (ii) 6 5 4 kurtosis 3 2 1 0 0.024 0.026 0.028 0.03 0.032 0.034 B 0.036 0.038 0.041 2 3 4 5 6 7 8 9 standard deviation Figure 3.5: Coexistence state stability boundaries for the caricature system in B − σ − K space, where σ is the standard deviation and K is the kurtosis of the non-local competition kernel. The blue (more horizontal) surface corresponds to instability of Turing type, the red (more vertical) surface to Turing-Hopf instability. 88 M odel (i) 0.5 0.45 0.4 0.35 q 0.3 0.25 0.2 0.15 0.1 0.05 00.1 0.105 0.11 B 0.115 0.12 1 0.125 2 3 4 5 6 7 8 9 σ M odel (ii) 1 0.9 0.8 0.7 q 0.6 0.5 0.4 0.3 0.2 0.1 0 0.024 0.026 0.028 0.03 0.032 0.034 B 0.036 0.038 0.041 2 3 4 5 6 7 8 9 σ Figure 3.6: Wavenumbers at the Turing instability boundaries for the coexistence state of the caricature system in B − σ − q space. Shown here is the value of the leading mode q = k 2 from linearized stability analysis at the critical value of kurtosis K. 89 We see here that for any fixed B and for both sets of local kinetics, near the Turing instability in the kurtosis parameter, the wavenumber decreases as the standard deviation σ increases. This suggests that in the full system the wavelength of stationary spatially periodic patterns might in general be expected to increase with the standard deviation σ of the non-local kernel. We defer testing this hypothesis until the next section, where we measure the wavelength of spatially periodic patterns observed behind invasion and compare these with the predictions here. 3.4 Non-local competition and predator invasions The linear stability analysis of the homogeneous states above gives us some indication of the fate of such a system starting from the homogeneous state. However, since the evolution and selection of patterns behind predator invasions can be complex (see for instance [22, 30]), it remains to be determined if and how these instabilities play a role in determining this pattern. We study here the dependence of this spatiotemporal pattern on the non-locality of the prey intraspecific competition. To this end, we performed finite difference numerical simulations of predator invasions for a range of the parameters B and σ for the models (i) and (ii) with non-local competition kernels 1-3 on a spatial domain [−L, L] with periodic boundary conditions. The domain length L was chosen so that the predator invasion remained away from the boundary for all time and also so that for the competition kernel we have G(L) < 10−9 , since for the convolution we truncate the kernel beyond this point. The simulations used forward time and centered spatial differences with a mesh interval of ∆x = 0.25 and timestep of ∆t = 0.02. Convolutions were computed using fast Fourier transforms following the algorithms in [26, 34]. Initial conditions were taken to be as in (3.6). Each simulation was run for a total time of 5000 non-dimensional time units. For each set of parameter values we plotted the prey and predator distributions for the final 5 time units of the simulation and attempted to determine the qualitative spatiotemporal pattern emerging behind the invasion. We observed four types of spatiotemporal patterns: spatially homogeneous 90 coexistence, stationary spatially periodic patterns, wave trains and “very irregular” oscillations. The spatially homogeneous coexistence state was observed behind invasion for all parameter regions for which this state is linearly stable (to the right of all stability boundaries in Figures 3.9-3.11). Stationary periodic patterns, an example of which is shown in Figure 3.7, are a new solution form in the non-local system and were observed for model (i) with the Gaussian and uniform kernels and for model (ii) with any of 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 time → p(x,t) the three kernels. 0 100 200 300 400 x 500 600 700 800 Figure 3.7: Example of a stationary spatially periodic pattern in the wake of predator invasion. The solution for p(x, t) is plotted at times t = 2420, 2428, · · · , 2500 as a function of space x, with vertical separation proportional to the time interval between solutions. Only the right side of the spatial domain is shown. From numerical simulation of model (i) with a Gaussian kernel and B = 0.11, σ = 8.44. Wave trains were observed behind the invasion for a large parameter space for both models and with any of the three kernels. In some cases these wave trains were unstable and became irregular near the origin, or transitioned to a different wave train as in Figure 3.1(c) and (d), respectively. For the assay here, we still considered these wave train solutions as long as a band of regular wave train was observed behind the invasion front, as in any of Figures 3.1(b)-(d). Finally, for both models with any of the kernels, there were parameters for which the spatiotemporal patterns in the wake of the invasion were so irregular we were unable to detect even a transient coherent pattern. A typical example is shown in Figure 3.8. In these cases we simply called the pattern “very irregular.” 91 1.2 time → p(x,t) 1 0.8 0.6 0.4 0.2 0 800 1000 1200 1400 1600 x Figure 3.8: Example of “very irregular” oscillations in the wake of predator invasion. The solution for p(x, t) is plotted at times t = 4980, 4982, · · · , 5000 as a function of space x, with vertical separation proportional to the time interval between solutions. Only the interval [750, 1750] of the spatial domain is shown. From numerical simulation of model (i) with a Gaussian kernel and B = 0.03, σ = 6.42. 92 The results of this assay are shown in Figure 3.9 for the Laplace kernel, Figure 3.10 for the Gaussian kernel and Figure 3.11 for the uniform kernel. We see in Figures 3.9-3.11 a number of trends of interest. First, at least locally in parameter space the instabilities of the spatially homogeneous coexistence state were indicative of the pattern observed behind the invasion. Specifically, if we consider the parameter B as a bifurcation parameter, we see that a Turing-Hopf instability resulted in the emergence of wave trains behind the invasion, and a Turing instability resulted in stationary spatially periodic patterns. Second, for both models increasing the standard deviation of the non-local kernel led to an increased propensity of the system for pattern formation (i.e. a larger parameter range of B that formed spatiotemporal patterns). Similarly, decreasing the kurtosis for either model also increased the parameter range for spatiotemporal patterns. In summary, at least for these models, wider range and lower kurtosis kernels led to an increased tendency to form non-homogeneous patterns behind invasions. In addition, when stationary spatially periodic patterns were observed behind invasion (for model (i) with Gaussian and uniform kernels only) we computed the approximate wavelength of the spatial pattern at the final timepoint using a simple numerical algorithm that measures the distance between successive troughs. Briefly, we found that for parameter values nearest the Turing instability boundary these stationary patterns had wavelengths close to those predicted by the linear stability analysis of the coexistence state. We also found that in this region the caricature system with the appropriate kurtosis value (K = 3 for the Gaussian kernel and K = 1.8 for the uniform kernel) was also a very good predictor of the wavelength. An example of these predicted and observed values is shown in Figure 3.12. In addition, for all values of B that exhibited stationary spatial patterns, we found that the wavelength of spatial pattern increased with the standard deviation σ of the kernel. This is also in agreement with the general pattern observed for the caricature system. 93 Model (i) 9 8 7 σ 6 5 4 3 2 1 0 0 0.02 0.04 0.06 B 0.08 0.1 0.12 Model (ii) 9 8 7 σ 6 5 4 3 2 1 0 0 0.005 0.01 0.015 0.02 B 0.025 0.03 0.035 0.04 Figure 3.9: Stability boundaries and observed patterns following predator invasion with the Laplace competition kernel for various B and σ. Solid (grey), dashed (red) and dot-dashed (blue) curves respectively indicate Hopf, Turing-Hopf and Turing instability boundaries. Solid (grey) circles, open (red) circles, ‘x’ (blue) and ‘+’ (green) symbols respectively denote spatially homogeneous, wave train, stationary spatially periodic, and “very irregular” oscillations behind the invasion. 94 Model (i) 9 8 7 σ 6 5 4 3 2 1 0 0 0.02 0.04 0.06 B 0.08 0.1 0.12 Model (ii) 9 8 7 σ 6 5 4 3 2 1 0 0 0.005 0.01 0.015 0.02 B 0.025 0.03 0.035 0.04 Figure 3.10: Stability boundaries and observed patterns following predator invasion with the Gaussian competition kernel for various B and σ. Solid (grey), dashed (red) and dot-dashed (blue) curves respectively indicate Hopf, Turing-Hopf and Turing instability boundaries. Solid (grey) circles, open (red) circles, ‘x’ (blue) and ‘+’ (green) symbols respectively denote spatially homogeneous, wave train, stationary spatially periodic, and “very irregular” oscillations behind the invasion. 95 Model (i) 9 8 7 σ 6 5 4 3 2 1 0 0 0.02 0.04 0.06 B 0.08 0.1 0.12 Model (ii) 9 8 7 σ 6 5 4 3 2 1 0 0 0.005 0.01 0.015 0.02 B 0.025 0.03 0.035 0.04 Figure 3.11: Stability boundaries and observed patterns following predator invasion with the uniform competition kernel for various B and σ. Solid (grey), dashed (red) and dot-dashed (blue) curves respectively indicate Hopf, Turing-Hopf and Turing instability boundaries. Solid (grey) circles, open (red) circles, ‘x’ (blue) and ‘+’ (green) symbols respectively denote spatially homogeneous, wave train, stationary spatially periodic, and “very irregular” oscillations behind the invasion. 96 26 full kernel caricature observed 24 spatial period 22 20 18 16 14 12 10 4 5 6 7 σ 8 9 10 Figure 3.12: Spatial period of stationary periodic patterns behind invasion. Shown here as solid circles are the spatial periods of observed patterns from simulations of model (i) with a uniform kernel and B = 0.12. The period predicted by the linear stability of the coexistence state using the full uniform kernel is shown as the solid (blue) curve, using the caricature system with K = 1.8 is the dashed (red) curve. 97 3.5 Discussion In this study we considered the effect of non-local intraspecific prey competition on the formation of spatiotemporal patterns behind predator invasions. In particular, we focused on two main properties of the non-local competition kernel: the standard deviation and kurtosis. With respect to the standard deviation, we found here that the qualitative spatiotemporal pattern observed behind predator invasions is unchanged from the local competition case when the kernel is sufficiently narrow (low standard deviation). If, however, the interaction has a large standard deviation we observe that new spatiotemporal patterns, in particular stationary periodic states, may emerge behind the invasion. Furthermore, we found that when stationary patterns emerged, larger standard deviation kernels produced a longer wavelength pattern than smaller standard deviation ones. This suggests that in natural systems prey species that compete intraspecifically over longer ranges may be more likely to form stationary periodic patterns behind invasions, and that these patterns may have a larger spatial period than for species with shorter ranges. Our consideration of three different non-local kernels and our study of the caricature system have also helped to underscore the significance of higher order moments, in particular the kurtosis, of the non-local interaction on the generation of these patterns. Here we find that more platykurtic kernels (lower kurtosis) have a much greater destabilizing effect than leptokurtic ones. Indeed, this is in accordance with studies of platykurtic kernels for other types of models, for instance [5, 16]. Ecologically, this means that systems where the strength of competition between individuals falls off initially slowly with distance are more likely to form patterns than when this competitive effect falls off rapidly away from the focal individual, and that including such an effect in the modelling framework may be more important in the platykurtic case. With respect to predator invasion we have also found here that the instabilities of the spatially homogeneous coexistence steady state were indicative of the spatiotemporal patterns observed following predator inva98 sions, at least in the local parameter region near the onset of the instability. This is not necessarily obvious, since the mechanism of selection of patterns behind predator invasions is complex and not yet fully understood even in the local system [22, 30]. In the case of stationary spatially periodic solutions, we found that the wavelength of these spatial patterns was well-approximated by the linear stability analysis at the coexistence state in both the full and caricature systems, and so the wavelength of these patterns increased with the standard deviation of the non-local competition. For the case of wave trains behind the invasion, however, even in the local system a one-parameter family of wave trains (parameterized by amplitude or wavelength) can exist with only a single member being selected [22, 30]. To understand the influence of non-local prey competition on the properties of wave train solutions, we therefore need to determine how wave train selection occurs in these systems, and this is in fact the topic of a separate study by the authors and is currently in preparation [21]. In summary, in our study here of predator invasion into prey with nonlocal intraspecific competition we find that if the competition is sufficiently non-local in the sense that the kernel has sufficiently high standard deviation or low kurtosis then this non-locality facilitates the formation of spatiotemporal patterns behind predator invasions. This result mirrors those found in previous studies of both single species and multi-species systems with non-local interactions and adds to the mass of literature suggesting that non-local interactions have a general destabilizing effect when the interaction is in some sense sufficiently non-local. An additional contribution of the current study is our consideration of the so-called “caricature” system of Section 3.3.3 that depends only on two empirically attainable parameters, standard deviation and kurtosis, of the non-local competition kernel. In terms of steady state instabilities and spatiotemporal patterns behind invasion this reduced system was a reasonable approximation of the full system and allows for kurtosis to be treated as a continuous parameter. In the current study this allowed us to study the influence of the kurtosis more easily, but the caricature system would also be useful for attaining predictions from the model when the full form of the non-local interaction is not 99 known. Most studies so far, including this one, have focused on the bifurcation and instabilities of spatially homogeneous steady states. A clear next step is to consider the influence of non-local competition on the stability of more complex spatiotemporal patterns. Of particular interest to us is the effect on the stability of wave trains selected behind predator invasions, as these were a frequently observed solution form in this study and the instability of wave trains is suggested to be a route to spatiotemporal chaos [25, 29, 32]. We believe that the stability of wave trains could be at least numerically established by adapting the methods of Rademacher et al (2007) to nonlocal systems [27]. In fact, the authors are currently studying this problem for model (ii) with the Laplace kernel, but have not been able yet to study this problem for more general non-local competition kernels [21]. Finally, we comment on the limitations of our results here and suggest more work to be done on these types of models in general. First, we have considered only fixed parameter transects and a few forms of the non-local competition, and so although we have no indications otherwise it is not guaranteed that our results hold for all parameter regions. Second, we have restricted ourselves here to consideration of invasion when the state ahead of the predator invasion is spatially homogeneous. However, as mentioned in Section 3.3, ahead of the invasion there may be a stationary periodic pattern. In this case, constant speed predator invasion fronts may not exist, but invasion can take place nonetheless. An example of this for model (i) along with the invasion in the corresponding local model is shown in Figure 3.13. From this example, we observe here that the invasion is occurring in a “patchy” manner due to the pre-patterned state of the prey, and in comparison to the local system the invasion occurs more slowly. Also note in this figure that the predator is nearly entrained into the pre-pattern set by the prey. There is some difference between the wavelengths and amplitudes ahead of versus behind the invasion, but this difference is relatively small. A variety of approaches have been considered for modelling the spread of a species into a heterogeneous environment [15]. We think that predator-prey models with non-local intraspecific prey com100 (a) non-local kinetics (i) 0.5 p(x,t) 0.4 0.3 0.2 0.1 0 0 50 100 150 200 250 300 350 200 250 300 350 200 250 300 350 200 250 300 350 h(x,t) x 4 3.5 3 2.5 2 1.5 1 0.5 0 0 50 100 150 x (b) local kinetics (i)’ 0.5 p(x,t) 0.4 0.3 0.2 0.1 0 0 50 100 150 h(x,t) x 4 3.5 3 2.5 2 1.5 1 0.5 0 0 50 100 150 x Figure 3.13: Example of a “patchy” predator invasion. Shown are the densities of predator and prey, respectively p(x, t) and h(x, t), from numerical simulation at times t = 1750, 1800, · · · , 2000 and increasing line thickness with time. Panel (a) is for the non-local system (i) with a uniform kernel and B = 0.11, σ = 10.9. Panel (b) is for the local system (i)0 also with B = 0.11. 101 petition could also be a useful approach, since in this case the heterogeneity of the resource (prey) is intrinsic to the system. A study of this dynamical regime, including a comparison with other modelling approaches (for instance reaction-diffusion systems with spatially heterogeneous growth rates, as in [33]) would be worth pursuing. In addition, there are many generalizations of the model that could be considered. For instance, we assumed here that the predator-prey interaction was purely local, but for some systems it may be desirable to also model this in a spatially averaged manner, such as in [11]. Also, if we wish to include a temporal delay we may want to consider spatiotemporal convolutions, as in [1, 4, 13, 14]. To limit our scope, we did not consider these here, but we expect that these generalizations could elicit even greater variety and complexity of spatiotemporal patterns behind predator invasions. 3.6 Conclusion We have studied here the effect of non-local intraspecific prey competition on spatiotemporal patterns behind predator invasion in two oscillatory reaction-diffusion models. In agreement with much of the literature on nonlocal interactions, we found that sufficiently non-local competition, i.e. large standard deviation and low kurtosis, can destabilize the spatially homogeneous coexistence state and lead to new spatiotemporal patterns, such as stationary spatially periodic states, behind the invasion. Our consideration of two models and three non-local kernel forms as well as a related “caricature” system illustrates that the effect of the non-locality is dependent on the particular model considered and on the standard deviation of the kernel, but that also low kurtosis of the kernel may be necessary for pattern formation. This suggests that non-local competition may be most important to include in models when the strength of competition between individuals varies over space in such a platykurtic manner. While this study adds to the body of literature showing that non-local competition can destabilize homogeneous steady states and lead to pattern formation, we believe that there is still much to explore, including the stability and selection of spa102 tiotemporal patterns behind invasion, and the invasion of a predator into a spatially patterned prey state. Acknowledgements This work was partially supported by funding from the Natural Sciences and Engineering Research Council of Canada (NSERC). As well, S.M. would like to acknowledge the Pacific Institute for the Mathematical Sciences (PIMS) for funding through the IGTC fellowship program. We would also like to thank Michael Doebeli for his suggestion that we consider non-local interactions in these models, and Mark Lewis for his suggestion of the momentexpansion that resulted in the caricature system. 103 References [1] P. Ashwin, M. V. Bartuccelli, T. J. Bridges, and S. A. Gourley. Travelling fronts for the kpp equation with spatio-temporal delay. Z. angew. Math. Phys., 53:103–122, 2002. [2] J. Billingham. Dynamics of a strongly nonlocal reaction-diffusion population model. Nonlinearity, 17:313–346, 2004. [3] N. F. Britton. Aggregation and the competitive exclusion principle. J. Theor. Biol., 136:57–66, 1989. [4] N. F. Britton. Spatial structures and periodic travelling waves in an integro-differential reaction-diffusion population model. SIAM J. Appl. Math., 50:1663–1688, 1990. [5] M. Doebeli, H. J. Blok, O. Leimar, and U. Dieckmann. Multimodal pattern formation in phenotype distributions of sexual populations. Proc. R. Soc. Lond. B, 274:347–357, 2007. [6] E. J. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, and X. Wang. AUTO97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont). Concordia University, 1997. [7] R. A. Fisher. The wave of advance of an advantageous gene. Annals of Eugenics, 7:355–369, 1937. [8] J. Furter and M. Grinfeld. Local vs. non-local interactions in population dynamics. J. Math. Biol., 27:65–80, 1989. [9] S. A. Gourley. Travelling front solutions of a nonlocal Fisher equation. J. Math. Biol., 41:272–284, 2000. [10] S. A. Gourley and N. F. Britton. A predator-prey reaction-diffusion system with nonlocal effects. J. Math. Biol., 34:297–333, 1996. 104 [11] S. A. Gourley, N. F. Britton, M. A. J. Chaplain, and H. M. Byrne. Mechanisms for the stabilisation and destabilisation of systems of reaction-diffusion equations. J. Math. Biol., 34:857–877, 1996. [12] S. A. Gourley, M. A. J. Chaplain, and F. A. Davidson. Spatio-temporal pattern formation in a nonlocal reaction-diffusion equation. Dynamical Systems, 16:173–192, 2001. [13] S. A. Gourley and J. W.-H. So. Dynamics of a food-limited population model incorporating nonlocal delays on a finite domain. J. Math. Biol., 44:49–78, 2002. [14] S. A. Gourley, J. W.-H. So, and J. H. Wu. Nonlocality of reactiondiffusion equations induced by delay: Biological modeling and nonlinear dynamics. J. Math. Sci., 124:5119–5153, 2004. [15] A. Hastings, K. Cuddington, K. F. Davies, C. J. Dugaw, S. Elmendorf, A. Freestone, S. Harrison, M. Holland, J. Lambrinos, U. Malvadkar, B. A. Melbourne, K. Moore, C. Taylor, and D. Thomson. The spatial spread of invasions: new developments in theory and evidence. Ecol. Lett., 8:91–101, 2005. [16] E. Hernandez-Garcia, C. Lopez, S. Pigolotti, and K. H. Andersen. Species competition: coexistence, exclusion and clustering. Phil. Trans. R. Soc. A, 367:3183–3195, 2009. [17] E. E. Holmes, M. A. Lewis, J. E. Banks, and R. R. Veit. Partial differential equations in ecology: Spatial interactions and population dynamics. Ecology, 75:17–29, 1994. [18] R. E. Kirk. Statistics: An Introduction. Holt, Rinehart and Winston, Inc., third edition, 1990. [19] M. A. Lewis. Spatial coupling of plant and herbivore dynamics: The contribution of herbivore dispersal to transient and persistent “waves” of damage. Theor. Popul. Biol., 45:277–312, 1994. 105 [20] Maplesoft. Maple 9: Learning Guide, 2003. [21] S. M. Merchant and W. Nagata. Selection and stability of wave trains behind predator invasion in a model with non-local prey competition. in preparation, 2009. [22] S. M. Merchant and W. Nagata. Wave train selection behind invasion fronts in reaction-diffusion predator-prey models. submitted, 2009. [23] J. D. Murray. Mathematical Biology I: An Introduction. SpringerVerlag, third edition, 2002. [24] S. Nadarajah. A generalized Normal distribution. J. Appl. Stat., 32:685–694, 2005. [25] S. V. Petrovskii and H. Malchow. Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics. Theor. Popul. Biol., 59:157–174, 2001. [26] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C: the Art of Scientific Computing. Cambridge University Press, second edition, 1992. [27] J. D. M. Rademacher, B. Sandstede, and A. Scheel. Computing absolute and essential spectra using continuation. Physica D, 229:166–183, 2007. [28] A. Sasaki. Clumped distribution by neighborhood competition. J. Theor. Biol., 186:415–430, 1997. [29] J. A. Sherratt. Unstable wavetrains and chaotic wakes in reactiondiffusion systems of λ − ω type. Physica D, 82:165–179, 1995. [30] J. A. Sherratt. Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave. Physica D, 117:145–166, 1998. [31] J. A. Sherratt, B. T. Eagan, and M. A. Lewis. Oscillations and chaos behind predator-prey invasion: mathematical artifact or ecological reality? Phil. Trans. Roy. Soc. Lond. B, 352:21–38, 1997. 106 [32] J. A. Sherratt, M. A. Lewis, and A. C. Fowler. Ecological chaos in the wake of invasion. Proc. Natl. Acad. Sci., 92:2524–2528, 1995. [33] N. Shigesada, K. Kawasaki, and E. Teramoto. Traveling periodic waves in heterogeneous environments. Theor. Popul. Biol., 30:143–160, 1986. [34] W. T. Vetterling, S. A. Teukolsky, W. H. Press, and B. P. Flannery. Numerical Recipes Example Book (C). Cambridge University Press, second edition, 1992. 107 Chapter 4 Selection and Stability of Wave Trains Behind Predator Invasions in a Model with Non-local Prey Competition 3 4.1 Introduction A number of studies have recently considered oscillatory reaction-diffusion systems as models for predator invasions in predator-prey systems. In one spatial dimension these systems have the general form ht = Dh hxx + f (h, p), pt = Dp pxx + g(h, p), (4.1) where h(x, t) is the prey density and p(x, t) is the predator density at position x and time t, Dh and Dp are the prey and predator diffusivities and the functions f (h, p) and g(h, p), referred to as the “kinetics,” describe the local interaction of prey and predator. To model invasion, we assume that the system without diffusion has a prey-only steady state (h, p) ≡ (1, 0) and a coexistence steady state (h, p) ≡ (h∗ , p∗ ), where the predator and prey 3 A version of this chapter will be submitted. S. M. Merchant and W. Nagata, Selection and stability of wave trains behind predator invasions in a model with non-local prey competition. 108 coexist at some non-zero levels. The system is considered “oscillatory” if in the diffusionless system the coexistence state undergoes a supercritical Hopf bifurcation to stable limit cycles as some parameter of the kinetics is varied. A common invasion scenario considered is the prey initially at carrying capacity everywhere on the spatial domain and the predator introduced to a localized region, for instance h(x, 0) = 1, , p(x, 0) = 0, |x| ≤ δ (4.2) |x| > δ where δ determines the spatial extent of an initial predator introduction at density . For the oscillatory kinetics studied here, numerical simulations of system (4.1) with initial conditions (4.2) evolve to a front of predator density that travels with constant speed away from the point of origin. Examples of such fronts in numerical simulations for this study are shown in Figure 4.1. The speed cinv at which fronts invade an unstable steady state has been the subject of much study. A comprehensive review of this topic is provided by van Saarloos (2003) [33]. More recent work by many authors has focused on the dynamics behind this invasion front, as it has been shown to exhibit a range of spatiotemporal patterns, including spatially homogeneous coexistence state solutions, irregular oscillations, and wave trains [28, 29, 32]. Examples of these spatiotemporal patterns for the model considered in this study are shown in Figure 4.1. Wave trains, also known as periodic travelling waves, are spatially periodic solutions that move through space at a constant speed while maintaining their shape. These solutions have a spatial period or wavelength L, a temporal period T and a phase speed given by c̄ = L/T . The existence of wave trains for an oscillatory reaction-diffusion system was first studied by Kopell and Howard (1973) for λ − ω systems and has since been established for generic oscillatory reaction-diffusion systems [13, 20]. In the past decade, a number of field studies have shown evidence of periodic travelling waves in natural populations, as summarized in [30]. In the same time period, theoretical work on wave train solutions behind invasion has also significantly 109 (a) B = 0.0280, σ = 2.0 time → p(x,t) 0.15 0.1 0.05 0 0 500 1000 1500 2000 1500 2000 1500 2000 1500 2000 x (b) B = 0.0110, σ = 3.0 time → p(x,t) 0.15 0.1 0.05 0 0 500 1000 x (c) B = 0.0045, σ = 4.0 time → p(x,t) 0.15 0.1 0.05 0 0 500 1000 x (d) B = 0.0180, σ = 7.5 time → p(x,t) 0.15 0.1 0.05 0 0 500 1000 x Figure 4.1: Spatiotemporal patterns behind invasion fronts. The solution for p(x, t) is plotted as a function of space x at ten fixed times, with vertical separation proportional to the time interval between solutions. From simulation of (4.4) with kernel (4.5) and parameter values in the text and as indicated. In (a) the invasion front decays to the spatially homogeneous coexistence state and remains there. (b)-(d) show the generation of wave trains in the invasion wake. In (b) the front first decays to the coexistence state followed by a subsequent transition to the wave train, whereas in (c) the front transitions directly to a wave train. In (d) the selected wave train is unstable and breaks up into irregular oscillations. 110 advanced, with two main areas of focus: (1) selection and (2) stability of such wave trains [30]. Both elements are necessary for an understanding the dynamics behind predator invasion, and we will study both topics for the system considered here. First, however, we provide some background on selection and stability for oscillatory reaction-diffusion systems. For the system (4.1) wave trains may exist as a one-parameter family beyond the Hopf bifurcation in the kinetics. To see this, transformation of the system to a moving frame of speed c̄ and subsequent reduction of the resulting ODE system to first order via y = h0 (ξ) and z = p0 (ξ) yields the four-dimensional ODE system dh dξ dp dξ dy dξ dz dξ = y, = z, = = 1 Dh 1 Dp (−c̄y − f (h, p)) , (4.3) (−c̄z − g(h, p)) . Wave train solutions of speed c̄ correspond to limit cycle solutions of system (4.3). We therefore have that for given fixed kinetic parameters wave trains may exist for a range of possible speeds c̄ and corresponding wavelengths L and temporal periods T . In numerical simulations, however, only a single member of this family evolves behind a predator invasion, as shown in Figures 4.1(b) and 4.1(c). Sherratt (1998) first studied the mechanism by which a single member of this family is generated behind invasion, in the case where the primary invasion front first decays to the spatially homogeneous coexistence steady state, as in 4.1(b) [26]. In chapter 2, we also studied this problem and extended the prediction of Sherratt (1998) to a more general selection criterion that provides a higher order approximation and can also be applied to the case in Figure 4.1(c) where the coexistence state is not approached. Wave trains selected behind predator invasions are not always stable, as shown in Figure 4.1(d). In the case where an unstable wave train is selected it can only be observed transiently, and numerical simulations suggest that the long term solution in such a case can have irregular spatiotemporal oscil- 111 lations. A number of studies have suggested that these irregular oscillations are in fact spatiotemporal chaos [17, 25, 29]. There are few analytical results on the stability of wave train solutions to oscillatory reaction-diffusion systems. For λ − ω systems, Kopell and Howard (1973) have shown that wave trains above a threshold in amplitude are stable and those below are unstable [13]. Since near the Hopf bifurcation in the kinetics any oscillatory reaction-diffusion system with equal diffusion coefficients can be approximated by such a λ − ω system, we also have that high amplitude wave trains are stable in this kinetic parameter regime. In fact, this approximation was used by Sherratt (2001) to predict whether predator invasions would generate stable periodic waves or spatiotemporal chaos for a number of natural populations [27]. In the more general case when the kinetics are not near the Hopf bifurcation or the diffusion coefficients are not sufficiently similar, however, there have been no such analytical results that divide the wave train family into stable and unstable bands. Methods have recently been developed for the numerical study of wave train stability using numerical continuation software such as auto [3, 5, 19, 22]. Smith and Sherratt (2007) applied these methods to both the λ − ω system and to an oscillatory predator-prey system [32]. The authors used numerical continuation to separate wave train solutions of the full predatorprey system into stable and unstable bands. Far from the kinetic bifurcation, the stability predictions for the full system deviated substantially from those predicted from the λ − ω system, and this knowledge of wave train stability was valuable in understanding observed behaviour of simulations, particularly for unstable selected wave trains. In this chapter, we study the selection and stability of wave trains behind invasion in a model derived from oscillatory reaction-diffusion predator-prey systems. In a previous study, we considered the influence of non-local prey competition on the spatiotemporal patterns formed following predator invasions [15]. In these systems intraspecific competition for the prey takes place over the entire spatial domain, rather than only at the same point. In this case the effective prey density hef f (x, t) for the competition term at a 112 given location is modelled as a spatial average Z ∞ G(x − y)h(y, t)dy, hef f (x, t) := (G ∗ h)(x, t) = −∞ where the non-local competition kernel G(z) here describes how the strength of competition varies with the distance z between individuals. Non-local competition has been considered in a number of other models for both single species and interacting populations [2, 4, 6–8, 24]. A general consensus from these studies is that non-local interactions can destabilize spatially homogeneous steady states and generate new spatiotemporal patterns, including stationary spatially periodic, standing wave, and wave train solutions. Our study in chapter 3, [15], mirrored these results, finding that sufficiently nonlocal prey competition, characterized by the standard deviation and kurtosis of the competition kernel, could generate stationary spatially periodic states behind predator invasions as well as extend the parameter regions for which wave train solutions were observed behind invasion. The study [15] was purely qualitative and did not consider how the various wave train parameters and stability were influenced by the non-locality. We therefore study this problem further here for one of the predator-prey systems considered in [15]. In Section 4.3 we numerically study wave train stability using some modification of the methods in [19]. We illustrate that for the family of wave trains, the region of stable waves shrinks dramatically as the degree of non-locality (standard deviation of the non-local competition kernel) increases, suggesting that irregular oscillations may be more prominent in such systems. To predict the dynamics for given parameters, we also need a prediction of the wave train selected behind predator invasion. In Section 4.4, we reconsider and generalize the selection criterion from [16] to also provide an approximation for such non-local systems. This revised selection criterion is then tested by comparison with observations from numerical simulations. While it does provide more accurate predictions than the local criterion, we examine why it clearly fails in some cases. Finally, in Section 4.5 we bring together the selection and stability results to provide a prediction of the dynamics that may be expected behind invasion 113 and suggest future work in this area. 4.2 The model and wave train families We study here one of the non-local reaction-diffusion predator-prey systems considered in [15] given by hp ht = Dh hxx + h(1 − G ∗ h) − A h+C , p pt = Dp pxx + Bp 1 − h , (4.4) where A, B, C > 0 are biological parameters. For the non-local competition kernel, we take the Laplace, or double-exponential probability density function with standard deviation σ given by [35] √ G(x) = 2 − e 2σ √ 2|x| σ . (4.5) We will sometimes refer to the corresponding “local system”, obtained by setting G(x) = δ(x), and to the “local kinetics” f (h, p) and g(h, p) given by p g(h, p) = Bp 1 − . h hp , f (h, p) = h(1 − h) − A h+C There are a number of reasons to consider kernel (4.5). For one, Gourley (2001) has shown how this form of kernel can arise naturally by considering resource-consumer dynamics [8]. Another reason is one of mathematical 2 convenience. If we define w(x, t) = (G∗h)(x, t), we then have wxx = 2 (w − σ h) and periodic solutions of (4.4) are in direct correspondence with periodic solutions of the three-dimensional system hp ht = Dh hxx + h(1 − w) − A h+C , p pt = Dp pxx + Bp 1 − h , 0 = wxx + 2 (h σ2 (4.6) − w). We may then compute wave train solutions of a given speed c̄ for the original system by converting (4.6) to a six-dimensional first order ODE system by 114 making the coordinate transformation ξ = x − c̄t and taking a = h0 (ξ), b = p0 (ξ) and d = w0 (ξ). The resulting system is thus h0 = a, p0 = b, w0 = d, hp a0 = D1h −c̄a − h(1 − w) + A h+C , p b0 = D1p −c̄b − Bp 1 − h , d0 = 2 (w σ2 (4.7) − h). In our study here we compute wave train solutions of fixed speed c̄ numerically by continuation in the kinetic parameter B from the Hopf bifurcation point of the coexistence steady state (h, p, w, a, b, d) ≡ (h∗ , p∗ , h∗ , 0, 0, 0) of system (4.7). For this, we use the numerical continuation software auto [5]. The system (4.4) with competition kernel (4.5) is one of the systems studied in [15]. For convenience, we take the same parameter choices as in [15]: we fix Dh = Dp = 1, A = 6.92 and C = 0.64 and vary the parameter B and the standard deviation σ of the non-local kernel. For these parameters, Figure 4.2 shows the primary instability of the spatially homogeneous coexistence state as either the kinetic parameter B or the standard deviation σ of the kernel are varied. We note that the nature of the primary instability changes at a critical value σ ∗ of the standard deviation, with σ ∗ ≈ 6.47. For small standard deviations, σ < σ ∗ , the primary instability is a Hopf instability at the kinetic parameter value B ≈ 0.02268, as for the local system. For σ > σ ∗ , however, the instability is instead of Turing-Hopf type, with onset at larger values of the kinetic parameter B. Before we proceed to the selection and stability of wave trains, we first comment here on the influence the non-local competition term has on the family of wave train solutions. In Figure 4.3 we show the wavelength L at the Hopf bifurcation point for system (4.7) in the kinetic parameter B for a range of standard deviations σ. For a given wavelength L, the curves in Figure 4.3 indicate the parameter value B at which wave trains of wavelength L bifurcate from the spatially 115 10 8 unstable σ 6 stable 4 2 0 0.01 0.015 0.02 0.025 B 0.03 0.035 0.04 Figure 4.2: Stability boundaries for the spatially homogeneous coexistence state. The solid curve corresponds to onset via complex eigenvalues at zero wavenumber (Hopf instability). The dashed curve corresponds to onset via complex eigenvalues at non-zero wavenumber (Turing-Hopf instability) 116 350 local σ=2 σ=4 σ=6 σ=8 σ=10 300 250 L 200 150 100 50 0 0 0.005 0.01 0.015 0.02 B 0.025 0.03 0.035 0.04 Figure 4.3: Wavelengths L at Hopf bifurcation points for various σ. Point symbols at the left edge of curves denote the edge of existence of the Hopf bifurcation curve. 117 homogeneous solution as B is decreased. Wave trains therefore exist to the left of these Hopf bifurcation points. We note that for sufficiently small L in each case the Hopf bifurcation curve ceases to exist because the speed c̄ at bifurcation approaches zero (indicated by filled circles in Figure 4.3). We also see from Figure 4.3 that increasing the standard deviation σ of the nonlocal competition kernel leads to Hopf bifurcation at shorter wavelengths L in addition to new parameter regions for B at which wave trains exist. We find here that the rightmost point in this case corresponds to onset of Turing-Hopf instability shown in Figure 4.2, suggesting that this instability leads also to the bifurcation of wave trains. 4.3 Wave train stability Here we consider the spectrum of the operator L obtained by linearizing the integro-differential system (4.4) with Dh = Dp = 1 in the comoving frame of speed c about a stationary spatially periodic solution (h0 (x), p0 (x)). We begin by taking h = h0 + u, p = p0 + v. Expanding in Taylor series where necessary and retaining only first order terms, we obtain the linear system 0 p0 ut = uxx + cux − h0 G ∗ u + 1 − w0 + (hAh 2 − 0 +C) 2 0 vt = vxx + cvx − B hp00 u + B 1 − 2p h0 v, Ap0 h0 +C u− Ah0 h0 +C v, where w0 = G ∗ h0 . We therefore have the linear operator L defined by L ! u v = ! uxx + cux − h0 G ∗ u vxx + cvx +A ! u v , where A(x) = a1,1 a1,2 a2,1 a2,2 ! 0 p0 1 − w0 + (hAh 2 − 0 +C) 2 = −B hp00 Ap0 h0 +C 0 − hAh 0 +C . 0 B 1 − 2p h0 118 4.3.1 Computation of essential spectra for wave trains In their study of wave train solutions of reaction-diffusion predator-prey systems, Smith and Sherratt (2007) numerically computed the continuous part, or essential spectrum, of the operator linearized about the wave train to determine wave train stability [32]. Methods for the numerical computation of essential spectra for reaction-diffusion operators linearized about a spatially periodic steady state are described by Rademacher et al. (2007) [19]. While the system (4.4) is not exactly of the form studied there, we can apply here their methods with some small modifications, and refer the reader to [19] for details of the method and its background. To derive the eigenvalue problem, we assume perturbations of the form u = eλt+νx ū(x), v = eλt+νx v̄(x) for ν = iγ, γ ∈ R. The eigenvalue problem is thus (upon dropping the bars) λu = (∂x + ν)2 u + c(∂x + ν)u + (a1,1 − w0 )u + a1,2 v − h0 q, λv = (∂x + ν)2 v + c(∂x + ν)v + a2,1 u + a2,2 v, where here we have set q = e−νx G ∗ (eνx u). As for the computation of wave trains, q satisfies the relation 0 = (∂x + ν)2 q − 2 (q − u), σ2 (4.8) and for periodic u(x) this uniquely determines q to be the convolution above. We may therefore compute essential spectra for wave trains using the methods of [19] with (4.8) as an additional equation in the continuation. Finally, we cast this as a first-order problem and have the following system for the 119 computation of essential spectra for wave train solutions of (4.4): u0 = e − νu, v 0 = f − νv, q 0 = g − νq, e0 = −νe − ce + (w − a1,1 )u + λu − a1,2 v + hq, (4.9) f 0 = −νf − cf − a2,1 u − a2,2 v + λv, g 0 = −νg + 2 (q σ2 − u), where (h(x), p(x), w(x)) is the spatially periodic stationary solution in the comoving frame of speed c. We note that in practice, we must therefore simultaneously solve the system (4.7) for the wave train (h(x), p(x), w(x)). In addition, for a wave train with wavelength L the system (4.9) is solved as a boundary value problem in real variables on a periodic domain of length L with periodic boundary conditions for all components. Computation of initial solutions and continuation of system (4.9) for the essential spectra of wave trains are otherwise the same as described in [19], and so we do not reproduce this here. We computed essential spectra for various sets of parameters (B, σ) and wave train wavelengths L. Both stable and unstable spectra were observed, and examples of the rightmost curve of eigenvalues are shown in Figure 4.4(a). We note that due to translation invariance this curve always passes through the origin in the complex plane. We found that, as in the study [32], in all cases examined the transition from stable to unstable wave trains occurs through an instability to a band of wavenumbers γ near zero, as shown in Figure 4.4(b). For the local system and near the kinetic Hopf bifurcation, this instability is known to correspond to an Eckhaus instability of the associated complex Ginzburg-Landau amplitude equation [1, 32]. Eckhaus instabilities have been extensively studied and are known to have various effects, including amplitude modulation, transitions to new wavelength solutions, or to spatiotemporal chaos [9, 10, 14]. While we do not study the nature of this instability further here, we do note that Rademacher et al. (2007) also outline methods for continuing the onset of such instabilities. 120 We study this stability boundary further in the next section and note that, as in [32], for convenience we continue to refer to sideband instability of the form in Figure 4.4 as “Eckhaus” instability even well away from the Hopf bifurcation. (a) (b) 0.3 2e-04 (i) (ii) (i) (ii) 0.2 1e-04 Re(λ) Im(λ) 0.1 0 0 -0.1 -1e-04 -0.2 -0.3 -2e-04 -1e-04 0 1e-04 -0.04 -0.02 0 0.02 0.04 γ Re(λ) Figure 4.4: Essential spectra for wave train solutions of (4.4). Panel (a) shows the essential spectra in the complex plane. Panel (b) shows the Re(λ) against the wavenumber γ. Solid (blue) curves are for σ = 3.0, B = 0.01268 and L = 203.224 and represent a spectrally stable wave train. Dashed (red) curves are for σ = 3.0, B = 0.02118 and L = 360.694 and represent an unstable wave train. 4.3.2 Computation of the Eckhaus instability boundary Following the same process as [19], we develop here the system to be used for continuation of the curve in B − L space that describes onset of wave train instability of Eckhaus type. We define the two quantities λ| , λ|| ∈ R by λ| := dλ0 dν , ν=0 λ|| := d2 λ0 dν 2 . ν=0 The sign of λ|| determines whether the rightmost curve λ0 (ν) of eigenvalues extends into the positive half of the complex plane, and hence allows us to detect Eckhaus instabilities, with onset when λ|| = 0. To compute the quantity λ| we differentiate system (4.9) with respect to ν = iγ and evaluate 121 the resulting equations at γ = 0. We thereby obtain the system u0| = −u + e| , v|0 = −v + f| , q|0 = −q + g| , e0| = −e − ce| + (w − a1,1 )u| − a1,2 v| − hq| + λ| u, (4.10) f|0 = −f − cf| − a2,1 u| − a2,2 v| + λ| v, g|0 = −g + 2 (q σ2 | − u| ), for computation of λ| , where here u| := ∂ν u, v| := ∂ν v, · · · ∈ R. For λ|| , we differentiate (4.9) with respect to ν twice and evaluate at γ = 0. u0|| = −2u| + e|| , v||0 = −2v| + f|| , q||0 = −2q| + g|| , e0|| = −2e| − ce|| + (w − a1,1 )u|| − a1,2 v|| − hq|| + λ|| u + 2λ| u| , (4.11) f||0 = −2f| − cf|| − a2,1 u|| − a2,2 v|| + λ|| v + 2λ| v| , g||0 = −2g| + 2 (q σ 2 || − u|| ), where u|| := ∂ν2 u, v|| := ∂ν2 v, · · · ∈ R. We note that for computation of the boundary of Eckhaus instability systems (4.10) and (4.11) are also cast as boundary value problems on a periodic domain of length L. We then continue in the parameters (B, L) solutions of systems (4.9), (4.10) and (4.11) at which λ = 0 and λ|| = 0 together with system (4.7) for the wave train solution (h, p, w). As for essential spectra, computation of initial solutions and continuation with auto is performed following the methods of [19]. For each of the values of σ for the non-local kernel we computed this stability boundary in the (B, L) plane to study how the stability of the wave train family varies with the degree of non-locality of the prey nonlocal competition. In Figure 4.5 we show this boundary for representative values of σ. From Figure 4.5 we see that for small values of σ and for all B relatively short wavelength wave trains are unstable, whereas sufficiently long wavelength wave trains are linearly (spectrally) stable. As the stan122 350 350 (a) σ = 2.0 300 (d) σ = 7.0 300 250 200 200 150 150 100 100 50 50 L 250 0 0 0.005 0.01 0.015 0.02 0.025 0.03 350 0.005 0.01 0.015 0.02 0.025 0.03 350 (b) σ = 4.0 300 (e) σ = 7.15 300 250 200 200 150 150 100 100 50 50 L 250 0 0 0.005 0.01 0.015 0.02 0.025 0.03 350 0.005 0.01 0.015 0.02 0.025 0.03 350 (c) σ = 6.0 300 (f) σ = 7.25 300 250 200 200 150 150 100 100 50 50 L 250 0 0 0.005 0.01 0.015 0.02 B 0.025 0.03 0.005 0.01 0.015 0.02 0.025 0.03 B Figure 4.5: Eckhaus instability boundary for wave trains in B − L space. Horizontal axis is the kinetic parameter B. Vertical axis is the wavelength L. Dashed lines indicate the Hopf bifurcation in the corresponding system (4.7). Shaded areas indicate linearly stable wave trains. 123 dard deviation σ of the non-local competition kernel is increased, the lower boundary, i.e. minimal wavelength L, for stable wave trains decreases, but a region of long wavelength wavetrains becomes unstable for small values of B. As the standard deviation is increased still further this region grows until eventually for all B there are unstable wave trains for large values of the wavelength L. At this point there is a closed region of stable wave trains for small wavelengths L that continues to shrink as the standard deviation σ is increased still further. 4.4 Wave train selection In [16] we conjectured a selection criterion for wave trains generated behind predator invasions in oscillatory reaction-diffusion systems. The essential idea behind this criterion is that in the comoving frame of a particular speed ccoh (assumed positive for convenience) the selected wave train has the frequency ωcoh of the linear unstable oscillatory mode of the coexistence state, termed the “pacemaker” frequency. Specifically, if we let L > 0 be the spatial period, cwtrain be the phase speed and T = L/cwtrain be the (signed) temporal period of the selected wave train then the frequency of the wave train in the comoving frame of speed ccoh is (2π/L)ccoh − (2π/T ). The criterion for the selected wave train is therefore that it satisfies 2π 2π ccoh − = ωcoh , L T (4.12) where the frequency ωcoh is the imaginary part of the eigenvalue of the linearization of the kinetics system about the coexistence state. The speed ccoh is the speed of the “relevant” front. In Figure 4.1(b) the invasion appears to involve two travelling fronts, a primary front invading the unstable prey-only state, and a secondary front invading the unstable coexistence state. These two fronts do not generally travel at the same speed. In this case, which we refer to as case I of the selection criterion the “relevant” front is the secondary one, from the coexistence state to the selected wave train. Figure 4.1(c) illustrates a different dynamic, referred 124 to as case II. In this case the primary invasion front does not decay to the coexistence state and the “relevant” front connects directly from the preyonly state to the selected wave train. The speed of fronts propagating into unstable steady states are grouped into two classes: pulled fronts that travel at the so-called “linear spreading speed” ν ∗ (defined below) and are in some sense generic, and pushed fronts that travel at a speed ν > ν ∗ . In [33], van Saarloos (2003) defines the linear spreading speed ν ∗ ∈ R, which is given by solving the saddle point equations ∗ ∗ 0 = S(k , ω ) 0 = (∂k + ν ∗ ∂ω )S|(k∗ ,ω∗ ) 0 = =(ω ∗ − ν ∗ k ∗ ), (4.13) for (k ∗ , ω ∗ , ν ∗ ) where S(k, ω) is the characteristic equation for the linearization about the unstable steady state ahead of the front, obtained by assuming the exponential ansatz h, p ∼ ei(kx−ωt) , k, ω ∈ C for perturbations of the linearized system. While there may be multiple solutions to (4.13), only those for which <(D) > 0, where D= −i(∂k + ν ∗ ∂ω )2 S 2∂ω S (k∗ ,ω ∗ ) are relevant and when there are several dynamically relevant saddle points we take the one with the largest corresponding ν ∗ for the linear spreading speed. For pulled fronts, it is argued that for initial conditions that decay sufficiently rapidly in space, a front invading the unstable state forms with asymptotic speed ν ∗ [33]. In the selection criterion (4.12) we assume that fronts emanating from both the prey-only state and those invading the coexistence state are of the pulled variety and so travel at speeds cinv and csec given by the linear spreading speeds obtained by linearization about the prey-only and coexistence states, respectively. We argue that case I of the selection mechansim occurs if and only if csec < cinv , since otherwise the secondary front would eventually catch up to the primary invasion front. We therefore take for the 125 comoving frame speed in the criterion ccoh := min{csec , cinv }. 4.4.1 Selection in the non-local system Here we modify the selection criterion (4.12) for the non-local system (4.4). For this, we require the two quantities ωcoh and ccoh for the non-local system. The linear spreading speed ccoh is as before obtained by linearization of the system about the “relevant” spatially homogeneous steady state (h0 , p0 ). The linearization of the system (4.4) about a spatially homogeneous steady state (h0 , p0 ) is ht = Dh hxx + a1,1 h + h0 h − h0 G ∗ h + a1,2 p, pt = Dp pxx + a2,1 h + a2,2 p, where a= a1,1 a1,2 a2,1 a2,2 ! = ∂f ∂h ∂g ∂h ∂f ∂p ∂g ∂p ! (h0 ,p0 ) is the linearization matrix for the local system about the steady state. Substitution of the ansatz h, p ∼ ei(kx−ωt) then gives the characteristic equation iω − k 2 Dh + a1,1 + h0 − h0 H(k) S(k, ω) = −a1,2 a2,1 , Z iω − k 2 Dp + a2,2 (4.14) ∞ G(z)e−ikz dz. We note that since k here is complex for √ the kernel G considered here this integral exists only for |=(k)| < 2/σ and 1 in this case we have H(k) = 1 2 2. 1 + 2σ k Unlike the local system, we are unable to solve the saddle point equawhere H(k) = −∞ tions (4.13) analytically for this characteristic equation. However, they may be solved numerically, and here we solve the system (4.13) numerically by continuation with auto. We note that, at least for the parameter range considered here, linear spreading speeds for fronts emanating from the preyonly state were found to be identical to those for the local system. For fronts 126 emanating from the coexistence state, however, the linear spreading speed varied with σ, in fact decreasing with increasing σ for all cases considered here. We verified this for a number of parameter sets by comparison of the linear spreading speed obtained by continuation with speeds observed in simulations, and found good agreement. This suggests that even for the non-local system fronts continue to be of the pulled variety. The other main component of the criterion (4.12) is the “pacemaker” frequency ωcoh . For the non-local system, the diffusionless system is still spatially dependent and we find that continuing to use the frequency from the local kinetics provides inaccurate predictions of the selected wave train. In [16], the linear unstable oscillatory mode of the coexistence state, ωloc , was chosen for this frequency because it extended predictions arising from more rigourous work on the complex Ginzburg-Landau equation, as well as the selection mechanism described in [26], both valid near the Hopf bifurcation in the kinetics. It was in fact somewhat surprising that this frequency continued to provide accurate predictions well away from the Hopf bifurcation, particularly in case II where the coexistence state is not even approached. We termed this frequency the “pacemaker” frequency for this reason, it appeared to serve as a global mode in the selection problem. Further study of the linear spreading speed analysis in the local system, however, illustrates why this mode plays a fundamental role. As shown in [16], for the local kinetics and with equal diffusion coefficients the linear spreading speed equations (4.13) can be explicitly solved. The solution for fronts emanating from the coexistence steady state is always of the form ω ∗ = ωloc + iβ, k ∗ = iγ, ν∗ = β γ (4.15) for some β, γ ∈ R. In the comoving frame of speed ν ∗ with coordinates ξ = x − ν ∗ t, the front emanating from the coexistence state takes the form ∼ eikξ−i(ω−ν ∗ k)t [33]. For the linear spreading speed saddle point solution (4.15) this reduces to e−γξ−iωloc t . The so-called “pacemaker” frequency is thus the frequency with which the solution leaves the steady state in the comoving frame of the linear spreading speed, and so the criterion (4.12) 127 stipulates that the wave train in the comoving frame matches this frequency. We note that while the frequency is based on linear analysis, it does not require that the system is near the kinetic Hopf bifurcation. It does, however, require that the coexistence steady state is approached behind the primary invasion front. For the non-local system, the linear spreading point k ∗ is not purely imaginary and so the front emanating from the coexistence state in the comoving frame of speed ν ∗ instead has the temporal oscillation given by ei(α−ν ∗ η)t , where α = <(ω ∗ ), η = <(k ∗ ). For the selection criterion in the non-local system we therefore take (4.12) with the frequency ωcoh = α − ν ∗ η, (4.16) where ν ∗ , ω ∗ = α + iβ, and k ∗ = η + iγ, α, β, η, γ ∈ R are solutions of the linear spreading speed system (4.13), with linearization about the spatially homogeneous coexistence state. We note that this frequency may provide a more accurate approximation to the selected wave train also in the case of unequal diffusion coefficients, as studied in [16] since this case also results in η 6= 0. We have recomputed the selection criterion for the models and cases considered in [16] and did indeed find better agreement with the observations from numerical simulations. We now apply this revised criterion to the non-local system (4.4) and study its performance for a range of standard deviations σ for the non-local kernel. 4.4.2 Methods We consider the performance of the selection criterion for a range of standard deviations σ =1, 2, 3, 4, 5 and 6, as well as for the local system. Standard deviations σ > 6 were not considered because in [15] we found that in this case the pattern behind invasion was very irregular for a large parameter range B, and we were unable to distinguish even transient wave trains. For the parameter B we restrict to the range B ∈ [0.002, 0.02268], since to the left of this interval for some σ the Hopf bifurcation points in Figure 4.3 do not exist and wave trains cannot be continued in this case. 128 To study the validity of the revised selection criterion, for each of the standard deviation values σ =1, 2, 3, 4, 5 and 6 we computed the predicted selected wave train for cases I and II of the criterion by continuation in the kinetic parameter B of system (4.7) together with the system (4.13) for the linear spreading speed ccoh and subject to the constraint (4.12). The system continued is posed as a boundary value problem on a domain of length L and with periodic boundary conditions. We note that for the algebraic system (4.13) and the criterion (4.12) we must solve these in tandem with the spatially varying wave trains, and so these are cast as boundary value problems also with periodic boundary conditions and initialized at a spatially homogeneous solution and verified that throughout the continuation these remained spatially homogeneous. To test these predictions, we also performed numerical simulations of predator invasions and computed the wavelength of wave trains observed in the wake of the invasion. For each of the values of σ above, we performed finite difference numerical simulations of predator invasions for B = 0.0020, 0.0025, · · · , 0.0225 for the system (4.4) on a spatial domain [−l, l] with periodic boundary conditions. The domain length l was chosen so that the predator invasion remained away from the boundary for all time and so for the competition kernel we have G(l) < 10−9 , since for the convolution we truncate the kernel beyond this point. The simulations used forward time and centered spatial differences with a mesh interval of ∆x = 0.125 and timestep of ∆t = 0.005. Convolutions were computed using fast Fourier transforms following the algorithms in [18, 34]. Initial conditions were taken to be as in (4.2) with = p∗ and δ = 50. Simulations were run for a total time of 6000 non-dimensional time units for B < 0.017. For larger values of B the evolution of the selected wave train was very slow since near the kinetic Hopf bifurcation at B = 0.02268 the speed csec of the secondary invasion front approaches zero. We therefore ran simulations for 12000 time units for B ≥ 0.017. For each simulation we recorded the distribution of predator and prey densities at times t = 1000, 2000, · · · . In nearly all cases, the final time profile exhibited a wave train evolving behind the invasion front. These profiles were then pro129 cessed by a simple routine that compared locations and magnitudes of local maxima to determine the wavelength of wave trains behind the predator invasion front at the final time. In the case that a wave train could not be distinguished in the final time series (B > 0.021 with any σ and B = 0.002, 0.0025 with σ = 6) we omit these points from our results. 4.4.3 Results The results for this study are summarized in Figure 4.6 where we show the predicted wavelengths L of wave trains selected behind invasion as well as observed wavelengths from final times of simulations for competition kernels with standard deviation σ = 2, 3, 4, 5, 6 and for the local system (panel (a)). We omit the results for σ = 1, but note that they are quite similar to the local case. Also in Figure 4.6 are the Eckhaus stability boundaries for each case, computed using the methods of Section 4.3. From this figure we observe several trends. First, for case I of the selection criterion there is close agreement with observations, in particular for regions near the kinetic Hopf bifurcation. Farther from this point, for small B, the observations deviate from the prediction with increasing deviation for larger σ. Case II of the criterion does not appear to apply at all for kernels with large standard deviation, and in fact wave trains satisfying the criterion in this case do not even exist for small B (the case II curve intercepts the Hopf bifurcation curve, ex. σ = 5 or 6). This is not surprising, since there is no reason to expect the frequency (4.16) to continue to play a role in selection when the coexistence state is not approached. In terms of the stability regions shown in Figure 4.6, we note that for all σ here the selected wave trains are linearly stable for small B, but unstable for B near the kinetic Hopf bifurcation. In fact, however, in nearly all cases even these unstable wave trains appear regular in the final time-series of our simulations, and we discuss this further in Section 4.5. A final trend of note is that, at least for regions near the local kinetic Hopf bifurcation the wavelengths of wave trains, both observed and predicted, decrease with increasing standard deviation of the non-local competition kernel. We find 130 350 350 300 300 300 250 250 250 200 200 200 150 150 150 100 100 100 L 350 (a) local 50 (b) σ = 2.0 50 0 0.005 0.01 0.015 0.02 0.025 (c) σ = 3.0 50 0 0 0.005 0.01 0.015 0.02 0.025 0.005 350 350 300 300 300 250 250 250 200 200 200 150 150 150 100 100 100 0.015 0.02 0.025 L 350 0.01 (d) σ = 4.0 50 (e) σ = 5.0 50 0 0 0.005 0.01 0.015 B 0.02 0.025 (f) σ = 6.0 50 0 0.005 0.01 0.015 B 0.02 0.025 0.005 0.01 0.015 0.02 0.025 B Figure 4.6: Predicted and observed wave trains in B − L space. Horizontal axis is the kinetic parameter B. Vertical axis is the wavelength L. Solid lines indicate case I of the criterion, dashed lines indicate case II. Point symbols are wavelengths observed behind simulated invasions. The shaded region indicates stable wave trains. 131 this interesting because it is in contrast to observations of stationary patterns formed behind invasion; in [15] the wavelength of stationary patterns formed behind invasion increased with the standard deviation of the kernel. Turing-Hopf instability and failure for large σ For the non-local competition standard deviations shown in Figure 4.6 the revised selection criterion continues to perform well for B close to the Hopf bifurcation in the kinetics. We did not perform simulations for even greater values of σ for two reasons. The first is that the study [15] found that large standard deviation resulted in very irregular behaviour behind the primary invasion, and so we would be unable to measure wave train parameters in this case. A second reason is that the revised selection criterion cannot provide a prediction for σ > σ ∗ , where σ ∗ ≈ 6.47, because it is no longer single-valued for given B. An example of the criterion for σ = 7.0 is shown in Figure 4.7. As discussed in [15], for σ > σ ∗ the primary instability of the coexistence steady state is no longer a Hopf instability, but a Turing-Hopf instability. In this case, wave trains may exist even for B > 0.02268 and may be limited to a band of wavelengths L rather than for all wavelengths above a minimum, as shown in Figures 4.3 and 4.7. Our numerical study suggests that the selection criterion becomes multi-valued at this same value, σ ∗ , and so we propose that it can only be expected to apply when the primary instability of the coexistence state is the Hopf instability. 4.5 Discussion In Section 4.4 we studied the selection of wave trains behind predator invasions in the non-local system (4.4) that proved to provide quite accurate predictions for a substantial range of the kinetic parameter and the non-local kernel standard deviation. In addition, in Section 4.3 we studied the stability of wave train solutions and were able to divide up the wave train family into stable and unstable parts. It is perhaps of greater interest to biologists 132 200 L 150 100 50 0 0.019 0.02 0.021 0.022 B 0.023 0.024 0.025 Figure 4.7: Selection criterion for σ = 7.0. Horizontal axis is the kinetic parameter B. Vertical axis is the wavelength L. Solid curve indicates case I of the selection criterion. Dashed line is the Hopf bifurcation curve. The shaded region indicates stable wave trains. 133 to know for what parameter values the selected wave train is stable. We can answer this question by coupling these two components together and solving for the intersection of the stability boundary with the selection criterion (4.12). We can do this by simultaneously continuing solutions to (4.7), (4.9), (4.10) and (4.11) with λ = 0, λ|| = 0 subject to the constraint (4.12). Using this method we have a prediction for the stability of wave trains selected behind invasion in B − σ space. Figure 4.8 shows the boundary of stability for selected wave trains behind invasion. We note that here we compute this boundary only for case I of the selection criterion since case II does not in many cases provide a valid prediction. 6 Unstable 5 σ 4 3 Stable 2 1 0 0.005 0.01 0.015 0.02 B Figure 4.8: Eckhaus instability boundary for selected wave trains in B − σ space In this figure, the most noticeable effect is that larger values of σ have a smaller parameter region, in terms of the parameter B, with stable selected wave trains. Since unstable selected wave trains are believed to be a route to spatiotemporal chaos, this suggests that spatiotemporal chaos may be more common behind predator invasions when the prey species competes 134 non-locally over a large range. While here we have studied only a single predator-prey model, the destabilizing effect of the non-locality on selected wave trains is in keeping with the overall destabilizing effects observed in [15]. Further study on other models is needed, but we hypothesize that highly non-local predator-prey systems will generally exhibit more unstable wave trains and resulting irregular oscillations. The study here provides two hypotheses of predator-prey systems with non-local prey competition that are empirically verifiable. The first, mentioned above, is that spatiotemporal chaos may be more prevalent in systems with higher non-locality of prey competition. The second is that the wavelength of wave trains selected behind invasion decreases with increasing non-locality of the prey competition. Unfortunately, since the idea of nonlocal interactions is relatively recent, we know of no studies to date that can be used to verify these predictions. Empirical research in this area in the future would therefore be very beneficial. There are also some clear paths forward on the theoretical side, which we describe here. The stability results we have presented here are based on spectral stability on the infinite line [21]. Natural systems, however, are usually finite in extent and in this case a number of studies have shown that so-called “convective” instabilities may not manifest in simulations. The idea behind this theory is that perturbations may propagate in space as they grow, and thereby convect through the boundary before becoming significant [23]. The study of such convective instabilities has been used by Sherratt et al. (2009) to explain the appearance of “bands” of regular wavetrains such as that shown in Figure 4.1(d) [31]. In fact, since our simulations here are for relatively short times nearly all of such supposedly “unstable” wave trains remained regular throughout the time-series for our simulations, but we expect that for longer times these would break down into irregular oscillations farther back, as in Figure 4.1(d). This brings to mind the important question of the biological relevance of the results found here. Whether or not a given system will exhibit the patterns and instabilities described here will depend on natural parameters, including the size of the domain and the magnitude of noise in the system. 135 Here, we have deliberately avoided boundary effects by keeping the domain boundaries well ahead of the primary invasion front. It would be interesting to also consider the more biologically relevant scenario of a fixed boundary and studying the dynamics after the invasion reaches the boundary. This scenario has been considered for local reaction-diffusion systems by several authors [11, 17]. Notably, when the boundary conditions are no-flux, in the long term regular wave trains have been shown to relax to a spatially homogeneous temporal oscillation whereas irregular oscilaltions persist. With regard to the effect of noise, the only perturbations our system is subjected to are numerical round-off errors and error from the discretization of the system. Natural systems are likely to also have a much greater level of noise, and these larger perturbations can have either a stabilizing or a destabilizing effect on the system. One way to model “noise” is to allow a small spatial variation of parameters. This approach has been shown in local systems to promote the generation of spatial and spatiotemporal patterns as well as to lead to persistence of wave trains on finite domains [12, 24]. The study of this form of noise in the non-local system here would also be useful for understanding more realistic “noisy” natural systems. A final area of future work we propose could be the more rigorous study of wave train selection in non-local systems. The criterion applied here did prove quite accurate for the system studied, but it is purely a heuristic extension of the criterion for the local case and fails to explain wave train selection in the case of a Turing-Hopf instability. We have so far been unsuccessful in studying this problem in more generality, but a more rigorous study of the selection mechanism, even in a reduced or “toy” non-local system, would be extremely helpful in further understanding the generation of wave trains in the wake of predator invasions. 4.6 Conclusion We have computed here predicted selected wave trains following predator invasion in a system with non-local prey competition as well as the boundary for onset of Eckhaus instability for the family of wave trains for various val136 ues of the non-local competition kernel standard deviation. This knowledge allows us to better understand the dynamics observed following predator invasions in simulations, and has the potential to help understand spatiotemporal patterns observed in nature. However, there is a need for a great deal of future work on stability and selection in such non-local systems, with two main areas of focus. The first is a need for empirical work on non-local competition in general to validate the model and with regard to wave train generation specifically to test the hypotheses we make here. Second, there is still a great deal of room for theoretical work on predator-prey systems with non-local prey competition, particularly on the wave train selection mechanism, as it does not yet have a sound mathematical basis. Acknowledgements This work was partially supported by funding from the Natural Sciences and Engineering Research Council of Canada (NSERC). As well, S.M. would like to acknowledge the Pacific Institute for the Mathematical Sciences (PIMS) for funding through the IGTC fellowship program. We would also like to thank Jens Rademacher for suggesting computation of essential spectra by continuation and answering our many questions on implementing the methods of [19]. 137 References [1] I. S. Aranson and L. Kramer. The world of the complex GinzburgLandau equation. Rev. Mod. Phys., 74:99–143, 2002. [2] J. Billingham. Dynamics of a strongly nonlocal reaction-diffusion population model. Nonlinearity, 17:313–346, 2004. [3] G. Bordiougov and H. Engel. From trigger to phase waves and back again. Physica D, 215:25–37, 2006. [4] N. F. Britton. Aggregation and the competitive exclusion principle. J. Theor. Biol., 136:57–66, 1989. [5] E. J. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, and X. Wang. AUTO97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont). Concordia University, 1997. [6] S. A. Gourley. Travelling front solutions of a nonlocal fisher equation. J. Math. Biol., 41:272–284, 2000. [7] S. A. Gourley and N. F. Britton. A predator-prey reaction-diffusion system with nonlocal effects. J. Math. Biol., 34:297–333, 1996. [8] S. A. Gourley, M. A. J. Chaplain, and F. A. Davidson. Spatio-temporal pattern formation in a nonlocal reaction-diffusion equation. Dynamical Systems, 16:173–192, 2001. [9] E. Hernandez-Garcia, J. Vinals, R. Toral, and M. San Miguel. Fluctuations and pattern selection near an Eckhaus instability. Phys. Rev. Lett., 70:3576–3579, 1993. [10] B. Janiaud, A. Pumir, D. Bensimon, V. Croquette, H. Richter, and L. Kramer. The Eckhaus instability for traveling waves. Physica D, 55:269–286, 1992. 138 [11] A. L. Kay and J. A. Sherratt. On the persistence of spatiotemporal oscillations generated by invasion. IMA J. Appl. Math., 63:199–216, 1999. [12] A. L. Kay and J. A. Sherratt. Spatial noise stabilizes periodic wave patterns in oscillatory systems on finite domains. SIAM J. Appl. Math., 61:1013–1041, 2000. [13] N. Kopell and L. N. Howard. Plane wave solutions to reaction-diffusion equations. Stud. Appl. Math., 52:291–328, 1973. [14] L. Kramer and W. Zimmermann. On the Eckhaus instability for spatially periodic patterns. Physica D, 16:221–232, 1985. [15] S. M. Merchant and W. Nagata. Instabilities and spatiotemporal patterns behind predator invasions in systems with non-local prey competition. in preparation, 2009. [16] S. M. Merchant and W. Nagata. Wave train selection behind invasion fronts in reaction-diffusion predator-prey models. submitted, 2009. [17] S. V. Petrovskii and H. Malchow. Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics. Theor. Popul. Biol., 59:157–174, 2001. [18] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C: the Art of Scientific Computing. Cambridge University Press, second edition, 1992. [19] J. D. M. Rademacher, B. Sandstede, and A. Scheel. Computing absolute and essential spectra using continuation. Physica D, 229:166–183, 2007. [20] J. D. M. Rademacher and A. Scheel. Instabilities of wave trains and turing patterns. Int. J. Bif. Chaos, 17:2679–2691, 2007. [21] B. Sandstede. Stability of travelling waves. In B. Fiedler, editor, Handbook of Dynamical Systems, volume 2, pages 984–1055. Elsevier Science B. V., 2002. 139 [22] B. Sandstede and A. Scheel. Absolute versus convective instability of spiral waves. Phys. Rev. E, 62:7708–7714, 2000. [23] B. Sanstede and A. Scheel. Absolute and convective instabilities of waves on unbounded and large bounded domains. Physica D, 145:233– 277, 2000. [24] A. Sasaki. Clumped distribution by neighborhood competition. J. Theor. Biol., 186:415–430, 1997. [25] J. A. Sherratt. Unstable wavetrains and chaotic wakes in reactiondiffusion systems of λ-ω type. Physica D, 82:165–179, 1995. [26] J. A. Sherratt. Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave. Physica D, 117:145–166, 1998. [27] J. A. Sherratt. Periodic travelling waves in cyclic predator-prey systems. Ecol. Lett., 4:30–37, 2001. [28] J. A. Sherratt, B. T. Eagan, and M. A. Lewis. Oscillations and chaos behind predator-prey invasion: mathematical artifact or ecological reality? Phil. Trans. Roy. Soc. Lond. B, 352:21–38, 1997. [29] J. A. Sherratt, M. A. Lewis, and A. C. Fowler. Ecological chaos in the wake of invasion. Proc. Natl. Acad. Sci., 92:2524–2528, 1995. [30] J. A. Sherratt and M. J. Smith. Periodic travelling waves in cyclic populations: field studies and reaction-diffusion models. J. R. Soc. Interface, 5:483–505, 2008. [31] J. A. Sherratt, M. J. Smith, and J. D. M. Rademacher. Locating the transition from periodic oscillations to spatiotemporal chaos in the wake of invasion. Proc. Natl. Acad. Sci., 106:10890–10895, 2009. [32] M. J. Smith and J. A. Sherratt. The effects of unequal diffusion coefficients on periodic travelling waves in oscillatory reaction-diffusion systems. Physica D, 236:90–103, 2007. 140 [33] W. van Saarloos. Front propagation into unstable states. Phys. Rep., 386:29–222, 2003. [34] W. T. Vetterling, S. A. Teukolsky, W. H. Press, and B. P. Flannery. Numerical Recipes Example Book (C). Cambridge University Press, second edition, 1992. [35] D. D. Wackerly, W. Mendenhall, and R. L. Scheaffer. Mathematical Statistics with Applications. Duxbury Press, fifth edition, 1996. 141 Chapter 5 Conclusion 5.1 Summary Since biological invasions are nowadays so frequent and nearly all ecosystems on the planet are at risk for the introduction of exotic species, it is important to be able to predict the effects of invasions in such systems [23]. Mathematical models have been shown to play a vital role in understanding the spread of invasive species and a classical modelling framework for this is reaction-diffusion systems [8]. A recent focus has been on so-called oscillatory reaction-diffusion systems, because these have been shown to generate complex spatiotemporal patterns behind predator invasions [15, 18–20]. Of course, these models are quite idealized in that they do not take into account a lot of biological reality such as interactions with other species, spatial and temporal variation of the environment, age structure and finite population sizes. Unlike the actual systems, however, they are simple enough to study mathematically, and so we can first attempt to understand these models with the goal of including more realism later on. In addition, even these greatly simplified models may help us identify and understand the characteristics of prey and predator species that give rise to spatiotemporal patterns, for instance oscillations of large amplitude and spatiotemporal chaos, that pose significant management challenges. This thesis therefore considers predator invasions in oscillatory reactiondiffusion predator-prey models as well as the non-local extensions of Chapters 3 and 4. The general goal is to understand, at least in these simple models, the generation of spatiotemporal patterns behind the invasion. The work of this thesis contributes to this understanding in three major areas: (1) wave train selection, (2) the influence of non-local prey competition and 142 (3) wave train stability. I now discuss each of these contributions in turn. 5.1.1 Wave train selection As outlined in the introduction, behind predator invasions somehow a single wave train solution is selected from a continuous one-parameter family of such solutions. Sherratt (1998) developed an approximation for the selected wave train using the λ − ω system that is valid near the supercritical Hopf bifurcation in the kinetics and when prey and predator diffuse at equal rates [14]. However, this prediction loses accuracy when these conditions are not met, and the development of a more accurate criterion was the objective of the study in Chapter 2. In this study we were indeed able to develop a more general wave train selection criterion for two-component reactiondiffusion predator-prey systems. We tested this criterion using the three oscillatory reaction-diffusion systems in the introduction and found that it retains accuracy farther away from the Hopf bifurcation and also when predator and prey diffusion rates are quite different. The selection of wave trains behind predator invasion was also studied for the system with non-local prey competition in Chapter 4. In a part of this study, we generalized the selection criterion from Chapter 2 to such non-local systems. This new criterion was also found to be more accurate for local reaction-diffusion systems with unequal diffusion coefficients than the criterion of Chapter 2. 5.1.2 The influence of non-local prey competition Non-local interactions have been considered in a range of modelling frameworks and have been added to reaction-diffusion models for single-species, competitive and predator-prey systems [1, 2, 5–7, 12]. Although the specific results vary, a common theme in these studies is that non-local interactions can be destabilizing and under various conditions can lead to complex spatiotemporal patterns. To my knowledge, non-local interactions have not so far been considered for oscillatory predator-prey systems and the studies of Chapters 3 and 4 are the first concerning spatiotemporal patterns formed 143 behind predator invasions in such systems. The objective of the study in Chapter 3 was to determine if and how nonlocal intraspecific competition of the prey species will qualitatively change the spatiotemporal patterns that evolve in the wake of predator invasions. We also sought to establish the dependence and sensitivity of these patterns to the standard deviation and kurtosis of the non-local interaction kernel. Therefore, in this study we considered three commonly used kernel forms as well as developed a caricature system to describe the influence of the standard deviation and kurtosis of the non-local competition on the patterns observed. We found that new patterns, notably stationary spatially periodic states, can emerge behind the invasion. We also found that the instabilities of the coexistence state were indicative of the pattern observed behind the invasion. In addition, the caricature system illustrates how both large standard deviation and low kurtosis facilitate the formation of spatiotemporal patterns behind predator invasions. Small standard deviation and highly kurtotic kernels, however, give qualitatively similar dynamics to the local system, as would be expected. 5.1.3 Wave train stability For both the local and non-local systems studied here, wave trains selected behind predator invasions can be unstable. In this case, the wave train may be observed only transiently immediately behind the invasion and break up into irregular oscillations farther back. Careful numerical studies suggest that these irregular oscillations are in fact spatiotemporal chaos [9, 13, 19]. Since chaotic systems present obvious management challenges, a natural question to ask is when the wave train selected behind invasion will be unstable. Although analytical results on wave train stability do not exist for general oscillatory reaction-diffusion systems, methods have recently been developed to compute their stability numerically [10]. Smith and Sherratt (2007) have recently applied these methods to an oscillatory reactiondiffusion system and were able to divide the family of wave train solutions into stable and unstable parts [22]. Coupled with a wave train selection 144 criterion such as that developed in Chapter 2, this provides a prediction of when irregular oscillations may be expected behind invasion. One of the goals of the study in Chapter 4 was to extend the methods of Rademacher et al. (2007) and thereby numerically study the stability of wave train solutions for one of the non-local systems in the introduction [10]. In this chapter, we first studied the effect of the non-local prey competition on the stability properties of the family of wave train solutions by dividing the family into stable and unstable parts, as in [22]. We found that as the standard deviation of the prey competition kernel increased, the band of wavelengths corresponding to stable wave trains shrank. We then coupled this with the selection criterion also developed in this chapter to get predictions for when the wave train selected behind a predator invasion would be unstable. We found that increased non-locality (higher standard deviation kernel) of prey competition resulted in a greater kinetic parameter range that selected unstable wave trains, suggesting that spatiotemporal chaos might be expected more frequently for highly non-local systems. 5.2 Discussion and future directions The results of the studies in Chapters 2-4 allow for more accurate predictions of wave trains selected behind predator invasions in these systems and an understanding of the effects of including non-local prey competition. However, this by no means completes the study of spatiotemporal patterns behind predator invasions. There are still many questions remaining in this field, and some of the results found here require further investigation as well as pose a few new problems. In the following I describe the limitations of the work shown in this thesis and outline some future directions for some of the remaining open problems. Wave train selection The wave train selection criterion presented in Chapter 2 and generalized in Chapter 4 gives relatively accurate predictions of selected wave trains for 145 some parameter regions, but loses accuracy in all cases far from the Hopf bifurcation in the kinetics. The study of the complex Ginzburg-Landau system in Chapter 2 illustrates how the wave train selection problem can be considered as the selection in a 2-parameter family (parameterized by frequency and speed) of point-to-periodic heteroclinic connections from the spatially homogeneous steady state to a wave train. For more general systems, the criterion we developed used the linear spreading speed and associated frequency to select a member of this family. Since these are determined from linearization about the homogeneous state they can provide only an approximation to the selected wave train, and this is a likely source of error for the selection criterion. A possible way forward in this respect may lie in the numerical study of the family of point-to-periodic heteroclinic connections, rather than the wave trains and homogeneous steady states in isolation. Numerical methods have recently been devised for the computation of these solutions and I think a study of the selection for the full point-to-periodic solutions could allow for a greater understanding of the selection mechanism [3]. There are also two other situations in which the selection criterion developed in this thesis fails. The first case has been studied by Dunbar (1996) and is when a wave train of the same speed as the primary invasion front exists [4]. In this case the invasion instead takes place as a point-to-periodic heteroclinic connection in the single comoving frame of the primary invasion front. The second case has not yet been studied, and is when the onset to instability of the coexistence state is a Turing-Hopf instability, as for large σ in the non-local systems of Chapters 3 and 4. In this case the coexistence state has a non-zero leading linear unstable mode, and the selection criterion developed clearly fails because it becomes multi-valued. A study of this case, in particular in a simple reduced system as for the complex GinzburgLandau equation of Chapter 2, would be a large step towards understanding wave train selection more generally. Such a study would be a substantial undertaking, but a good place to start may be with the caricature system of Chapter 3. 146 Wave train stability The stability results discussed here in Chapter 4 as well as in Smith and Sherratt (2007) are for the spectral stability of wave train solutions on an infinite domain [22]. As recently studied by Sherratt et al. (2009), on finite domains so-called “convective” instabilities, where growing perturbations are transported towards a boundary, can lead to a band of regularappearing oscillations behind the invasion, followed by irregular oscillations farther back [11, 21]. The instability may not manifest for a long distance behind the invasion, and for short times the wave train appears regular, as noted in Chapter 4. A more useful indicator of irregular oscillations in this case may be whether the wave train is “absolutely” unstable, meaning that perturbations grow in time at every fixed point in the domain [11, 21]. This absolute stability can be determined by computing the so-called “absolute spectrum” of the corresponding operator, numerical methods for which are described in Rademacher et al. (2007) [10]. In fact, Sherratt et al. (2009) use this absolute spectrum to predict the width of regular wave train bands behind invasion in λ−ω systems [21]. Since the λ−ω system is an approximation to general oscillatory reaction-diffusion systems with nearly equal diffusion coefficients near a supercritical Hopf bifurcation in the kinetics, this helps us understand why we observe these regular wave train bands, and to predict their width when the approximation is valid. However, better results could be found by considering the absolute spectrum for wave train solutions of the full system. To my knowledge this has not yet been studied, although methods for this are outlined in [10]. This is likely because the geometry of the absolute spectrum is more complex and less well understood than for the essential spectrum, and in particular initial solutions for the continuation procedure can be difficult to find [10]. Still, I think that future work in this area could be productive. Non-local systems Since the studies of Chapters 3 and 4 are the first to consider spatiotemporal patterns behind predator invasions in these non-local systems, there are 147 many directions to continue this work in. Some are obvious extensions, such as spatiotemporal convolution for the prey competition, other choices for the local predator-prey kinetics, unequal diffusion coefficients (in particular stationary prey) and other kernel forms (especially for the selection and stability of wave trains). While studies of these problems would be helpful, there are two main directions that I think have the greatest potential for increasing our understanding of this topic. The first is the study of predator invasions in these models where the predator-prey interaction is modelled non-locally. Such a non-local interaction could be expected to influence the spatiotemporal pattern observed behind the invasion, but in addition may alter the rate of spread of the predator, since the speed of the predator invasion front is dependent on the growth rate of the predator at the leading edge. The study of this problem would ideally begin with the development of an appropriate form for the non-local predation interaction, similar to Gourley (2000) for non-local intraspecific prey competition [5]. There would then be a number of interesting problems to consider, including analysis of linear spreading speeds and front speeds, and selection and stability of patterns behind invasion. Second, wave trains are also known to be generated in oscillatory reactiondiffusion systems by boundaries, for instance Dirichlet boundary conditions of the form u = 0 or Robin conditions ∂u ∂x = δu, applied at either end of the domain and to both or only one of the two species. Either of these may be appropriate when the habitat outside of the spatial domain is hostile, and for the Robin condition the magnitude of δ (δ > 0) determines the degree of this hostility [20]. This mechanism for the generation of wave trains is mathematically quite different from generation by predator invasions, and it would be interesting to consider the influence of non-local interactions on this selection [16, 17]. Since studies of this kind are performed on finite domains, this problem would also involve a consideration of how to appropriately model non-local interactions on finite domains, a question of interest on its own. 148 Concluding remarks The studies in this thesis contribute much to our understanding of the generation and stability of spatiotemporal patterns in the wake of predator invasions, and as this discussion illustrates there are still many exciting directions in which to go. One major direction I have not yet mentioned is the validation of these theoretical models with empirical data. Since natural systems are far more complex than these models account for, it is important to establish when these models are adequate representations, and to determine the most promising directions for incorporating further biological reality. In fact, although there is evidence in natural populations of wave train patterns, I know of no definite evidence in nature of such solutions being generated by invasion. However, the statistical techniques for the analysis of spatiotemporal data and the identification of wave train patterns have only been recently developed. This is therefore a key time for collaboration to design and conduct empirical studies on spatiotemporal patterns behind invasions, and thereby inform future directions for the theoretical models studied here. 149 References [1] J. Billingham. Dynamics of a strongly nonlocal reaction-diffusion population model. Nonlinearity, 17:313–346, 2004. [2] N. F. Britton. Aggregation and the competitive exclusion principle. J. Theor. Biol., 136:57–66, 1989. [3] A. R. Champneys and B. Sandstede. Numerical computation of coherent structures. In B. Krauskopf, H. M. Osinga, and J. Galan-Vioque, editors, Numerical Continuation Methods for Dynamical Systems, pages 331–358. Springer, 2007. [4] S. R. Dunbar. Traveling waves in diffusive predator-prey equations: periodic orbits and point-to-periodic heteroclinic orbits. SIAM J. Appl. Math., pages 1057–1078, 1986. [5] S. A. Gourley. Travelling front solutions of a nonlocal Fisher equation. J. Math. Biol., 41:272–284, 2000. [6] S. A. Gourley and N. F. Britton. A predator-prey reaction-diffusion system with nonlocal effects. J. Math. Biol., 34:297–333, 1996. [7] S. A. Gourley, M. A. J. Chaplain, and F. A. Davidson. Spatio-temporal pattern formation in a nonlocal reaction-diffusion equation. Dynamical Systems, 16:173–192, 2001. [8] S. V. Petrovskii and B. L. Li. Exactly solvable models of biological invasion. Mathematical biology and medicine series. Chapman & Hall, 2006. [9] S. V. Petrovskii and H. Malchow. Wave of chaos: new mechanism of pattern formation in spatio-temporal population dynamics. Theor. Popul. Biol., 59:157–174, 2001. [10] J. D. M. Rademacher, B. Sandstede, and A. Scheel. Computing absolute and essential spectra using continuation. Physica D, 229:166–183, 2007. 150 [11] B. Sanstede and A. Scheel. Absolute and convective instabilities of waves on unbounded and large bounded domains. Physica D, 145:233– 277, 2000. [12] A. Sasaki. Clumped distribution by neighborhood competition. J. Theor. Biol., 186:415–430, 1997. [13] J. A. Sherratt. Unstable wavetrains and chaotic wakes in reactiondiffusion systems of λ − ω type. Physica D, 82:165–179, 1995. [14] J. A. Sherratt. Invading wave fronts and their oscillatory wakes are linked by a modulated travelling phase resetting wave. Physica D, 117:145–166, 1998. [15] J. A. Sherratt. Periodic travelling waves in cyclic predator-prey systems. Ecol. Lett., 4:30–37, 2001. [16] J. A. Sherratt. Periodic travelling wave selection by Dirichlet boundary conditions in oscillatory reaction-diffusion systems. SIAM J. Appl. Math., 63:1520–1538, 2003. [17] J. A. Sherratt. A comparison of periodic travelling wave generation by Robin and Dirichlet boundary conditions in oscillatory reactiondiffusion equations. IMA J. Appl. Math., 73:759–781, 2008. [18] J. A. Sherratt, B. T. Eagan, and M. A. Lewis. Oscillations and chaos behind predator-prey invasion: mathematical artifact or ecological reality? Phil. Trans. Roy. Soc. Lond. B, 352:21–38, 1997. [19] J. A. Sherratt, M. A. Lewis, and A. C. Fowler. Ecological chaos in the wake of invasion. Proc. Natl. Acad. Sci., 92:2524–2528, 1995. [20] J. A. Sherratt and M. J. Smith. Periodic travelling waves in cyclic populations: field studies and reaction-diffusion models. J. R. Soc. Interface, 5:483–505, 2008. 151 [21] J. A. Sherratt, M. J. Smith, and J. D. M. Rademacher. Locating the transition from periodic oscillations to spatiotemporal chaos in the wake of invasion. Proc. Natl. Acad. Sci., 106:10890–10895, 2009. [22] M. J. Smith and J. A. Sherratt. The effects of unequal diffusion coefficients on periodic travelling waves in oscillatory reaction-diffusion systems. Physica D, 236:90–103, 2007. [23] M. Williamson. Biological invasions, volume 15 of Population and community biology series. Chapman & Hall, 1996. 152
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Spatiotemporal patterns in mathematical models for...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Spatiotemporal patterns in mathematical models for predator invasions Merchant, Sandra M. 2009
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Spatiotemporal patterns in mathematical models for predator invasions |
Creator |
Merchant, Sandra M. |
Publisher | University of British Columbia |
Date Issued | 2009 |
Description | Much attention has been given to oscillatory reaction-diffusion predator-prey systems recently because, in the wake of predator invasions, they can exhibit complex spatiotemporal patterns, notably wave trains and associated irregular spatiotemporal oscillations, thought to occur in natural systems. This thesis considers the generation and stability of spatiotemporal patterns behind invasion in these models and an extension that includes non-local intraspecific prey competition. In the first part, we study the mechanism by which a single member is selected from a continuous family of wave train solutions behind the invasion. This was first studied by Sherratt (1998), where the author develops a selection criterion that is valid near a supercritical Hopf bifurcation in the kinetics and when the predator and prey diffuse at equal rates. We formulate a ``pacemaker" selection criterion that generalizes the criterion of Sherratt (1998), but does not depend on these assumptions. We test this pacemaker criterion on three sample systems and show that it provides a more accurate approximation and can apply to unequal diffusion coefficients. In the second part of the thesis, we study the effect of including non-local intraspecific prey competition in these systems. We first study the qualitative effect of non-local competition on the spatiotemporal patterns behind predator invasions in these models, and in a related caricature system. We find that non-local prey competition increases the parameter range for spatiotemporal pattern formation behind invasion, and that this effect is greater for lower kurtosis competition kernels. We also find that sufficiently non-local competition allows the formation of stationary spatially periodic patterns behind invasion. Second, we revisit the selection and stability of wave train solutions. We modify the selection criterion from the first part, also applying it to the non-local system, and study how the properties of selected wave trains vary with the standard deviation of the non-local prey competition kernel. We find that the wavelength of selected wave trains decreases with the standard deviation of the non-local kernel and also that unstable wave trains are selected for a larger parameter range, suggesting that spatiotemporal chaos may be more common in highly non-local systems. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0068997 |
URI | http://hdl.handle.net/2429/17988 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2010-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2010_spring_merchant_sandra.pdf [ 2.79MB ]
- Metadata
- JSON: 24-1.0068997.json
- JSON-LD: 24-1.0068997-ld.json
- RDF/XML (Pretty): 24-1.0068997-rdf.xml
- RDF/JSON: 24-1.0068997-rdf.json
- Turtle: 24-1.0068997-turtle.txt
- N-Triples: 24-1.0068997-rdf-ntriples.txt
- Original Record: 24-1.0068997-source.json
- Full Text
- 24-1.0068997-fulltext.txt
- Citation
- 24-1.0068997.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0068997/manifest