Analysis of an Integrodifference Model for Biological Invasions With a Quasi-Local Interaction by Sandra M . Merchant B . S c , The University of Northern British Columbia, 2001 A THESIS S U B M I T T E D IN PARTIAL F U L F I L M E N T OF THE R E Q U I R E M E N T S FOR THE D E G R E E OF M A S T E R OF SCIENCE in The Faculty of Graduate Studies (Department of Mathematics, Institute of Applied Mathematics) We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRUTISH C O L U M B I A September 7, 2003 © Sandra M . Merchant, 2003 In presenting this thesis in partial fulfilment of the requirements for an ad-vanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study, f further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics, Institute of Applied Mathematics The University Of British Columbia Vancouver, Canada Date Abstract 11 Abstract The behaviour of a new model for the spatial spread of biological invasions with non-overlapping synchronous generations and well-defined dispersal and seden-tary stages is examined. In this integrodifference model, competition between conspecifics takes the form of a quasi-local interaction, where the strength of competition between two individuals depends on their physical distance from each other. Both the deterministic model and a stochastic analogue are exam-ined by numerically simulating the spread of a localized initial population over several generations. By modelling intraspecific competition with a quasi-local interaction, the shape of the travelling waves changed significantly from that of the classical model with only local competition, creating more variable and com-plex wavefront shapes than are possible with the classical model. The addition of quasi-local competition was also found to alter several aspects of the initial be-haviour of this model, including the invasion speed and spatial structure, although in the deterministic case the asymptotic invasion speed and population density behind the front of the wave agreed with those of the classical model. In the stochastic analogue, however, the rate of spread of the invasion was found to be considerably lower than that of the classical model, both initially and asymptot-ically. Furthermore, the speed achieved by the stochastic invasions was found to depend on the parameters of the quasi-local interaction kernel. Contents 111 Contents Abstract i i Contents i i i List of Figures v Acknowledgements vii 1 Introduction to Invasion Models and Quasi-Local Interactions . . . 1 1.1 Historical Invasion Models 1 1.2 Travelling Waves and Invasion Speed 4 1.3 Quasi-Local Interactions 6 1.4 Invasion Models Including Quasi-Local Interactions: Open Ques-tions 7 2 Methods 11 2.1 Model System Chosen for Simulations 11 2.2 Deterministic Simulations 16 2.3 Stochastic Simulations 18 3 Results 22 3.1 Deterministic Simulations 22 3.1.1 Simulations Without a Quasi-Local Interaction 22 3.1.2 Simulations With a Quasi-Local Interaction 24 3.2 Stochastic Simulations 31 3.2.1 Simulations Without a Quasi-Local Interaction 31 3.2.2 Simulations With a Quasi-Local Interaction 34 4 Discussion and Conclusions 48 4.1 Summary and Discussion 48 4.2 Conclusion and Recommendations for Future Work 50 Contents iv Bibliography 52 A Source Code for Deterministic Simulations 55 B Source Code for Stochastic Simulations 70 List of Figures v List of Figures 3.1 Travelling Waves for the Deterministic Simulation Without a Quasi-Local Interaction 22 3.2 Wavefront Distance from Origin over Time for the Deterministic Simulation Without a Quasi-Local Interaction 23 3.3 Population Density Waves for the First 25 Generations of the De-terministic Simulation With a Quasi-Local Interaction 25 3.4 Simulation Invasion Speeds over Time for the Deterministic Sim-ulations 26 3.5 Simulation Invasion Speeds vs. Time for the First 30 Generations of the Deterministic Simulations 27 3.6 Density Waves at Generation 50 for the Deterministic Simulation With a Quasi-Local Interaction 28 3.7 Height of the Wavefront Peak over Time for the Deterministic Simulation With a Quasi-Local Interaction 29 3.8 Formation of a Wavefront Peak for the Deterministic Simulation With a Quasi-Local Interaction 30 3.9 Spread of a Single Population for a Stochastic Simulation Without a Quasi-Local Interaction 32 3.10 Travelling Waves for the Stochastic Simulation Without a Quasi-Local Interaction 34 3.11 Difference in Simulation Invasion Speeds from the Theoretical Deterministic Value for the Stochastic Simulation Without Quasi-Local Interactions at Various Mesh Sizes 35 3.12 Initial Invasion Speeds for the Stochastic Simulation With a Quasi-Local Interaction 41 3.13 Invasion Speed at each Generation for the Stochastic Simulation With a Quasi-Local Interaction 42 3.14 Average Location of the Furthest Individual from the Origin at each Generation for the Stochastic Simulation With a Quasi-Local Interaction 43 List of Figures vi 3.15 Average Survival Probability of the Furthest individual from the Origin at Generation 50 for the Stochastic Simulation With a Quasi-Local Interaction 44 3.16 Invasion Speed at Generation 50 for Various Mesh Sizes for the Stochastic Simulation With a Quasi-Local Interaction 45 3.17 Average Population Density Waves at Generations 100 and 200 for the Stochastic Simulation With a Quasi-Local Interaction . . . 46 3.18 Average Population Density Wave at Generation 100 for the Stochas-tic Simulation with a Quasi-Local Interaction 47 Acknowledgements Vll Acknowledgements I would like to thank my supervisor Michael Doebeli for his continual support throughout this project, as well as the entire Doebeli lab group for their many sug-gestions on the construction and analysis of these invasion models. In particular, I would like to thank Rik Blok for providing me with several technical suggestions and the source code for drawing Poisson distributed random deviates. I would also like to thank the Zoology 502 class for their review and suggestions on the introductory chapter. Finally, I would like to thank Cameron Rowe for his 24-hour caffeine provision and technical support and for having so much patience with me. Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 1 Chapter 1 Introduction to Invasion Models and Quasi-Local Interactions Over the past half century the mobility of human societies, and the scale over which we alter our environment, has changed dramatically. One of the conse-quences of this increased globalization is an increase in the frequency of biolog-ical invasions [28]. The introduction of an exotic species to a novel environment often has negative impacts on the persistence of local flora and fauna, resource productivity, and human health. Common textbook examples of these impacts include i) the introduction of the Nile Perch into Lake Victoria in East Africa, resulting in the extinction of many endemic cichlid species [25] ii) the introduc-tion of rabbits to Australia in 1859, causing millions of dollars in lost livestock production each year [25, 28] and iii) the intentional introduction of African bees to South America. The African bees hybridized with local species, producing the infamous "killer bees" that have now killed hundreds of people [25]. Since biological invasions can have such great impacts on humans and our en-vironment, models to predict the rate of spread, density of the population, and shape or pattern of an ecological invasion have been studied since the 1950s. Nonetheless, some key components required for their application, such as a means for incorporating empirically observed dispersal patterns and spatially-explicit density-dependent effects, have remained largely unresolved. This chapter pro-vides a brief review of the development of invasion models over the past 50 years, and illustrates how integrodifference equations and quasi-local interaction theory have been developed to address these two issues and create more versatile and applicable invasion models. 1.1 Historical Invasion Models The traditional way of modelling ecological invasions is through reaction-diffusion equations. The first application of reaction-diffusion equations to ecological inva-sions was by Skellam in 1951. In his paper, Skellam applied the reaction-diffusion Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 2 equation developed by Fisher (1937) to model the propagation of an advantageous allele, to the spread of muskrats across Europe [5, 23, 24]. Reaction-diffusion models combine a diffusion term that describes the movement of individuals with a growth term that describes the growth dynamics of a population. For example, the Fisher equation oW N d2N _ = r W ( 1__ ) + f l_ ( U ) combines Brownian random dispersal with logistic population growth. Here N(x,t) is the population density at the point x and at time t, r is the intrinsic rate of pop-ulation growth, K is the carrying capacity and D is a constant of diffusion [8]. Since Skellam's use of the most simple reaction-diffusion model for a biolog-ical invasion, reaction-diffusion invasion models have been modified to include density-dependent mortality (as in the Fisher equation above), heterogeneous en-vironments, allee effects and multi-species interactions, making them quite versa-tile and well-understood [8]. However, one major limitation of reaction-diffusion models is that they inherently assume that individuals disperse randomly and con-tinuously, resulting in a distribution of dispersal distances (from point of origin and over a fixed amount of time) that is normally-distributed [3]. In contrast, dis-persal during actual invasions may take many forms and it soon became desirable to develop models that could incorporate these other dispersal patterns [1, 14, 23]. In an effort to satisfy this need, new model formulations incorporating various dispersal distributions were developed and examined. One such model type, de-veloped by Kendall (1965) and Mollison (1972, 1977) to model the spread of epidemics, is especially effective at incorporating non-random dispersal assump-tions, because the dispersal kernel is an explicit part of the model formulation. These models are called spatial contact models [9, 18, 19, 20]. Although Kendall and Mollison's original papers included only continuous-time spatial contact models, both continuous-time and discrete-time spatial con-tact models are now common in the literature [9, 18, 19, 20]. This paper will focus on discrete-time contact models, also known as integrodifference equations. Like reaction-diffusion equations, integrodifference models have their origins in population genetics [26, 27], and they have only recently been applied to popula-tion ecology [11]. In population ecology, integrodifference equations are used to model populations of organisms with non-overlapping synchronous generations and well-defined dispersal and sedentary stages. A typical application of such models is to invasive annual plant populations [1, 6]. Al l movement is assumed to take place during the dispersal stage, and population growth occurs only during Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 3 the sedentary stage. In integrodifference models, the sedentary/reproductive stage is typically modelled by a standard difference equation Nl+l=f(N,) (1.2) such as the general fitness equation introduced by Maynard Smith and Slatkin, * '+' = T W ° ' 3 ) or the more simple Beverton-Holt stock-recruitment relationship « - = T ^ which is simply the case b = 1 in the Maynard Smith equation [2]. In these equa-tions N, is the population density at time step or generation /, X is the intrinsic rate of growth in the absence of competition, b defines the form and strength of intraspecific competition and a scales the carrying capacity, which in biological terms is the maximum sustainable population density. [2, 4]. As it is currently written, equation (1.2) does not explicitly contain a spatial component, and so it cannot allow for dispersal. In order to address this, we instead consider N, (x) to represent the population density in generation / at the point x in space. Note that x will be a vector if more than one spatial dimension is under consideration. If we assume that the environment is homogeneous and that there is no dispersal, population growth can now be described by the relationship Nl+i(x)=f(Nl(x)). (1.5) If one wanted to include environmental heterogeneity in this model, one could do so by making the function f depend explicitly on the point in space and time, ie., f(x,t,Ni(x)), but for the purpose of this paper i will assume that resource availability is constant across space and time and so / will be a function of N,(x) alone. Next, to include dispersal in the model, a dispersal kernel, k(x,y), that de-scribes the dispersal pattern of the population is introduced. More precisely, k(x,y) is the probability density function for propagules dispersing from the point y, and so k(x,y)dy is the probability of an individual dispersing from a small neighbourhood of length dy about y to a neighbourhood of equal length about x [6, 11, 12]. The population density at x after dispersal is therefore the integral of the post-reproduction population densities at all points in space multiplied by the Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 4 probability density function value at each point. The resulting integrodifference equation is Nl+](x) = J k(x,y)f(Nt(y))dy (1.6) where the domain of integration is the domain of the probability density function k(x,y). As can be seen in (1.6), the dispersal kernel k(x,y) is an explicit part of this model, and so it is possible to incorporate knowledge of the dispersal pattern of an organism into integrodifference invasion models. However, a few comments are in order regarding the dispersal kernel k(x,y). First, i f we assume that all individuals survive the dispersal process, then all individuals must go somewhere and it is necessarily true that V x, J k(x,y)dy = 1. As well, since k(x,y) is a probability density function Vx,y , k(x,y)>0. Finally, although such dispersal kernels may depend on both the points x and y, it is common to assume that the probability of an individual dispersing from a point x to the point v i s a function only of the distance between them. That is, k(x,y) = k(\x-y |). In keeping with the assumption that the environment is homogeneous, I wi l l also assume this form of dispersal kernel in this paper. We now have the integrod-ifference equation Nt+l(x) = J k(\x-y\)f(Nt(y))dy, (1.7) which is the starting point for many invasion models [1, 11, 12, 13, 14, 15, 16,17, 26]. 1.2 Travelling Waves and Invasion Speed There are a number of results from various analyses of integrodifference equations of the form (1.7), perhaps the most well-known, and most ecologically relevant of Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 5 which relates to the expected rate of spread of a population into a new environ-ment. It is known from Henry Weinberger's (1974, 1984) studies of integrodiffer-ence equations in population genetics that if (1.7) satisfies: (i) f'(N) > 0 V7V>0 (ii) f(N) < f(0)N VN > 0, and (iii) jfc(£) = fk(\z\)(&dz exists for some interval of E, about 0 (Note that the 3rd condition is that the moment generating function of k(\z\) must exist, which is equivalent to requiring that k(\z\) have exponentially bounded tails). then there will exist travelling wave solutions to the integrodifference equation (1.7) (see Figure 3.1 for an example of travelling waves) [26, 27]. These travelling waves have minimum wave speed (ie. distance between succesive wavefronts) c* = mm{l-ln[f(0)k(£,)}}- (1-8) Furthermore, if the initial data NQ(X) are compact, then c* is the asymptotic speed of propagation of the population. This means that the spread of the population will have speed tending asymptotically toward c*. This result is relevant to eco-logical invasion models because most population growth functions without allee effects satisfy conditions (i) and (ii) and many dispersal kernels, Gaussian dis-persal kernels in particular, have exponentially bounded tails. As well, founder populations always have bounded spatial support. The asymptotic speed of prop-agation defined by equation (1.8) has been calculated for a number of dispersal kernel forms. For instance, if the dispersal kernel is a Gaussian NormaliO^j) density function, hereafter abbreviated N(0, Oj), the asymptotic invasion speed is known to be c* = Gcj\/2ln(f'(0)) [12]. Interestingly, reaction-diffusion models also predict the same constant rate of spread for a population undergoing diffusive dispersal and with the same conditions (i) and (ii) on the growth function [8]. Of course, not all dispersal kernels have moment generating functions. In such a case, Kot et al. (1996) have outlined two possibilities for the invasion speed: 1. The invasion speed grows polynomially with time: This occurs when the dispersal kernels have "fat tails" or more precisely k(\z\) has finite moments of all orders but no moment generating function. Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 6 2. The invasion speed grows exponentially with time: This occurs when the dispersal kernels have extremely fat tails. That is, k(\z\) has moments that are infinite [12]. 1.3 Quasi-Local Interactions Quasi-local interactions, also known as neighbourhood effects, are a recent ad-dition to spatial population and metapopulation theory. Sasaki (1997) provided the first analysis of neighbourhood effects in continuous-space population models, and Doebeli & Killingback (2003) developed the concept for metapopulation/discrete-space models [4, 22]. Recall from section 1.1 that integrodifference models are built around an assumed growth function of the form Nl+l(x)=f{N,(x)). Here population growth at a point x is assumed to depend only on the popula-tion density at that point. This means that only individuals at x have an influence on the growth rate at x, whereas individuals at other coordinates have no impact whatsoever. For many organisms this is not a biologically plausible assumption, since neighbouring individuals can be expected to influence survival and/or repro-ductive output to some degree. In such a case, rather than as a delta function, the strength of competition between individuals may in fact better be described by a more gradual (usually decreasing) function of the distance between them. That is to say that individuals far from x can be expected to have less strong an im-pact on the reproductive output or mortality of individuals at x than would nearby conspecifics, but that intermediate degrees of competition are possible. Many pro-cesses would be expected to lead to such a dynamic, including: accumulation of waste products, antibiotics directed to conspecific competitors, behavioural inter-ference/harassment of neighbours, or nutrient absorption/foraging by neighbours [4, 22]. Since there are so many ways for quasi-local interactions to exist within a population, we can expect numerous organisms to exhibit them. To take into account these quasi-local interactions in spatial population mod-els, a competition function C(d), assumed to depend on the distance d between two points, is introduced. This function describes the strength of competition between individuals at a distance of d from each other. Although any form of function can be used, as with the dispersal kernel it is common to use standard probability density functions for this competition function, which I will hereafter refer to as the quasi-local interaction kernel. Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 7 If the quasi-local interaction kernel C(d) is known, we can then consider the effective population density, Nefftt(x) at the point x and in generation / to be the weighted sum (weighted by the quasi-local interaction kernel) of the population densities at all points in space: NeffM = Jc(\x-z\)N,(z)dz, (1.9) where the domain of integration is the domain of the quasi-local interaction kernel. We can then substitute this effective population density for the actual population density in the density-dependent component of the growth function of any spatial population model, and we would have a spatial model that incorporates intraspe-cific competition in the form of a quasi-local interaction. For example, in Sasaki's paper (1997) the classical Lotka-Volterra competition equation with diffusion is extended to a neighbourhood competition model using a Gaussian N(0,o) quasi-local interaction kernel, resulting in dN a r (*-y)2 d2N %-=R(x)N--—{ e~ ^ N(t,y)dy}-N + D— (1.10) at V2TCO" J OXL where D is the diffusion constant, a is the mortality due to unit population den-sity and R(x) is the per capita net growth rate at x [22]. (In this model resource availability varies in space.) With this model, Sasaki demonstrated a number of interesting consequences of adding quasi-local interactions to spatial population models, including how the addition of a quasi-local interaction could allow the pattern of spatial distribution to become strongly clumped when only a negligibly small variation in resource availability is present. In addition, when this model was applied to a population expanding in space, he was able to produce shapes for the travelling waves that were not possible with the basic model, including a quasi-periodic travelling wave with strongly clumped spatial structure stably maintained behind the front of the wave [22]. 1.4 Invasion Models Including Quasi-Local Interactions: Open Questions We could also add a quasi-local interaction to integrodifference models of inva-sion processes. If we replace N, (y) with the effective population density function Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 8 (1.9) evaluated at y into the density-dependent growth component of the integrod-ifference equation (1.7) we have a discrete-time spatial model #,+,(*) = / k{\x-y\)f{N,(y),Neff.t(y))dy (1.11) that describes the invasion of continuous space by a population of organisms that exhibits intraspecific competition in the form of a quasi-local interaction. The dynamics of this model need to be examined so as to determine whether adding quasi-local interactions to integrodifference invasion models can significantly al-ter the resulting dynamics. That is, do the salient features of such invasion models change when a quasi-local interaction is added? Based on known behaviours of the precursor models, several questions need to be addressed. The properties of this new model that I intend to examine fall into four categories: 1. Invasion Speed To ecologists, the rate of spread of an invading species is often the most critical element of a biological invasion [7, 23]. Since many invasions are undesirable, predictionsof the rate of spread can be valuable when devising containment plans, performing damage assessments, and for forecasting when the invasion will reach a particular location. In addition, the travelling wave solutions and the constant asymptotic speed of propagation for traditional integrodifference models also ren-ders the rate of spread of the invasion a feature of mathematical interest. The primary questions I would like to answer regarding the invasion speed of an inte-grodifference model including a quasi-local interaction are: • Does the addition of a quasi-local interaction speed up or slow down the invasion initially? • Is the asymptotic invasion speed different from that expected without a neighbourhood effect, and how does the invasion speed achieved relate to Weinberger's formula? 2. Wavefront Shape Although the shape of the wavefront is not often very important from the ecolog-ical perspective, particularly because real invasions do not usually create smooth and well-defined wavefronts, the wavefront shape is important mathematically be-cause changes in the shape of the wavefront in the deterministic model can often Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 9 reflect fundamental changes in the structure of the underlying model. For instance, in diffusion models, making the dispersal component density-dependent has been shown to significantly alter the shape of the wavefront [8]. 1 am therefore in-terested in examining changes in the shape of the wavefront when a quasi-local interaction is added to an integrodifference invasion model. I will attempt to an-swer the following questions about the wavefront shape: • Does the addition of a quasi-local interaction alter the shape of the travelling waves? • If so, how does the new wavefront shape depend on the parameters of the quasi-local interaction kernel? 3. Spatial Patterns The spatial pattern of biological invasions can be quite simple or complex, de-pending on the properties of the invading organism, the distribution of resources through space, and the distribution of interacting species (predators, competitors, pathogens and parasites). Moreover, the spatial pattern formed by an invasion is important to predict and observe because the level of impact the invader will have on the environment, humans, and other species will depend on the density of the invader. Invasion models to date have had relatively little success reproducing the complex spatial patterns formed by some biological invasions, and I am interested in seeing whether the addition of a quasi-local interaction to integrodifference models can produce more complex spatial patterns [16]. I therefore pose the fol-lowing questions regarding the effect of introducing a quasi-local interaction on the spatial pattern formed: • Do quasi-local interactions alter the spatial patterns formed by the invading population? For instance, is there more or less clumping and is the carrying capacity the same? • If the spatial pattern differs, what are the relationships between the parame-ters of quasi-local interaction kernels and spatial pattern formation? 4. Comparison with a Stochastic Individual-Based Analogue Although deterministic models such as this are typically used to predict spatial patterns and invasion speeds for ecological invasions, the dynamics observed in the deterministic model are not always present in empirical studies, or even in Chapter 1. Introduction to Invasion Models and Quasi-Local Interactions 10 stochastic and individual-based models. Therefore, this research will also attempt to determine whether a stochastic and individual-based model of this same system would display the same dynamics, and if it responds to the addition of a quasi-local interaction in the same way. In particular: • Is the speed of the invasion the same in the stochastic simulations as for the deterministic simulations? • Is the spatial pattern (carrying capacity and homogeneity) that is formed different when the simulation is stochastic and individual-based? Is the shape of the wavefront the same in the stochastic simulations as in the deterministic simulations? Chapter 2. Methods 11 Chapter 2 Methods To answer the questions posed in the previous section, I decided to principally examine this model by numerical simulations. There are a number of reasons for approaching the problem in this way. First, numerical simulations would al-low me to examine all of the properties of interest (invasion speed, wavefront shape, and spatial pattern). Secondly, it would provide me with the opportunity to examine the initial behaviour of the model, as well as approximate the asymp-totic behaviour. Third, simulating the population numerically would allow me to examine its behaviour as the parameters of the quasi-local interaction kernel are changed, making it possible to determine how the various features depend on these parameters. Finally, numerical simulations were necessary for constructing and analyzing a stochastic, individual-based analogue to the deterministic system. Two basic types of simulations were constructed: a purely deterministic inte-grodifference invasion process described by the recursion (1.11), and a stochastic individual-based system that will be described in detail below. I should note that to maximize computational efficiency all simulations were 1-dimensional in space. 2.1 Model System Chosen for Simulations For the numerical simulations, it was necessary to choose a particular functional form for the growth function, the dispersal kernel and the quasi-local interaction kernel. Because of its simplicity and its frequent use in biological models, the growth function I have chosen to use for my simulations is the Beverton-Holt stock-recruitment relationship [2]: W / N N XN,(x) + aN,(x) The parameters of this model are X > 0 which is the intrinsic rate of increase, and a > 0 which scales the carrying capacity (the carrying capacity/equilibrium den-sity is I should note here that although difference equations for population Chapter 2. Methods 12 growth can potentially exhibit complex dynamics, including periodic and chaotic behaviour, there is no complexity in the Beverton-Holt model. Under this model, provided the intrinsic rate of growth X is greater than 1, an isolated population wil l always grow or shrink monotonically towards the carrying capacity which is a globally stable equilibrium. Correspondingly, classical integrodifference mod-els with the growth component described by the Beverton-Holt model wil l have very stable dynamics and wil l not exhibit complexity in spatial pattern formation or the shape of the travelling waves. In addition, the Beverton-Holt function (without a quasi-local interaction) can easily be shown to satisfy the requirements (i) and (ii) for Weinberger's wave speed theorems: Verification of(i): f'(N) > 0 V N > 0. X-{\+aN)-XN-a (\+aN)2 (1+aN)2 > 0 V J V > 0 (since X>0) Verification of(ii): f[N) < f(0)N V N > 0. 1 + aN > 1 (since a > 0) => - L - < i 1+aN -XN < XN (since X > 0) => f(N) < XN / (AO < f(0)N (since f(0) = X) If we write the Beverton-Holt growth function in the form f(Nt(x)) =Nl(x)-g(Nt(x))=Nl(x) 1 \+aN,{x) Chapter 2. Methods 13 we can see that the density-dependent component, or in biological terms the per capita growth rate, of this function is g{N,(x)) = Therefore, to add a quasi-local interaction to the Beverton-Holt growth function we would simply replace A (^JC) in this function with the effective population density at x, Neffj{x), and so the form of the growth function that was used in the simulations including a quasi-local interaction is: AN,(x),NeW(x))- l+aNefL,(x)' The quasi-local interaction kernel chosen for the simulations was the N(0, cc) function l C(d) = e '2noc which leads to the following definition of the effective population density: NeffM = -^=~ r N ' ^ e 2 a I dy-Therefore, the effective population density Neffj{x) at a point x depends on a single parameter ac that indicates how rapidly the competitive effect of neighbours falls off with distance. A low value of oc would mean that the the impact of neighbours falls of quickly with distance, so only neighbours very nearby really matter. In contrast, a high ovvalue would indicate a more gradual decrease in impact with distance, so neighbours quite far away have nearly as great an effect as neighbours nearby. Of course, it should be noted here that the quasi-local interaction kernel is normalized in order to achieve the same effective population density when all local populations are at carrying capacity, regardless of the value of o c . It should also be noted that with this quasi-local interaction kernel the effective population density at a point x in a population that is not at carrying capacity everywhere will differ depending on the value of ac. For example, i f , . f 1 i f | x | < 1, N,(x)=l 1 ' 0 otherwise. Chapter 2. Methods 14 then, 1 f] v2 oc = l -> NeffJ{0) = -== / e'^dy = 0.6826 V27w-i rj c = 10 -> NefLt(0) = -= / e'^dy = 0.0796. V271 • 10 J-\ It is therefore expected that the behaviour of this model, and hence of the simula-tions, will vary as the width of the quasi-local interaction kernel, Gc, is varied. Any dispersal kernel with exponentially bounded tails satisfies condition (iii) for Weinberger's theorem on the asymptotic invasion speed. I chose a N(0,Oj) dispersal kernel for my simulations because it is the most frequently used in inva-sion models, including reaction-diffusion models. As well, it is easy to show that the N(0, Gd) dispersal kernel m 1 ~-2 '2KG d satisfies c o n d i t i o n ( i i i ) . Verification of (iii): k(t,) = f^kdz^e^dz exists for some interval oft, about 0 . %&) = r k{\z\)e^dz J —oo 1 - £ e 2%Gd 2^e^dz i r -+ttf-*?M e d dz 2%Gd • V((*-<&)2)+<#2) i r —— / e 2 c 2%Gd 1 e 2 71 kit,) = e+ which exists for t, G (—°°,°°). 2%Gd 1 irt n e 2 • y/% e~" -\f2Gddu - c o Chapter 2. Methods 15 From this simple form for k(E), which is just the moment generating func-tion for a N(0,Gd) distribution, we can see that the constant asymptotic speed of propagation 1 itt = M I ? < F + - p > >^o £ 2 for this integrodifference invasion model without a quasi-local interaction is quite easy to determine. If we define then Setting this equal to zero we have that °d giving cj* = "^g"^ as the only positive critical point. We can verify that there is a minimum at cj* by noting that 2ln{\) Thus, the asymptotic speed of propagation for this integrodifference system with a N(0,Od) dispersal kernel and without a quasi-local interaction is = Gdy/ln(k) | od^ln(X) c* = od^2ln{\). (2.1) - P - ( W = ^ > 0 . Chapter 2. Methods 16 2.2 Deterministic Simulations The goal of the deterministic simulations is to simulate as closely as possible the model including a quasi-local interaction with the function forms described above for a specified number of generations, as well as to estimate the asymptotic speed of propagation. Therefore, I wrote a C++ computer program to calculate the population densities via this recursion relation for a fixed number of generations and on a fixed grid of spatial coordinates. A copy of the source code for this program can be found in Appendix A. The gen-eral flow of the program is quite simple and is described below. Given the following: • JVo(jc): the initial (/ = 0) population densities at all spatial coordinates • n: the total number of generations • Oj: the standard deviation of the dispersal kernel • o"c: the standard deviation of the quasi-local interaction kernel • a and X: the parameters for the Beverton-Holt stock-recruitment relation-(2.2) ship Then, for each generation t = 0,1 ,2,3, . . . , n: for each spatial coordinate x: 1. Calculate the effective population density atx: 2. Calculate the growth/post-reproduction population density atx: f(Nt{x),Neff,{x)) XNt(x) 1 +a/V e / / ) / (x) Chapter 2. Methods 17 3. Calculate the population density atx after redistribution/movement: 1 , 0 0 _ ( ± Z ± £ Nl+l(x) = - = ^ e M f(N,(z),NeW{z))dz this is the population density at x in generation t+\. The composite trapezoidal rule rb m~l h (h — d) / f(x)dx^h-Jjf(a + ih) + -(f(a)+f(b)): h = — to approximate the integral of the function f(x) over the interval [a, b] using the function values at m + 1 equally-spaced points, was used in all cases for the nu-merical integrations required in steps 1 and 3 [10, 21]. Since this is intended to be a means for examining the features of interest and comparing them to the features of the model excluding the quasi-local interac-tion, I have also written routines to track the following variables throughout the simulations: • Distance of the wavefront from the origin at each time step: It should be noted that since the dispersal kernel has infinite tails, after the first gener-ation there is always a non-zero population density at every spatial coordi-nate. The "front" of the invasion is therefore considered to be the point fur-thest from the origin at which a preset, arbitrary density threshold has been exceeded. In all my simulations this density threshold was set to 0.001. • Invasion Speed for each value of ac: Ideally, the speed of the invasion at the end of a simulation is the difference between the location of the wavefront in the final generation, and its location in the previous generation. However, when mesh size is large, this can provide only a very coarse estimate of the invasion speed (for instance only to 1 decimal place for most of my simu-lations). Therefore, unless otherwise noted, invasion speed at the end of a simulation is therefore calculated as the mean displacement per generation of the invasion wave over the last half of the simulated generations. Simulations were also run assuming there was no quasi-local interaction, so that I could compare the behaviour to the simulations with a quasi-local interaction. In this case, the program above was modified to have 7Veyy,(x) = Nt(x) in step 1. Due to limited time and computing power, it was necessary to limit the scope of the deterministic simulations. The deterministic simulations were run with the following assumptions and limitations: Chapter 2. Methods 18 1. The parameters X and a in the Beverton-Holt growth function were fixed at X = 2 and a = 1 for most simulations and a retrospective sensitivity analysis was performed to ensure that changing these parameters would not qualita-tively alter the nature of any trends observed. 2. The parameter Gj was fixed at 1.0 for most simulations, and later runs with Gj = 0.5,0.75,1.25,1.5,1.75 and 2.0 were conducted to verify that the qual-itative behaviour of the model was independent of Gj. 3. Al l simulations were run for a minimum of 200 generations. Some simula-tions (o c = 1,10 and 100) were run for 500 generations in order to examine the persistence of initial trends. 4. The spatial grid over which the population densities were calculated and the integrations were performed was fixed with a mesh size (spacing between grid points) of 0.1 for most simulations. The effect of changing the mesh size was later examined by running several simulations with a mesh size of 0.01. 5. Only two different initial population density configurations were examined: (A) NQ(X) = 1.2 (greater than carrying capacity density) on [—0.1,0.1] and (B) No(x) = 0.1 (below carrying capacity density) on [—0.1,0.1]. 6. Finally, the simulations were run for 100 equally-spaced values of Gc, from Gc = 0.1 to Gc = 10.0, as well as the extreme case Gc = 100.0. 2.3 Stochastic Simulations The deterministic invasion model is somewhat unrealistic in a number of ways. First, only population densities are tracked, and it is assumed that population size is infinite and that population density is a continuous variable. However, most populations occur as discrete individuals and fractions of individuals cannot mi-grate, be born, or suffer from mortality. Also, the deterministic model incorporates no variation in individual behaviour and fitness. The addition of these factors to other integrodifference models has been shown to greatly alter the dynamics of the population [14, 16]. It is therefore of interest to determine how the effects of quasi-local interactions seen in the deterministic model will translate to a stochas-tic individual-based model. Chapter 2. Methods 19 To explore this question, I constructed a computer simulation, again written in C++, intended to be the stochastic and individual-based analogue of the deter-ministic simulations from section 2.2, again for a fixed number of generations and on a fixed grid of spatial coordinates. Note that in this case the grid points of the program correspond in some sense to "bins", with each bin corresponding to a section of the real line and the width of this section determined by the mesh size h. It is also important to note that this discretization of the real line means that distances between individuals are only approximations, with accuracy determined by the mesh size. In particular, individuals that move to only slightly different locations can end up in the same "bin" and therefore be treated as if they were at the exact same spatial coordinate. A copy of the source code for this program can be found in Appendix B. The general flow of the program is described below. Given the following: • NQ(X): the initial (/ = 0) population sizes at all spatial coordinates • n: the total number of generations • 0^ : the standard deviation of the dispersal kernel • G C : the standard deviation of the quasi-local interaction kernel • a and X: the parameters for the Beverton-Holt stock-recruitment relation-ship Then, for each generation t = 0,1,2,3,...,n: for each spatial coordinate/grid point x currently occupied: 1. Calculate the effective population atx: 1 (*-.y) 2 Nef/, (x) = h •'VNtiy) • —f=—e 2(j2 , h = mesh size y V 27t0c 2. Calculate the survival probability at x: Chapter 2. Methods 20 3. Determine how many of the individuals atx survive: - For each individual at x a random number is drawn from a Uniform [0,1] distribution, and compared to the survival probability at x. If it is greater than the survival probability the individual dies, otherwise it survives. 4. Calculate the number of offspring at x: - For each individual still alive at x, a random number is drawn from a Poisson distribution with a mean of X 5. Remove the parents and relocate the offspring: - The number of adults at all spatial coordinates is reset to zero. Then, for each offspring at x, a random number, z, is drawn from a JV(0, ov) distribu-tion and the number of new adults (kept distinct temporarily from the old adults so they do not get zeroed) at x + z (rounded to nearest mesh unit) is incremented by one. For the stochastic simulations 1 have tracked the following variables throughout the simulations: • Location of the furthest individual from the origin: This is used as a proxy for the "front" of the invasion wave. The reasons for choosing the furthest individual as a means to track the invasion wave are two-fold. First, it is extremely convenient computationally as it requires no other knowledge of population size or distribution. Secondly, from an applied perspective, the location of the furthest individual in the population is often seen as the lead-ing edge of the invasion. • Invasion speed for each value of o c : Invasion speed at the end of each stochastic simulation was also calculated as the mean displacement per gen-eration of the invasion wave's "front" over the last half of the simulated generations. • Population size at each time step: To identify whether the quasi-local inter-action affects population size. • Locations of all individuals at all time steps: To be used for comparing the average density at various time steps and spatial coordinates. • Whether a particular run/population went extinct before completion: To be used for comparing extinction rates between various simulations. Chapter 2. Methods 21 Since the purpose of the stochastic simulations was only to examine the dif-ferences between the deterministic model and its stochastic analogue, the scope of the stochastic simulations was kept quite limited. In particular, the following assumptions and limitations were made in the stochastic simulations I ran: 1. The parameters X and a in the Beverton-Holt growth function were fixed at A. = 2.0 and a = 1.0 for all simulations with a quasi-local interaction. 2. A l l simulations were started with a single individual at the origin. 3. The dispersal kernel standard deviation was fixed at O"^ = 1.0 for all simu-lations. 4. The performance of the stochastic simulations was found to be sensitive to the mesh size used. I therefore ran the stochastic simulations with mesh sizes of 0.1, 0.05, 0.025 and 0.01. 5. A l l simulations with a mesh size of 0.1 were run for 200 generations. Sim-ulations with a mesh size of 0.05 were run for 100 generations and all other simulations were run for 50 generations only due to processor time limita-tions. 6. As with the deterministic case, most simulations were run for, 100 equally-spaced values of o c , from o c = 0.1 to oc = 10.0. Due to the high computa-tion time required, however, it was only possible to conduct the mesh size 0.01 runs for o c = 0.1,1.1,2.1,3.1,4.1,5.1,6.1,7.1,8.1, and 9.1. As with the deterministic simulations, I constructed a version of the stochastic simulation without the quasi-local interaction, in this case by replacing the effec-tive population size at x calculated in step 1 with the actual number of individuals at the grid point x. Unfortunately, in the simulations without a quasi-local interac-tion, using the same parameter values in the growth function as was used for most of the deterministic simulations created a carrying capacity of 1 individual per grid point, resulting in very high mortality whenever two or more individuals dis-persed to the same grid point. Consequently a large proportion of the populations eventually went extinct. I therefore also ran the simulations that did not include a quasi-local interaction with an intrinsic growth rate of X = 6.0 and a mesh size of 0.01, as these two changes would greatly increase the population size at any given time, and therefore decrease the probability of extinction due to stochastic factors. Chapter 3. Results 22 Chapter 3 Results 3.1 Deterministic Simulations 3.1.1 Simulations Without a Quasi-Local Interaction The simulations with no quasi-local interaction behaved as predicted, producing a density wave spreading across space, as shown in Figure 3.1. 1.2 i 1 1 1 1 1 1 1 1 x ( s p a t i a l c o o r d i n a t e ) Figure 3.1: Travelling Waves for the Deterministic Simulation Without a Quasi-Local Interaction. These travelling waves are produced by an in-tegrodifference equation with Beverton-Holt local dynamics (a = 1, X = 2), initial condition (B) and a N(0,1) dispersal kernel. The den-sity curves shown are at 5 generation intervals. Chapter 3. Results 23 (a ) /= 0 - 5 0 0 (b) 7 = 0 - 2 0 Generation Generation Figure 3.2: Wavefront Distance from Origin over Time for the Deterministic Sim-ulation Without a Quasi-Local Interaction. Produced by the determin-istic simulation of the integrodifference equation with Beverton-Holt local dynamics (a = 1, X = 2), initial condition (A) and a 7V(0,1) dis-persal kernel. There are a number of features to note about the behaviour of this simulation, and consequently of the underlying model. First, we can see that after only a few generations the population density near the origin reaches the carrying capac-ity, which is in the Beverton-Holt equation. Once carrying capacity at the origin is reached, the density wave splits into two symmetric waves moving in opposite directions, behind each of which the density remains at carrying capac-ity thereafter. It is easy to see from Figure 3.1 why these waves are referred to as "travelling waves", since the wavefront shape remains constant through time and is simply translated along the x-axis each generation. We can also see from this figure that, after some initial redistribution, the distance between the density waves at equally spaced time intervals is approximately the same. The distance between two successive density waves is in fact the "asymptotic speed of propa-gation" referred to in Weinberger's theorem, and can be thought of as the invasion speed in biological invasion models. For the deterministic simulations, the distance between successive waves did indeed converge asymptotically toward a constant value, as can be seen in Figure 3.2 (a). In Figure 3.2 (b), we can also see that the invasion speed is initially faster than the asymptotic invasion speed in these simulations, causing the graph to be slightly concave down. This is due to the initial population density surrounding the origin Chapter 3. Results 24 being dispersed, and immediately causing many nearby grid points to exceed the threshold density limit. The result is that the invasion appears to be spreading more quickly initially and slowing down as it progresses. The asymptotic inva-sion speed was approached from below for all simulations without a quasi-local interaction, and the speed approached does appear to be the same as that theoreti-cally predicted by the equation (2.1) for each of these simulations. For instance, in the simulation with a = 1, X = 2, initial conditions (A), a 7V(0,1) dispersal kernel and a mesh size of 0.1, after 500 generations the calculated average distance the wavefront moved per time step was within 0.005 units of the value predicted by (2.1). 3.1.2 Simulations With a Quasi-Local Interaction The deterministic simulations that included a quasi-local interaction also produced symmetric density waves expanding out from the origin. Simulations with a stan-dard deviation of the quasi-local interaction kernel that was less than 1.2 times that of the dispersal kernel had a wavefront shape that was qualitatively indistinguish-able from that produced by the simulations without a quasi-local interaction. Sim-ulations with o c > 1.2 • Gj, however, produced travelling waves with a distinctly different shape, and correspondingly differed in spatial pattern and initial invasion speeds. Instead of the travelling waves shown in Figure 3.1, simulations with pa-rameter values that exceeded this threshold of ^ « 1.2 produced travelling waves with a shape similar to that of the quasi-periodic waves found by Sasaki (1997) [22]. A typical example of such travelling waves is shown in Figure 3.3. Similarly to my threshold, Sasaki (1997) found a more precise threshold of « 1.14 between the ratio of the standard deviation of the quasi-local interaction kernel and of the mean lifetime diffusion distance. If his parameters exceeded this ratio, Sasaki also found that the population expanded through space as a quasi-periodic travelling wave rather than as a simple monotonic wavefront [22]. From Figure 3.3 we can see that as well as the difference in wavefront shape, there are a number of ways in which the behaviour of this model differs from that without a quasi-local interaction. I wil l address each difference in turn, as well as discuss how the various properties change depending on the width of the quasi-local interaction kernel. Chapter 3. Results 25 1.6 I 1 1 r -30 -20 -10 0 10 20 30 x (spatial coordinate) Figure 3.3: Population Density Waves for the First 25 Generations of the Deter-ministic Simulation With a Quasi-Local Interaction. Produced by the deterministic simulation of the integrodifference equation with Beverton-Holt local dynamics (a = I, X = 2), initial condition (B), a7V(0,1) dispersal kernel and a 7V(0,5) quasi-local interaction kernel. 1. Invasion Speed For all values of o"c, and for both initial conditions, the invasion speed for simula-tions including a quasi-local interaction still asymptotically approached a constant value. In addition, as shown in Figure 3.4, for all o c values and for both initial con-ditions, the constant invasion speed approached appears to be the same constant as that approached by the deterministic simulations without a quasi-local interaction. Although it does not appear in this figure that the constant approached is exactly equal to that predicted by Weinberger's theorem for the model without the quasi-local interaction, the difference is most likely due to the limitation in precision from using a mesh size of 0.1 and averaging over only 250 wavefront translations. We can see in Figure 3.5, however, that for some values of a c , oc > 1.2 • 0^ in my simulations, the speed of the invasion actually differs considerably from that Chapter 3. Results 26 1.18 Figure 3.4: Simulation Invasion Speeds over Time for the Deterministic Simu-lations. Produced by the deterministic simulation of the integrodif-ference equation with Beverton-Holt local dynamics (a = 1, X = 2), initial conditions (A), a 7V(0,1) dispersal kernel and a N(0,oc) quasi-local interaction kernel. The simulation invasion speeds are calculated as the average distance between successive wavefronts for the previ-ous 250 generations. of the simulation with no quasi-local interaction for several generations near the beginning of the invasion. For these values of o c the rate at which new territory is invaded is lower than that of simulations with ac < 1.2 • Oj or without a quasi-local interaction. The source of this difference in initial invasion speeds for different 0 C values is the differential effect of the quasi-local interaction kernel in the initial generations of the invasion. In these initial generations, very little total territory is significantly occupied, and populations with a "wide" quasi-local interaction kernel, ie. a high value of o c , benefit from the empty space outside of the invasion wave. The result is that population density increases much more quickly near the origin when o c is large than when oc is small. For values of 0 C > 1.2 • Oj this increased initial growth causes the population density near the origin to overshoot Chapter 3. Results 27 the carrying capacity, as shown in Figure 3.3, with larger values of o 6 correspond-ing to a greater degree of overcompensation. This overshooting effect in fact negatively impacts population growth in the tails of the invasion wave, and results in decreased invasion speeds for larger o"c values. As we can also see in Figure 3.3, the increased population density near the origin eventually decreases to the expected carrying capacity as the invasion progresses and more territory is occu-pied, and the invasion speeds for all values of o c approach those of the simulation without a quasi-local interaction. Figure 3.5: Simulation Invasion Speeds vs. Time for the First 30 Generations of the Deterministic Simulations. Produced by the deterministic simu-lation of the integrodifference equation with Beverton-Holt local dy-namics (a — I, X = 2), initial conditions (B), a N(0,1) dispersal ker-nel and a 7V(0,oc) quasi-local interaction kernel. Note that for these simulations mesh size was 0.005 and in this figure invasion speed in generation / is calculated as the absolute difference between the loca-tion of the wave's "front" in generation / and its location in generation t - l . Chapter 3. Results 28 2. Wavefront Shape From Figure 3.3 we can see that the shape of the wavefront for the simulations with the quasi-local interaction with oc > 1.2 • ov differs from the shape of the wavefront when no quasi-local interaction is included or when G C < 1.2 Oj. Rather than a simple, Gaussian-shape monotonic descent from carrying capac-ity to zero population density, the density waves produced by these simulations have a damped oscillatory wavefront, which is more clearly seen in Figure 3.6. Figure 3.6: Density Waves at Generation 50 for the Deterministic Simulation with a Quasi-Local Interaction. Produced by the deterministic simulation of the integrodifference equation with Beverton-Holt local dynamics (a = 1, X = 2), initial condition (A), a N(0,1) dispersal kernel and a N(0,oc) quasi-local interaction kernel. Also, the amplitude of the oscillations at the wavefront differs considerably de-pending on the value of oc. As we can see in Figure 3.6, larger values of o c result in higher amplitude oscillations at the wavefront. With regards to the long-term behaviour of this model, it should be noted that the oscillations produced by adding a quasi-local interaction to the simulation are Chapter 3. Results 29 maintained, after some initial variation, from generation to generation and are not damped in time. The asymptotic height of the wavefront peak depends on the width of the quasi-local interaction kernel, with smaller values of 0 C correspond-ing to lower wavefront peaks. Figure 3.7 illustrates the fact that peak height does settle to a constant value, as well as how the asymptotic peak height depends on O) CD X CO CD 0_ CO CD Q E x CO 3.5 3 2.5 2 1.5 1 0.5 o c < 1.2 o c = 2.5 a c = 5.0 o c = 7.5 c c = 10.0 50 100 G e n e r a t i o n 150 200 Figure 3.7: Height of the Wavefront Peak over Time for the Deterministic Simu-lation With a Quasi-Local Interaction. Produced by the deterministic simulation of the integrodifference equation with Beverton-Holt lo-cal dynamics (a = 1, X = 2), initial condition (A), a7V(0,1) dispersal kernel and a 7V(0, o"c) quasi-local interaction kernel. Heuristically, it is not difficult to understand why this altered wavefront shape appears when oc is large relative to c</. If the wavefront shape was the same monotonic decrease as when there is no quasi-local interaction then coordinates near the top edge of wavefront would have an effective population density much lower than the actual density at these coordinates, as shown in Figure 3.8 (a). This decreased effective population density would result in a greater population growth at the coordinates near the top edge of the wavefront. Finally, since o"^ is small Chapter 3. Results 30 (a) (b) 140 Figure 3.8: Formation of a Wavefront Peak for the Deterministic Simulation With a Quasi-Local Interaction, (a) Actual population density and effective population density at generation 0. (b) Population density near the front of the invasion wave. Produced by the deterministic simulation of the integrodifference equation with Beverton-Holt local dynamics (a = 1, A, = 6), a JV(0,1) dispersal kernel, a JV(0,10) quasi-local in-teraction kernel, a mesh size of 0.1 and an initial No(x) population density distribution created by the 100th generation of the determinis-tic simulation without a quasi-local interaction. relative to o c the dispersal kernel will only redistribute this additional growth to the relatively local area, resulting in a greater population density in this area at the start of the next generation. Figure 3.8 (b) shows an example of the growth of the wavefront peak in a simulation that was started with an initial density curve created by the 100th generation of the deterministic simulation without a quasi-local interaction. As the wavefront peak forms, the increased population density at this peak causes the effective population density to increase, and after sufficiently many time steps the height of the peak stabilizes in the sense that the decrease in the effective population density from being near the edge is offset by the increased local population density. Of course, the damped oscillations in the wave shape are caused by the quasi-local interaction as well because this spatial coupling causes the effect of this peak of very high density to be felt by neighbouring points, resulting in lower growth at these points. This effect propagates itself down the invasion wave toward the origin, however it is dissipated in part by the dispersal kernel at each generation. It is therefore not surprising that the wavefront peak is higher and the oscillations have a higher amplitude when o c is increased. Chapter 3. Results 31 3. Spatial Pattern Since the wavefront shape is oscillatory i f a quasi-local interaction is included, it is apparent that the spatial pattern formed near the "front" of the invasion is a more patchy distribution than when no quasi-local interaction is added, with some areas having higher densities than others. Furthermore, this patchiness depends on the magnitude of oc, with a more patchy distribution when rj c is large. Addition-ally, since the wavefront shape was maintained in all simulations, the population density near the front of the wave would always be greater than inside the wave. Since in his model Sasaki (1997) was able to produce a quasi-periodic trav-elling wave with stably maintained and strongly clumped periodical structure be-hind the front of the wave, I thought it may be possible that a value of G C large enough could produce undamped oscillations in the density wave, and thereby a permanently maintained patchy distribution [22]. However, values of o c up to 100 were tested, and all produced oscillations that were spatially damped toward the origin. Therefore, when the simulations were run for enough generations that the wavefront was sufficiently distant, the density at any x-coordinate and for any value of <5C settled to the same constant carrying capacity as in the simulation with no quasi-local interaction. It is possible however that simulations were simply not run for sufficiently many generations, or the ratio between the parameters oc and was not sufficiently high to generate a permanently patchy distribution. 3.2 Stochastic Simulations 3.2.1 Simulations Without a Quasi-Local Interaction There are two ways in which I examined the stochastic simulations. The first method was to examine the behaviour of single runs of the simulation, since this represents the behaviour of a single invasion event under this model. A n example of a typical single simulation is shown in Figure 3.9. As can be seen in Figure 3.9, the population is spreading across space, and in a manner somewhat similar to the deterministic simulations. Two key differences should be noted. The first and most obvious difference is that only a discrete number of individuals exist at each spatial location, making both the wavefront's shape and location much more difficult to distinguish. Secondly, most spatial coordinates at any time step do not contain the number of individuals predicted by the carrying capacity for the Beverton-Holt equation. As we can see in Figure 3.9, even at grid points/spatial coordinates that are far behind the wavefront and Chapter 3. Results 32 a)/ = 10 • / = 20 c) -80 -60 -40 -20 0 20 40 60 80 x (spatial coordinate) / = 30 -80 -60 -40 -20 0 20 40 60 80 x (spatial coordinate) d) -80 -60 -40 -20 0 20 40 x (spatial coordinate) t = 40 -80 -60 -40 -20 0 20 40 60 80 x (spatial coordinate) Figure 3.9: Spread of a Single Population for a Stochastic Simulation Without a Quasi-Local Interaction. Produced by the stochastic simulation of the integrodifference equation with Beverton-Holt local dynamics (a = 1, X = 6), a/V(0,1) dispersal kernel and a mesh size of 0.01. Chapter 3. Results 33 therefore would be at carrying capacity in the deterministic case had as few as 0 individuals and as many as 17 in a single "bin" (carrying capacity in Figure 3.9 is 5 individuals per grid point). In addition, due to the relatively great degree of variability from simulation to simulation, it is difficult to ascertain whether differences between the deterministic and stochastic simulations are caused by fundamental differences in the invasion process, or if they are merely a reflection of using a non-representative stochastic run in the comparison. Since individual runs of the simulation are so difficult to analyze for similari-ties and differences in the shape, location and hence speed of the invasion waves, f also examined the average behaviour of a large number of runs of each type of simulation. The density graphs of the averaged population sizes would give a bet-ter defined invasion wave and using the mean behaviour of the stochastic process would give a more useful illustration of how the invasion would be expected to progress. Unless indicated otherwise, all average values for the stochastic simu-lations are the mean of 200 iterations of each stochastic simulation, and extinct populations are excluded from all averages. In fact, when the number of individuals at each spatial coordinate and at each time step is averaged over a large number of simulation runs, the resulting pop-ulation density distribution looks much more like the deterministic simulations without a quasi-local interaction. As can be seen in Figure 3.10, travelling density waves do appear for the stochastic simulations, and the wavefront shape and spa-tial pattern created by the stochastic simulations without a quasi-local interaction is in fact the same as that predicted by the deterministic simulations. In addition, the speed of the travelling waves, measured as the mean distance between the furthest individual from the origin in consecutive time steps, does in fact appear to approach a constant value, a result consistent with other stochas-tic models and simulations [14]. Interestingly, however, the constant value ap-proached for the stochastic simulations is significantly different from that pre-dicted by the invasion speed equation (2.1) and differs considerably also from the speed of the travelling waves in the deterministic simulations without a quasi-local interaction. Specifically, the invasion speed for the stochastic simulations is considerably lower than that predicted for the deterministic model, an attribute of stochastic models long conjectured and borne out in a number of other stochastic simula-tions [14, 20]. Also, in Figure 3.11 we can see that the difference in the stochastic simulation invasion speed from the deterministic expected value decreases as the mesh size is decreased. This is not surprising because decreasing the mesh size allows for more local populations per unit length, thereby creating both a greater Chapter 3. Results 34 x ( s p a t i a l c o o r d i n a t e ) Figure 3.10: Travelling Waves for the Stochastic Simulation Without a Quasi-Local Interaction. Produced by the stochastic simulation of the in-tegrodifference equation with Beverton-Holt local dynamics (a = 1, X = 6), a 7V(0,1) dispersal kernel and a mesh size of 0.01. The pop-ulation densities at each time step are calculated as the average of 1000 runs of the simulation (extinct populations are excluded from this average). number of individuals per unit length, and a greater total number of individuals in the population at any time step. Thus, increasing the mesh size brings the stochas-tic simulation closer in line to the deterministic assumption of infinite population size. 3.2.2 Simulations With a Quasi-Local Interaction As with the stochastic simulations that did not include a quasi-local interaction, I examined both the single run and expected/average behaviour of the stochastic simulations with a quasi-local interaction. In the single runs it was again difficult to find differences in the invasion wave shape because of the large amount of Chapter 3. Results 35 3 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.I h ( m e s h s i z e ) Figure 3.11: Difference in Simulation Invasion Speeds from the Theoretical De-terministic Value for the Stochastic Simulation Without Quasi-Local Interactions at Various Mesh Sizes. Produced by the stochastic sim-ulation of the integrodifference equation with Beverton-Holt local dynamics (a = 1, X = 6), a N(0,1) dispersal kernel and a mesh size of h. The invasion speed at each mesh size is calculated as the mean distance between furthest individuals in generations 50 to 100. variation in population size at each coordinate, but when a large number of runs were averaged for each value of oc, the resolution of the invasion waves improved. When examining this mean behaviour, the stochastic simulations that included a quasi-local interaction were found in many ways to be qualitatively quite similar to the corresponding deterministic simulations, but did differ in some way for each of the features of interest. 1. Invasion Speed As with the simulations without a quasi-local interaction, the most significant dif-ference in the behaviour of the stochastic simulations with a quasi-local interac-Chapter 3. Results 36 tion from that of the deterministic simulations is that the invasion speeds achieved by the stochastic simulations are significantly less than both those predicted by the equation (2.1) and those achieved in the deterministic simulations. When comparing the initial invasion speeds of stochastic simulations with dif-ferent values of 0 C , the invasion speeds achieved seem to follow a pattern not unlike that for the initial phases of the deterministic simulations. In particular, i f ac is greater than 1 the initial invasion speeds achieved decrease as oc is increased, just as in the deterministic case. We can see this similarity in initial pattern clearly in Figure 3.12(a) which shows the invasion speeds achieved during the first 30 generations for several values of 0 C . For a more detailed view of how the initial invasion speeds change with o"c, Figure 3.12(b) shows the relationship between oc and the invasion speed for the stochastic simulations after 50 generations. We can also see from Figure 3.12(a), that compared to the deterministic simu-lations, the differences in invasion speeds at different o"c values appear to be much greater in magnitude and to persist for many more generations. In fact, when ex-amining the long term invasion speeds of the stochastic simulations, as in Figure 3.13, I found that these differences in invasion speeds for various o"c values were still present even after 200 generations (the maximum number of generations for which the stochastic simulations were run). Furthermore, by examining the location of the wavefront as the invasion pro-gresses across the landscape, as in Figure 3.14, we can see that for each value of o"c the rate of spread does still appear to tend toward a constant value, and that for each value of 0"c there does indeed appear to be a different constant asymptotic rate of spread. (Note that the slope of the line for each value of ac is approximately equal to the invasion speed at that time step). The fact that the stochastic individual-based simulations achieve different asymp-totic invasion speeds depending on o"c is both a significant and interesting result. It is especially interesting because the deterministic simulations do not show this behaviour. Since the only stage of the stochastic simulation in which the value of o c is used is the survival stage, it is apparent that the difference in invasion speed at various 0 C values must be a reflection of differential survival rates at different oc values. As well, because the individuals very far from the wavefront have a very low probability of dispersing out past the wavefront into new territory, it is probable that the differential invasion speeds are the result of differential survival rates near the edge of the invasion wave. To investigate this possibility, I examined the survival rates at all points along the wavefront, and found that survival rates did indeed differ a significant amount and in the way expected for the various values of o c . These differences, however, Chapter 3. Results 37 Figure 3.12: Initial Invasion Speeds for the Stochastic Simulation With a Quasi-Local Interaction, (a) Invasion speeds for the first 30 generations at various o"6 values and (b) Invasion speeds after 50 generations for each oc value. Produced by the stochastic simulation of the inte-grodifference equation with Beverton-Holt local dynamics (a = 1, X = 2), a N(0,1) dispersal kernel and a A ^ O , ^ ) quasi-local interac-tion kernel. For (a) 5000 iterates of each simulation was performed and invasion speed at each generation was calculated as the mean difference in distance between furthest individuals in the current and previous generation. Chapter 3. Results 38 1.2 1.1 1 ac=1.0 oc=5.0 ac=10.0 o T5 CD CD CL OO 0.9 CO > o 0.4 0.3 20 40 60 80 100 120 140 160 180 200 G e n e r a t i o n Figure 3.13: Invasion Speed at each Generation for the Stochastic Simulation With a Quasi-Local Interaction. The distance travelled by the in-vasion wave in each generation is calculated as the mean difference in distance between furthest individuals in the current and previous generation. Produced by the stochastic simulation of the integrodif-ference equation with Beverton-Holt local dynamics (a = 1, X = 2), a N(0,1) dispersal kernel, a mesh size of 0.1 and a N(0, o"c) quasi-local interaction kernel. were only found in the very tails of the invasion waves, and were most prominent for the individuals furthest from the origin. Figure 3.15 shows the average survival rate for the furthest individual after 50 generations under various o c values. As we can see in this figure, the survival rate of the furthest individual follows a pattern very similar to that of the average invasion speed at generation 50 (shown in Figure 3.12(b)), which supports the conjecture that the invasion speed differences are in fact caused by differing survival rates near the wavefront edge. By examining the quasi-local interaction kernel used to calculate effective population densities, it is possible to see how differential survival for different G c values could come about in a finite population, and in particular in the tails of Chapter 3. Results 39 80 0 10 20 30 40 50 60 70 80 90 100 G e n e r a t i o n Figure 3.14: Average Location of the Furthest Individual from the Origin at each Generation for the Stochastic Simulation With a Quasi-Local Inter-action. Produced by the stochastic simulation of the integrodiffer-ence equation with Beverton-Holt local dynamics (a = 1, A. = 2), a N(0,1) dispersal kernel and a7V(0,o c) quasi-local interaction kernel. the invasion wave, where typically very few individuals are present. As a sim-ple example, consider the following population distribution on a grid of spatial coordinates with a mesh size of h: i f 1 i f x G [-100,100] o r x = 1 0 5 , 0 otherwise. Chapter 3. Results 40 0.9 0.85 CO "8 0.8 CO • | 0.75 =3 CO §> 0.7 CO i_ CD > < 0.65 + + . +++++-H-0.6 5 Or 10 Figure 3.15: Average Survival Probability of the Furthest individual from the Ori-gin at Generation 50 for the Stochastic Simulation With a Quasi-Local Interaction. Produced by the stochastic simulation of the in-tegrodifference equation with Beverton-Holt local dynamics (a = 1, X = 2), a N(0,1) dispersal kernel and a N(0, oc) quasi-local interac-tion kernel. then, ov = 0.1 Or = 1.0 ov = 50 AV/,,(105) 5(105) 1 271-0.1 •h 1 1 + 3.99/? 0.715 (if/? = 0.1) J W 1 0 5 ) 5(105) 1 271 •h 1 7Ve//,,(105) 1+0.399/; 1 0.962 (if A = 0.1) 5(105) « 0.66 (regardless of h). Chapter 3. Results 41 In essence, this simplification shows how individuals that disperse a long dis-tance in front of the remainder of the invasion wave have differing survival rates depending on o c , with the same pattern as found in differences in invasion speeds for different 0 C values. For very small a c values (a c < 1), the survival probabil-ity of such long distance dispersers increases with a c due to the increase in the denominator for the effective population size calculation. For larger 0 C values however, the probability of survival for such a long distance disperser decreases with oc because the effective population size felt by such an individual increases toward 0.5 as oc increases. In biological terms, long distance dispersers are more able to escape the effects of the central population if a c is small than if 0"C is large. It is important to note that this problem does not occur in the deterministic simula-tion because such long distance dispersers do not occur. Rather, very low densities are present at all spatial coordinates in the tails of the invasion wave and although these densities have increased growth when 0 C is small, this growth is limited by both the small population densities creating it and by the effects of the non-zero population densities at all coordinates surrounding them. Long-distance dispersers wil l be present in every individual-based stochastic simulation, and so one would expect to always see some difference in invasion speed at different values of o"c. However, we might expect that this effect could be decreased by moving closer to the deterministic assumptions of infinite population size. One method of increasing population size while maintaining the same carry-ing capacity is to decrease the mesh size. To examine whether changing mesh size would in fact diminish the differences in invasion speed observed when 0 C values were varied, I conducted the stochastic simulations with a quasi-local interaction for various mesh sizes. Figure 3.16 shows the effect of changing mesh size on the relationship between invasion speed and oc. As can be seen in this figure, although decreasing the mesh size increases in-vasion speeds for all values of 0 C , it did not noticeably alter the trend for 0 C > 1 of decreasing invasion speed as oc is increased. In retrospect, this could have been expected, since decreasing the mesh size both increases the number of individu-als in the tails and decreases the effect they have on their neighbours because of the mesh size multiple h in the quasi-local interaction kernel sum/approximation, resulting in essentially the same difference in survival rates for long distance dis-persers at different 0 C values. Decreasing the mesh size does have a significant effect on the trend of increasing invasion speed with increasing o c values for low 0 C . The differences in invasion speeds at low 0 C values were diminished consider-ably when the mesh size was decreased. Again, this response is predictable, since decreasing the mesh size decreases the effect a single individual has on itself by Chapter 3. Results 42 w c o -4—' co CD c CD CD o in k_ < •a CD CD Q. CO c o CO CO > 1.05 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 fc=0.1 /z=0.05 /?=0.025 /i=0.01 x x v > s < > x xx c x ^ x X x ++ + + + ++ + 10 Figure 3.16: Invasion Speed at Generation 50 for Various Mesh Sizes for the Stochastic Simulation With a Quasi-Local Interaction. Invasion speeds are calculated as the mean distance between furthest individ-uals in generations 25 to 50. Produced by the stochastic simulation of the integrodifference equation with Beverton-Holt local dynamics (a = 1, X = 2), a 7V(0,1) dispersal kernel, a mesh size of h and a N(0,oc) quasi-local interaction kernel. decreasing the multiple h in the effective population size calculation. When o c is small, the effective population size for an individual that is a long distance from the remainder of the population is almost solely determined by the individual's effect on itself. Decreasing the mesh size therefore results in a lower effective population size, and better survival, for long distance dispersers. With the evidence that the survival of long distance dispersers varies with o c , I must consider the possibility that this different trend in the asymptotic invasion speed for the stochastic simulations is an artifact of my technique used to estimate the location of the "front" of the invasion wave. Since in my stochastic simulations 1 used the furthest individual from the origin at each time step as an indicator of the location of the front of the invasion wave, it is possible that I have been Chapter 3. Results 43 (a) o c = 2.0 (b) o c = 8.0 -200 -150 -100 -50 0 50 100 150 200 -200 -150 -100 -50 0 50 100 150 200 .v (spatial coordinate) x (spatial coordinate) Figure 3.17: Average Population Density Waves at Generations 100 and 200 for the Stochastic Simulation With a Quasi-Local Interaction. Produced by the stochastic simulation with Beverton-Holt local dynamics (a = 1, A. = 2), ayV(0,1) dispersal kernel, a mesh size of 0.1 and &N(0,oc) quasi-local interaction kernel. Population densities shown are the average of 200 runs of the stochastic simulation (extinct populations are excluded from this average). merely tracking the progress of outliers across space. Since these outliers survive differently depending on Gc I may therefore have detected differences in invasion speed that may not actually be differences representative of the remainder of the invasion wave. For instance, perhaps the tails of the invasion waves are actually stretching over time as these outliers move faster than the remainder of the wave. However, from a simple visual inspection of the invasion waves, there does appear to be a difference in invasion speed that depends on o c irrespective of which part of the wave is tracked. For instance, in Figure 3.17 we can clearly see the difference in distance the invasion waves have travelled for the two different values of o c . It therefore seems that tracking the furthest individual did provide a reasonable approximation of the entire wave's movement. 2. Wavefront Shape In the single run simulations with a quasi-local interaction, the wavefront shape is quite variable and it is very difficult to find consistent differences between sim-ulations. Only for very large values of o"c were significant changes in the single run wavefront shape detected. However, as mentioned above, once they were av-eraged over a large number of iterations, the shape of the invasion waves for the Chapter 3. Results 44 stochastic simulations with a quasi-local interaction were quite similar to those produced by the deterministic simulations. As an example, Figure 3.18 illustrates the peaks of increased density at the front of the invasion waves, as well as the oscillations back toward the origin, that occur when the value of o"c is large. x ( s p a t i a l c o o r d i n a t e ) Figure 3.18: Average Population Density Wave at Generation 100 for the Stochas-tic Simulation with a Quasi-Local Interaction. Produced by the stochastic simulation of the integrodifference equation with Beverton-Holt local dynamics (a = 1, A. = 2), a iV(0,1) dispersal ker-nel and a7V(0,10) quasi-local interaction kernel. However, because of the inherent variation in the stochastic simulations, the exact height of the wavefront peaks is very difficult to estimate, although from visual inspection the peaks heights do appear to be approximately the same as those of the corresponding deterministic simulations. As well, the trend of peak height increase as the value of o c is increased is also maintained in the stochastic simulations. Finally, for smaller values of a c (less than 2.5) the peaks are virtu-ally impossible to distinguish visually from the other variations around carrying capacity, even when a large number of simulations have been averaged. It is there-fore too difficult in these stochastic simulations to estimate accurately the value of Chapter 3. Results 45 G c at which the peaks first appear, or to find the ratio between <5d and oc required for peaks to form. 3. Spatial Pattern As with the wavefront shape, the average population density at all spatial coordi-nates was similar to that of the deterministic simulations, with spatial coordinates sufficiently far behind the wavefront being quite close to the predicted carrying capacity. In addition, although single simulations showed large variations in pop-ulation sizes both at the wavefront edge and behind the wavefront where the pop-ulation was at "equilibrium", no distinguishable repeated patterns were formed in single simulations. Some single simulations did have long distance dispersal events that resulted in temporary local patches. However, although it seems likely that the additional spatial coupling created by a quasi-local interaction would af-fect the formation of local patches, I was not able to identify any trends in the frequency or size of these temporary local patches as G c was varied. Chapter 4. Discussion and Conclusions 46 Chapter 4 Discussion and Conclusions 4.1 Summary and Discussion The primary purpose of this study was to help establish whether adding a quasi-local interaction to integrodifference invasion models significantly alters the speed, wavefront shape, and spatial pattern formed by the invading population. For the particular system examined in this paper, an integrodifference model with a Gaus-sian dispersal kernel, Beverton-Holt local dynamics, and Gaussian quasi-local in-teraction kernel, it appears that adding a quasi-local interaction can alter all of these aspects in some way. Since adding a quasi-local interaction to the deterministic version of this par-ticular invasion model did not appear to alter the asymptotic speed of propagation of the travelling waves produced, this study does not show that it is necessary to include quasi-local interactions in deterministic invasion models in order to get accurate predictions of asymptotic invasion speeds. However, the addition of a quasi-local interaction to this model did cause significant changes in the initial rates of growth and spread for the corresponding invasion, which suggests that it may be beneficial to include quasi-local interactions in invasion models that are used to study initial behaviour. Furthermore, since the initial growth and spread of the invasion was found to depend on the standard deviation of the quasi-local interaction kernel in this system, adding quasi-local interactions to deterministic invasion models could allow them to express a larger range of initial behaviour that can be tuned by varying the parameters of the quasi-local interaction ker-nel, without affecting the asymptotic invasion speed. In this sense, quasi-local interactions could provide a valuable tool for models that focus on the speed of a biological invasion through time. In this study, the most dramatic difference caused by the addition of a quasi-local interaction to the integrodifference invasion model was in the shape of the travelling wave. The addition of a quasi-local interaction with a sufficiently large standard deviation relative to that of the dispersal kernel (oc > 1.2 • 0 ^ ) caused the wavefront shape to shift from a simple Gaussian wavefront to a damped os-Chapter 4. Discussion and Conclusions 47 dilatory wavefront. The new wavefront form consists of a peak of high density (above carrying capacity) at the front followed by spatially damped oscillations (oscillating around carrying capacity) behind the front. In addition, by varying the standard deviation of the quasi-local interaction the amplitude of these oscil-lations is varied. Since the addition of a quasi-local interaction to this invasion model greatly broadens the range of possible travelling wave shapes, it is proba-ble that they would also broaden the wave shapes in other invasion models as well. Therefore, quasi-local interactions should be considered as a possible component of invasion models that focus on wavefront shape. Finally, aside from the changes that arise from the different wavefront shapes, the addition of a quasi-local interaction to this invasion model was not found to alter the spatial pattern formed by the invading population. It is therefore still unknown whether the addition of a quasi-local interaction to integrodifference invasion models can increase the complexity of spatial patterns and whether is would be profitable to add quasi-local interactions to invasion models that focus on modelling the complex spatial patterns formed by biological invasions. A secondary goal of this study was to establish whether a more realistic stochas-tic and individual-based model of this system differed significantly from the deter-ministic model in the invasion speeds achieved, and the wavefront shapes and spa-tial patterns formed. The simulations of the analogous individual-based stochastic process revealed a number of interesting differences from the behaviour of the de-terministic system. With regards to the wavefront shape and the spatial pattern formed, there were certainly differences between the results of single runs and the deterministic simulations. However, no consistent significant differences were noted in the single simulations, and when averaged over a large number of sim-ulations the shape of the invasion wave and the carrying capacity achieved were nearly identical to those of the deterministic simulations for all stages of the inva-sion and for all values of G c and without a quasi-local interaction. The most dramatic difference in behaviour was seen when examining the inva-sion speeds of the stochastic simulations, and how the invasion speed was affected by the addition of a quasi-local interaction. First, the invasion speeds achieved by the stochastic simulations even without a quasi-local interaction were significantly lower than those achieved for the corresponding deterministic simulations, a dif-ference in stochastic invasion models that has long been conjectured and received much support from numerical studies [14, 20]. As well, as the population size was increased by the reduction of the mesh size, thereby approaching the deter-ministic model which assumes infinite population size, the invasion speeds of the stochastic simulations increased substantially. This result is akin to one found in Chapter 4. Discussion and Conclusions 48 the analysis of a stochastic nonlinear integrodifference invasion model by Lewis (2000a), who found similar increases in invasion speed by decreasing the radius of the local neighbourhoods in his model [14]. Finally, a result that was not ex-pected was that the invasion speeds achieved, both initial and asymptotic, for the stochastic simulations depended on the standard deviation oc of the quasi-local interaction kernel. In contrast, the asymptotic invasion speeds for the determinis-tic simulations were found to be the same for all values of o c as well as for the simulations without a quasi-local interaction. In addition, the dependence of the asymptotic invasion speed on a c in the stochastic simulations did not diminish as the mesh size was decreased, implying a fundamental difference in the function of the quasi-local interaction in the deterministic and stochastic models. This ex-ample serves to emphasize that stochastic models do not always behave as their deterministic analogues do, and although it may not be beneficial to add a quasi-local interaction to a deterministic model, it can produce interesting results in stochastic ones. 4.2 Conclusion and Recommendations for Future Work The addition of a quasi-local interaction to this simple integrodifference invasion model has been shown to produce interesting changes in some of the features of most interest to both theoretical and applied ecologists. Most significantly, this study has shown how the addition of a quasi-local interaction can produce more complex wavefront shapes in classical integrodifference invasion models, as well as alter both the spatial pattern and invasion speed of the invasion in the initial stages. With regards to stochastic models, this study has provided added evidence that stochastic and individual-based invasions proceed at significantly slower speeds than their deterministic analogues. As well, and perhaps more im-portantly, this study has demonstrated that the invasion speed achieved by such a stochastic model can be affected, both initially and asymptotically, by the addition of a quasi-local interaction. These results show that quasi-local interactions pro-vide a valuable means for adding spatial coupling to invasion models, and thereby to create more complex behaviour. Judging by the ability of this specific type of quasi-local interaction kernel to broaden the possible behaviours of this specific integrodifference model, this suggests that quasi-local interaction theory may have the potential to greatly enrich the range of behaviours available in other types of Chapter 4. Discussion and Conclusions 49 invasion models as well, and the field of biological invasion theory in general. Of course, this study has only considered the effects and importance of includ-ing a specific type of quasi-local interaction into an integrodifference model with a single functional form for the dispersal kernel and for the local growth dynamics. Future work should examine the effects of introducing quasi-local interactions to invasion models with other functional forms for the dispersal kernel and growth function, as well as different forms for the quasi-local interaction kernel. In addi-tion, as this study considered only a 1-dimensional biological invasion, it would be very interesting to examine how the results found in this study translate to the 2 and 3-dimensional cases. As well, analytical analysis of integrodifference equa-tions that include a quasi-local interaction should be undergone to determine the conditions under which such an invasion would have an asymptotically constant rate of spread, and how Weinberger's formula for the value of this constant would be changed by such an addition. More generally, introducing quasi-local interactions to invasion models with multi-species interactions, for instance with predator-prey dynamics, host-parasitoid relationships, or interspecific competition could be very rewarding. Since this study considered only the addition of quasi-local interaction via intraspecific com-petition and to a single-species invasion model, I cannot speculate how the dy-namics of multi-species models would change with the addition of a quasi-local interaction. However, as few biological invasions actually occur with no inter-species interactions, and quasi-local interactions can be expected to occur for a broad range of species, the examination of multi-species invasion models that in-clude a quasi-local interaction is a logical next step in this research. Bibliography 50 Bibliography [1] M . Andersen. Properties of some density-dependent integrodifference equa-tion population models. Mathematical Biosciences, 104:135-157, 1991. [2] R. J. H . Beverton and S. J. Holt. On the Dynamics of Exploited Fish Popu-lations, volume 19(2) of Fish Investigations. Ministry of Agriculture, Fish and Food, 1957. [3] J. Van den Bosch, A . J. Metz, and O. Diekmann. The velocity of spatial population expansion. J. of Mathematical Biology, 28:529-565, 1990. [4] M . Doebeli and T. Killingback. Metapopulation dynamics with quasi-local competition. Theoretical Population Biology, in press, 2003. [5] R. A . Fisher. The wave of advance of an advantageous gene. Annals of Eugenics, 7:355-369, 1937. [6] D. P. Hardin, P. Takac, and G. F. Webb. Dispersion population models dis-crete in time and continuous in space. J. of Mathematical Biology, 28:1-20, 1990. [7] R. Hengeveld. Dynamics of Biological Invasions. Chapman and Hall Ltd., London, 1989. [8] E. E. Holmes, M . A . Lewis, J. E. Banks, and R. R. Veit. Partial differential equations in ecology: Spatial interactions and population dynamics. Ecol-ogy, 75:17-29, 1994. [9] D. G. Kendall. Mathematical Models of the Spread of Infection, in: Mathe-matics and Computer Science in Biology and Medicine. Medical Research Council, 1965. [10] D. Kincaid and W. Cheney. Numerical Analysis. Brooks/Cole Publishing Co., Pacific Grove, second edition, 1996. Bibliography 51 [11] M . Kot. Discrete-time travelling waves: Ecological examples. J. of Mathe-matical Biology, 30:413^136, 1992. [12] M . Kot, M . A . Lewis, and P. van den Driessche. Dispersal data and the spread of invading organisms. Ecology, 77:2027-2042, 1996. [13] M . Kot and W. M . Schaffer. Discrete-time growth-dispersal models. Math-ematical Biosciences, 80:109-136, 1986. [14] M . A . Lewis. Spread rate for a nonlinear stochastic invasion. J. of Mathe-matical Biology, 41:430-454, 2000. [15] M . A . Lewis and P. Kareiva. Allee dynamics and the spread of invading organisms. Theoretical Population Biology, 43:141-158, 1993. [16] M . A . Lewis and S. Pacala. Modeling and analysis of stochastic invasion processes. J. of Mathematical Biology, 41:387^429, 2000. [17] R. Lui . Existence and stability of travelling wave solutions of a nonlinear integral operator. J. of Mathematical Biology, 16:199-220, 1983. [18] D. Mollison. Possible velocities for a simple epidemic. Advances in Applied Probability, 4:233-257, 1972. [19] D. Mollison. The rate of spatial propagation of simple epidemics. Proceed-ings of the Sixth Berkeley Symposium, 3:579-614, 1972. [20] D. Mollison. Spatial contact models for ecological and epidemic spread. J. of Royal Statistical Society, Series B, 39:283-326, 1977. [21] W. H . Press. Numerical Recipes in C++: The Art of Scientific Computing. Cambridge University Press, Cambridge, second edition, 2002. [22] A . Sasaki. Clumped distribution by neighbourhood competition. J. of Theo-retical Biology, 186:415-430,1997. [23] N . Shigesada and K . Kawasaki. Biological Invasions: Theory and Practice. Oxford Series in Ecology and Evolution. Oxford University Press, Oxford, 1996. [24] J. G. Skellam. Random dispersal in theoretical populations. Biometrika, 38:196-218, 1951. Bibliography 52 [25] C. Starr and R. Taggart. Biology: The Unity and Diversity of Life. Wadsworth Publishing Co., Toronto, seventh edition, 1995. [26] H . R Weinberger. Asymptotic Behavior of a Model in Population Genetics. Number 648 in Lecture Notes in Mathematics. Ed. J. M . Chadam, 1978. [27] H . F. Weinberger. Long-Time Behavior of a Class of Biological Models. Ed. W. E. Fitzgibbon. Pitman, London, 1984. [28] P. C. L . White and G. Newton-Cross. An Introduced Disease in an Invasive Host: The Ecology and Economics of Rabbit Calcivirus Disease (RCD) in Rabbits in Australia. Eds. C. Perrings, M . Williamson and S. Dalmazzone. Edward Elgar Publishing Inc., Cheltenham, 2000. Appendix A. Source Code for Deterministic Simulations 53 Appendix A Source Code for Deterministic Simulations /* This is the deterministic version of the invasion model including quasi-local interactions. Starting from a hard coded i n i t i a l population density distribution, the program iterates the integrodifference equation with a N(0,sigma-d) dispersal kernel and a N(0, sigma-c) quasi-local interaction kernel for a specified number of generations and on a fixed grid of spatial coordinates with a specified spacing "mesh-size" between grid points. The program does this for a specified range of values for sigma-c. The population density at each coordinate at the start of each generation is output to the f i l e "wave_graph_sigma_ disp_sigma-d_sigmacomp_sigma-c_generations_number-of-generations_mesh_mesh-size.dat" The maximum population density in each generation is output to the f i l e "density_limit_sigma_disp_sigma-d_sigmacomp_ sigma-c_generations_number-of- generations_mesh_mesh-size.dat" At each generation, the generation number and the rightmost coordinate to have exceeded the preset density threshold of 0.001 are output to the f i l e "maxdistance_determ.dat" For each value of sigma-c, the current sigma-c value and the calculated invasion speed achieved by the simulation (calculated as the average distance between successive wavefronts for the non-discarded points in the f i l e "maxdistance_determ.dat") are output to the f i l e "invasion_speed_sigma_disp_sigma-d_generations_number-of-Appendix A. Source Code for Deterministic Simulations 54 generations_mesh_mesh-size.dat" */ t i nc lude <iostream> t i nc lude <math.h> t i nc lude <fstream> t i nc lude <stdl ib .h> t i n c l u d e <string> using n a m e s p a c e s td ; v o i d i n i t i a l i z e _ d e n s i t y _ a r r a y ( ) ; //zeros the entries in the density array, and places the //appropriate values in the grid points where the i n i t i a l //population density is non-zero v o i d output_densi ty_array() ; //outputs the pop_density array to the screen. d o u b l e effect_of_neighbour ( i n t y_grid_no, d o u b l e x, d o u b l e y) ; //This function calculates the strength of competition //experienced at the point x from neighbours at the point //y. y_grid_no is the grid point number (and hence the //array index) corresponding to y d o u b l e ef fect ive_pop ( i n t x_grid__no, d o u b l e x ) ; //This function determines the e f f e c t i v e population density //experienced by an i n d i v i d u a l at a point x. x_grid_no is //the grid point number corresponding to x d o u b l e growth_function ( i n t gr id_no, d o u b l e x ) ; //This is the function that describes the growth stage of //the l i f e c y c l e . The Beverton-Holt stock-recruitment / / r e l a t i o n s h i p is assumed to govern the growth stage. d o u b l e d i s p e r s a l _ p r o b a b i l i t y ( d o u b l e x, d o u b l e y ) ; //This function determines the p r o b a b i l i t y of dispersing //from a point x to the point y d o u b l e grow_and_disperse ( i n t grid_no, d o u b l e x, d o u b l e y) ; //This is merely the product of the growth_function Appendix A. Source Code for Deterministic Simulations 55 //evaluated at y and the d i s p e r s a l _ p r o b a b i l i t y function //evaluated at (x,y). grid_no is the grid point number //corresponding to the point y d o u b l e i n tegra te ( d o u b l e (*func) ( i n t gr id_no, d o u b l e x, d o u b l e y ) , d o u b l e x) ; //This function integrates the function func over the / / e n t i r e mesh using the trapezoidal rule. Note this assumes //that x is fixed and we are i n t e g r a t i n g over y d o u b l e f ind_invasion__limit ( d o u b l e dens i ty_threshold) ; //finds the rightmost x-coordinate at which the //'density_threshold is exceeded d o u b l e average_wave_movement(); //using the file "maxdistance_determ.dat" as input, this //function returns the average distance the invasion wave //has t r a v e l l e d per generation. To eliminate i n i t i a l //behaviour, a number the data points at the beginning of //the file are discarded (the number is set by the discard / / v a r i a b l e ) . c o n s t d o u b l e mesh_size = 0 .1 ; //Spacing between grid points (ie. mesh coarseness) c o n s t d o u b l e max_expected_distance = 100.0; //Beyond this x-coordinate all population densities are //assumed to be zero and are ignored. For all simulations / / t h i s was set to at least twice the predicted invasion //speed m u l t i p l i e d by the number of generations. c o n s t i n t no_grid_points = int(max_expected_distance/mesh_size)+1; //the number of grid points to be tracked. For most runs //due to symmetry only the right-hand wave was tracked. c o n s t d o u b l e lambda = 2 .0 ; / / i n t r i n s i c rate of growth in the Beverton-Holt equation c o n s t d o u b l e a = 1.0; //scales carrying capacity in the Beverton-Holt equation c o n s t d o u b l e sigma_disp = 1.0; //the value of the parameter sigma_d in the gaussian //dispersal kernel c o n s t d o u b l e p i = 3.141592654; Appendix A. Source Code for Deterministic Simulations 56 c o n s t i n t generations =50; //the total number of generations for the simulation to be //run. c o n s t i n t d i s ca rd = 25; //The number of the first generations to discard when //calculating the invasion speed, (usually set to half of //"generations") c o n s t d o u b l e sigma_comp_min = 0 .1 ; //the minimum value to be used for the standard deviation //of the gaussian quasi-local interaction kernel c o n s t d o u b l e sigma_comp_max = 0 .1; //the maximum value for the standard deviation of the //quasi-local interaction kernel c o n s t d o u b l e sigma_comp_increment = 0 .1 ; //the spacing between consecutive values of the standard //deviation of the quasi-local interaction kernel d o u b l e sigma_comp; //will contain the current value of the quasi-local //interaction kernel standard deviation d o u b l e pop_dens i ty[no_gr id_poin ts ] [4] ; //The first column will contain all of the population //densities N(x) in the current time step. The second //column will contain the effective population Neff(x) at //the grid point. The third column will contain the growth //portion (f(N(x)) for each grid point. The fourth column //will contain the next year's density at the corresponding //grid point ofstream density_graph_stream; / / o u t p u t f i l e stream for the population density graphs ofstream invas ion_ l imi t_s t r eam; //output f i l e stream for the location of the invasion front //at each time step ofstream dens i ty_ l imi t_s t r eam; / / o u t p u t f i l e stream for the maximum density at each time //step ofstream invasion_speed_stream; //output fiel stream for the invasion speed at each //sigma_comp value Appendix A. Source Code for Deterministic Simulations 57 i n t main() { s t r i n g density_graph_stream_name; //the name for the population density graphs file. s t r i n g density_limit_stream_name; //the name for the maximum density file s t r i n g invasion_speed_stream_name; //the name for the invasion speed file s t r i n g i s s n l = "invasion_speed_sigma_disp_"; //first constant part of invasion speed file name s t r i n g i ssn2 = "_generations_"; //second constant part of invasion speed file name s t r i n g i ssn3 = "_mesh_"; //third constant part of invasion speed file name s t r i n g dat_ext = " .dat" ; //file name extension (same for all output files) sigma_comp = sigma_comp_min; //initialize the sigma_comp variable to the minimum //value c h a r buffer_sigma_disp[20] ; //temporary buffers for constructing file names c h a r buffer_generations[20] ; c h a r buffer_mesh_size[20]; c h a r buffer_sigma_comp[20] ; gcvt (sigma_disp, 6,buf f er__sigma_disp) ; //assigns the above char arrays the appropriate values gcvt (genera t ions ,6 ,buf fe r_genera t ions ) ; gcvt(mesh_size,6,buffer_mesh_size) ; invasion_speed_stream_name = i s s n l + buffer_sigma_disp + issn2 + buffer_generat ions + issn3 + buffer_mesh_size + dat_ext; invasion_speed_stream.open(invasion_speed_stream_name.c_str()); //creates and opens the file for invasion speed stream while(sigma_comp<=sigma_comp_max+0.01) Appendix A. Source Code for Deterministic Simulations 58 //outer loop. Iterates over all values of sigma_comp //from sigma_comp_min to sigma_comp_max by incrementing //sigma_comp by sigma_comp_increment each time { d o u b l e x_coord = 0.0; //variable to keep track of the x-coordinate //corresponding to a grid point d o u b l e i n v a s i o n _ l i m i t = 0.0; //the current location of the wave front d o u b l e invasion_speed = 0.0; //the invasion speed achieved by the current //sigma_comp value d o u b l e max_density = 0.0; //variable for the maximum density achieved in the //current generation gcvt(sigma_comp,6,buffer_sigma_comp); s t r i n g dgsnl = "wave_graph_sigma_disp_"; //first constant part of density graph file name s t r i n g dgsn2 = "_sigmacomp_"; //second constant part of density graph, maximum //density, and density limit file names s t r i n g d l s n l = "dens i ty_ l imi t_s igma_disp_" ; //first constant part of density limit file name density_graph_stream_name = dgsnl + buffer_sigma_disp + dgsn2 + buffer_sigma_comp + issn2 + buffer_generat ions + issn3 + buffer_mesh_size + dat_ext; density_graph_stream.open(density_graph_stream_name.c_str()); //creates and opens the file for the density graph //stream density_limit_stream_name = d l s n l + buffer_sigma_disp + dgsn2 + buffer_sigma_comp + issn2 + buffer_generat ions + issn3 + buffer_mesh_size + dat_ext; Appendix A. Source Code for Deterministic Simulations 59 densi ty_l imi t_s t ream.open(dens i ty_l imi t_s t ream_name.c_s t r ( ) ) ; //creates and opens the f i l e for the density limit //stream invasion_limit_strearn.open("maxdistance_determ.dat") ; //creates and opens the f i l e for the invasion limit //stream i n i t i a l i z e _ d e n s i t y _ a r r a y ( ) ; / / i n i t i a l i z e s the pop_density matrix output_densi ty_array() ; //sends the starting densities to the density graph //file f o r ( i n t i=0; i<generat ions; i++) / / t h i s loop i t e r a t e s for the number of generations //specified (so i=current generation) { max_density = 0.0; //resets max_density to zero each generation i n v a s i o n _ l i m i t = f a b s ( f i n d _ i n v a s i o n _ l i m i t ( 0 . 0 0 1 ) ) ; //determines the rightmost coordinate to exceed //the density threshold of 0.001 i nvas ion_ l imi t_s t r eam << i << " " << i n v a s i o n _ l i m i t << endl ; //send location of wavefront to regression f i l e f o r ( i n t j=0; j<no_grid_points; j++) //for each grid point (each row in the //pop_density array). so j= current grid point //This loop calculates the effective population //density at all grid points { x_coord = mesh_size*j; //finds the x-coordinate corresponding to / / t h i s grid p o i n t . Appendix A. Source Code for Deterministic Simulations 60 pop^densi ty t j ] [1] = e f f ec t ive_pop( j , x_coord); //calculates the effective population density //at this grid point } f o r ( i n t j=0; j<no_grid_points; j++) //for each grid point (j =current grid point) //This loop calculates the next generation's //population density at all grid points and //determines the maximum population density this //generation i x_coord = mesh_size*j; //finds the corresponding x-coordinate pop_densi ty[ j ] [3] = integrate(grow_and_disperse, x_coord); //calculates the new population density at x //and places the value in column 3 i f (pop_dens i ty [ j ] [0 ] > max_density) { max_density = pop_dens i ty [ j ] [0 ] ; } //if this grid point's population density is //greater than the maximum density at all //previous grid points, then the max_density //variable is updated to this pop. density } f o r ( i n t j=0; j<no_grid_points; j++) //for each grid point (j=current grid point) //This loop copies the next generation's //population densities in column 3 into column 0 { pop_densi ty[ j ] [0] = pop_dens i ty [ j ] [3 ] ; } output_dens i ty_array() ; //output to the density graph file Appendix A. Source Code for Deterministic Simulations 61 dens i ty_ l imi t_s t r eam << i << " " << max_density << endl ; //output to the maximum density file } //end generations loop invasion_speed = average_wave_movement(); //calculates the invasion speed for this sigma_comp //value invasion_speed_stream << sigma_comp << " " << invasion_speed << endl ; //output to the invasion speed file sigma_comp = sigma_comp + sigma_comp_increment; //increment sigma_comp density_graph_stream.close() ; //close files with sigma_comp specific file names i nvas ion_ l imi t_s t r eam.c lose ( ) ; dens i ty_ l imi t_s t r eam.c lose ( ) ; } //end sigma_comp loop invasion_speed_stream.close() ; //close the invasion speed file r e t u r n 0; } v o i d i n i t i a l i z e _ d e n s i t y _ a r r a y ( ) { f o r ( i n t i=0; i<no_grid_points ; //zeros all entries of the i++) pop_density array { pop_dens i ty[ i ] [0] = 0.0 pop_dens i ty[ i ] [1] = 0.0 pop_dens i ty[ i ] [2] = 0.0 pop_dens i ty[ i ] [3] = 0.0 Appendix A. Source Code for Deterministic Simulations 62 } //now assign any non-zero population densities desired //for the initial conditions pop_density[0][0] = 0 .1 ; pop_density[1][0] = 0 .1; } v o i d output_densi ty_array() { d o u b l e x_coord = 0.0; //variable for the x-coordinate corresponding to a grid //point f o r ( i n t i=0; i<no_grid_points ; i++) //for each grid point (i=current grid point) //outputs the correponding x-coordinate and population //density at the x-coordinate to the density graph file { x_coord = mesh_size*i ; //calculates the corresponding x-coordinate density_graph_stream << x_coord << " << pop_dens i ty[ i ] [0] << endl ; density_graph_stream << endl << endl ; //carriage returns for spacing in output file } d o u b l e ef fective_pop ( i n t x_grid_no, d o u b l e x) { d o u b l e pop_effect = 0.0; pop_effect = in tegrate(effect_of_neighbour , x) ; //integrates the effect_of_neighbour function over y //for this fixed x. r e t u r n pop_effect; } Appendix A. Source Code for Deterministic Simulations 63 d o u b l e effect_of_neighbour ( i n t y_grid_no, d o u b l e x, d o u b l e y) { d o u b l e e f fec t = 0.0; i f ( y _ g r i d _ n o >=0) //because only one side of the x-axis is explicitly //tracked in the pop_density array, two cases must be //considered. First, if the y-value is non-negative, //then the effect at x of neighbours at y can be //calculated using the population density at y { ef fec t = (1/(sigma_comp*sqrt(2*pi)))*exp( -(pow((x-y),2.0))/(2*sigma_comp*sigma_comp)) *pop_density[y_grid_no][0]; //this is the population density at y weighted by a //gaussian probability density } e l s e //otherwise- the y-value is negative and not explicitly //tracked in the population density array, so the //population density at -y is used in the calculation { ef fec t = (1/(sigma_comp*sqrt(2*pi)))*exp( - (pow((x-y),2.0))/(2*sigma_comp*sigma_comp)) *pop_density[-y_grid_no][0]; } r e t u r n e f f e c t ; } d o u b l e growth_function ( i n t x_grid_no, d o u b l e x_coord) { d o u b l e f_N_x = 0.0; i f (x_grid_no>=0) //again, must consider 2 cases. First if x is //non-negative then we can use the value in the //population density array to calculate growth at x { Appendix A. Source Code for Deterministic Simulations 64 f_N_x = lambda*pop_density[x_grid_no][0]/ (l+a*pop_density[x_grid_no][1]) ; / / t h i s i s j u s t the Beverton Holt growth equation } e l s e //otherwise x i s negative so we use the population //density at -x in the calculation { f_N_x = lambda*pop_density[-x_grid_no][0]/ ( l+a*pop_densi ty[-x_grid_no][1]) ; } r e t u r n f_N_x; } d o u b l e d i spersa l_prob ( d o u b l e x, d o u b l e y) { d o u b l e p r o b a b i l i t y = 0.0; p r o b a b i l i t y = (1/(sqr t (pi*2)*sigma_disp))*exp( -pow((x-y) ,2 .0) / (2*s igma_disp*s igma_disp)) ; / / t h i s i s j u s t the value of the N(0,sigma_disp) //probability density function evaluated at (x-y) r e t u r n p r o b a b i l i t y ; } d o u b l e grow_and_disperse( int y_grid_no, d o u b l e x, d o u b l e y) { d o u b l e product = growth_function(y_grid_no, y )* d i spersa l_prob(x , y ) ; //simply the product of the growth function evaluated //at y and the probability of dispersal from y to x . / / t h i s i s the function that must be integrated each //generation r e t u r n product; } d o u b l e i n tegra te ( d o u b l e (*func) ( i n t y_grid_no, d o u b l e x, Appendix A. Source Code for Deterministic Simulations 65 d o u b l e y ) , d o u b l e x) { d o u b l e i n t e g r a l = 0.0; d o u b l e y = -no_grid_points*mesh_size; / / s t a r t y at the leftmost grid point f o r ( i n t k=-no_grid_points + 1; k<no_gr id_poin ts - l ; k++) / / t h i s loop runs over the entire mesh except the //endpoints. This i s the composite trapezoidal rule. { y = y + mesh_size; //moves over one grid point each iteration i n t e g r a l = i n t e g r a l + mesh_s ize*func(k ,x ,y) ; } i n t e g r a l = i n t e g r a l +l /2*mesh_size*(func(no_grid_points-1,x,y+mesh_size) + func(0,x,mesh_size* ( -no_gr id_poin ts ) ) ) ; //for composite trapezoidal rule need to add only 1/2 //of the first and last grid points r e t u r n i n t e g r a l ; } d o u b l e f i n d _ i n v a s i o n _ l i m i t ( d o u b l e densi ty_threshold) { f o r ( i n t k=no_gr id_poin ts - l ; k >=0; k—) //runs over all grid points, starting at the right end { i f ( p o p _ d e n s i t y [ k ] [ 0 ] > densi ty_threshold) r e t u r n mesh_size* (k) ; //if the density threshold i s exceeded, the /'/x-coordinate for this grid point i s returned //and the function i s exited } r e t u r n 0 .0; //if no coordinate exceeds the threshold, this returns //an x-coordinate of 0 Appendix A. Source Code for Deterministic Simulations 66 } d o u b l e average_wave_movement() { i n t number_of_points = 1000; //the maximum number of data points expected //must be greater or equal to the number of generations d o u b l e xcoords[number_of_points]; //the array of the x-values d o u b l e ycoords[number_of_points]; //the array of the corresponding y-values i f s t ream in_stream; //the input file stream f o r ( i n t j=0; j<number_of_points; j++) //runs over the arrays xcoords and ycoords and //initializes them { xcoords[ j ] = 0.0; ycoords[ j ] = 0.0; } in_stream.open("maxdistance_determ.dat"); //assigns the input stream to and opens the input file i n t coun te r l = 0; //tracks the number of data points that have been read //in from the input file w h i l e (! in_stream.eof()) //this loop places the input data from the file into //the xcoord and ycoord arrays { in_stream >> xcoords [coun te r l ] ; //xcoords[counterl] = next; in_stream >> ycoords [coun te r l ] ; //ycoords[counterl] = next; counterl++; //note that this increments counterl one too many Appendix A. Source Code for Deterministic Simulations 67 //times, which is why I use counterl-1 below, //note also that there are two blank lines at the //end of the input file, and so this loop proceeds //for one iteration too many } d o u b l e sum_y = 0.0; f o r ( i n t j=d iscard ; j<counter l - l ; j++) //This loop runs over the data values for generations //that are not desired to be discarded and sums the //differences in the location of the invasion front { sum_y = sum_y + ycoords[ j ] - y c o o r d s [ j - 1 ] ; } in_s t ream.close() ; r e t u r n ( sum_y/double(counter l - l -d iscard) ) ; //returns the average distance between successive //wavefront locations for the non-discarded generations Appendix B. Source Code for Stochastic Simulations 68 Appendix B Source Code for Stochastic Simulations /* This is the stochastic version of the integrodifference invasion model including quasi-local interactions. Starting from a hard coded initial population density-distribution, the program follows the following sequence for each generation. 1) the effective population density at each point is calculated (according to the quasi-local N(0,sigma-c) interaction kernel) 2) each individual is determined to either survive or die (based on the density dependent component of the Beverton-Holt stock recruitment equation) 3) each survivor produces offspring (Poisson distributed random deviates with a mean of lambda) and 4) offspring move a distance drawn from a N(0,1) distribution. The population exists on a fixed grid of spatial coordinates with a specified spacing "mesh-size" between adjacent grid points. The process is repeated for the desired number of "iterations" and for a specified range of values for sigma-c The average number of individuals at each coordinate at the start of each generation is output to the file "density_wave_sigma_comp_sigma-c_gens_number-of-generations_mesh_mesh-size.dat" For each generation, the generation number and the average location of the furthest individual from the origin are output to the file "maxdistance.dat" For each generation, the generation number and the average cumulative population size (the sum of all surviving individuals) is output to the file "popsize_sigma_comp_ Appendix B. Source Code for Stochastic Simulations 69 sigma-c_gens_number-of-generation_mesh_mesh-size.dat" For each value of sigma-c, the current sigma-c value and the calculated average invasion speed achieved by the simulations for that sigma-c value (calculated as the mean of the differences in the average location of the farthest individual in successive generations for the non-discarded points in the file "maxdistance_determ.dat") are output to the file "invasion_speed_avg_gens_number-of-generations_mesh_mesh-size.dat" The random number generator used and the routines for drawing Poisson deviates were provided to me by Rik Blok (2001). The source code for these routines is in the file "deviate.c" and was compiled into the executable. This source code will not be provided here but is available upon request (email: merchant@math.ubc.ca) */ t i n c l u d e <iostream> t i n c l u d e <stdl ib .h> t i n c l u d e <time.h> t i n c l u d e <math.h> t i n c l u d e <fstream> t i n c l u d e "devia te .h" //header for random number generator and Poisson //deviate routines f i nc lude <string> us ing n a m e s p a c e s td ; c o n s t d o u b l e mesh_size = 0 .1 ; //Spacing between grid points (ie. mesh coarseness) c o n s t d o u b l e max_dist = 275.0; //Beyond this x-coordinate all population densities are //assumed to be zero and are ignored. For all simulations //this was set to at least 1.25 times the predicted //invasion speed multiplied by the number of generations. c o n s t i n t number_of__bins = i n t (2*max_dist/mesh_size) + 1; //the number of grid points to be tracked. Appendix B. Source Code for Stochastic Simulations 70 c o n s t d o u b l e a=1.0; //scales carrying capacity in the Beverton-Holt equation c o n s t d o u b l e lambda = 2 .0 ; //intrinsic rate of growth in the Beverton-Holt equation c o n s t d o u b l e sigma_disp = 1.0; //the value of the parameter sigma_d in the gaussian //dispersal kernel (set to 1.0 for all simulations) c o n s t d o u b l e p i = 3.141592654; c o n s t i n t d i s c a r d = 100; //The number of the first generations to discard when //calculating the invasion speed, (usually set to half of //"generations") c o n s t i n t number_iterations=200; //the number of times simulations are run for each set of //parameter values (used to examine the mean behaviour) c o n s t i n t maximum_generations = 200; //the total number of generations for each simulation. c o n s t d o u b l e sigma_comp_min=0.1; //the minimum value to be used for the standard deviation //of the gaussian quasi-local interaction kernel c o n s t d o u b l e sigma_comp_max=0 . 5; //the maximum value for the standard deviation of the //quasi-local interaction kernel c o n s t d o u b l e sigma_comp_increment = 0 .1 ; //the spacing between consecutive values of the standard //deviation of the quasi-local interaction kernel ofstream maxdistance_stream; //output stream for the farthest distance per generation ofstream invasion_speed_avg_stream; //output stream for the average invasion speed ofstream density_wave_stream; //output stream for the average density graphs ofstream popsize_stream; //output stream for the population size per generation ofstream no_fa i l s_s t ream; //output stream for the number of extinctions per sigma-c //value v o i d i n i t i a l i z e _ m a t r i x ( d o u b l e matrix[][number_of_bins] , i n t rows); Appendix B. Source Code for Stochastic Simulations 71 //initializes the 2-D array matrix that has number-of-bins //columns and "rows" number of rows (0.0 is placed in all //entries) v o i d i n i t i a l i z e _ m a t r i x ( d o u b l e matrix[number_of_bins]) ; //initializes the 1-D array matrix that has number-of-bins //rows (0.0 is placed in all entries) v o i d c a l c u l a t e _ x _ c o o r d s ( d o u b l e matr ix[][number_of_bins]) ; //determines the corresponding x-coordinate for each grid //point/bin number and places it in the first row of the //array "matrix" v o i d o u t p u t _ l o c a t i o n _ s c r e e n ( d o u b l e matr ix[][number_of_bins] , i n t rows); //outputs the entire array "matrix" with "number-of-bins" //columns and "rows" number of rows to the screen v o i d i n i t i a l _ f i l l ( d o u b l e matr ix[][number_of_bins]) ; //places the desired initial local population sizes in the //corresponding entries of the array "matrix" (hard-coded) v o i d c a l cu la t e_Nef f ( d o u b l e matr ix [] [number_of_bins], d o u b l e sigma_value); //calculates the effective population size at each grid //point/bin number v o i d d e t e r m i n e _ s u r v i v o r s ( d o u b l e matr ix[][number_of_bins]) ; //determines the number of survivors at each grid point/bin //number v o i d c a l c u l a t e _ o f f s p r i n g ( d o u b l e matrix[][number_of_bins])• //determines the total number of offspring produced at each //grid point/bin number v o i d move_off sp r ing ( d o u b l e m a t r i x l [ ] [number_of_bins ] , d o u b l e matrix2[number_of_bins]); //moves the number of offspring at each coordinate in //matrixl, they are placed at their new location/grid point //in matrix2 Appendix B. Source Code for Stochastic Simulations 72 v o i d copy_matrix ( d o u b l e m a t r i x l [ ] [number_of_bins ] , d o u b l e matrix2[number_of_bins], i n t row_to_copy); //copies the "row-to-copy" of matrix2 over the same row of //matrixl (they must have "number-of-bins" columns) d o u b l e GaussDev ( d o u b l e mean, d o u b l e s t d ) ; //returns a Normal(mean, std) deviate d o u b l e maxdistance ( d o u b l e matr ix[][number_of_bins]) ; //finds the distance from the origin of the furthest //individual in the array "matrix" d o u b l e average_wave_movement(); //calculates the average distance between the y-values of //the file "maxdistance.dat". i n t main() { s t r i n g density_wave_file_name; //these strings will be used for the file names s t r i n g popsize_fi le_name; s t r i n g invasion_speed_file_name; s t r i n g number_of_fails_fi le_name; s t r i n g density_name_part_one = "density_wave_sigma_comp_"; //first constant part of density wave file name s t r i n g popsize_name_part_one = "popsize_sigma_comp_"; //first constant part of pop size file name s t r i n g invasion_name_part_one = "invasion_speed_avg"; //first constant part of invasion speed file name s t r i n g fails_name_part_one = "number_of_fails"; //first constant part of number of fails file name s t r i n g name_part_two = "_gens_"; //second constant part of file names s t r i n g name_part_three = "_mesh_"; //third constant part of file names s t r i n g fi le_name_extension = " .dat" ; Appendix B. Source Code for Stochastic Simulations 73 //file name extension (same for all output files) c h a r buffer_sigma_comp[20] ; //temporary holders of variables as char arrays for //creating file names c h a r buffer_gens[20]; c h a r buffer_mesh[20] ; gcvt(maximum_generations,6, buffer_gens); //place the appropriate values in the char arrays gcvt (mesh_size, 6,buf f er_m'esh) ; invasion_speed_file_name = invasion_name_part_one + name_part_two + buffer_gens + name_part_three + buffer_mesh + f i le_name_extension; invasion_speed_avg_stream.open(invasion_speed_fi le_name.c_str()) ; //creates and opens the invasion speed file number_of_fails_file_name = fails_name_part_one + name_part_two + buffer_gens + name_part_three + buffer_mesh + f i le_name_extension; no_fa i l s_s t ream.open(number_of_fa i l s_f i le_name.c_s t r ( ) ) ; //creates and opens the number of extinctions file d o u b l e sigma_comp = sigma_comp_min; //sets sigma-comp to the minimum value desired d o u b l e maxdist_matrix[maximum_generations]; //This array will contain the average distance //travelled by the wave in each generation d o u b l e density_wave_matrix[maximum_generations][number_of_bins]; //This array will contain the average density in each //bin number at each generation Note: this array is //extremely large. For simulations with a large number //of generations, I actually only keep track of every //fifth or tenth generation d o u b l e locat ion[6][number_of_bins] ; Appendix B. Source Code for Stochastic Simulations 74 //the array contains the following information //row 0 = x-value corresponding to each grid point //row 1 = number of individuals in the bin //row 2 = effective population density at x //row 3 = survival probability of individuals at x //row 4 = number of survivors at x //row 5 = number of total offspring produced at x d o u b l e temploc[number_of_bins]; //Temporarily holds the new population sizes at each //grid point after offspring have moved SetDevSeed(- l ) ; //seed the random number generator while(sigma_comp<=sigma_comp_max+0.01) //this outer loop runs over all the values of //sigma-comp desired. the variable sigma-comp is //incremented at the end each time { c h a r buffer_sigma_comp[20]; //temporary char array for the sigma-comp variable //for making file names gcvt(sigma_comp,6,buffer_sigma_comp); density_wave_file_name = density_name_part_one + buffer_sigma_comp + name_part_two + buffer_gens + name_part_three + buffer_mesh + f i le_name_extension; popsize_file_name = popsize_name_part_one + buffer_sigma_comp + name_part_two + buffer_gens + name_part_three + buffer_mesh + f i le_name_extension; / / c r e a t e s the d e n s i t y wave and population size file //names for this sigma-comp value d o u b l e population[maximum_generations]; Appendix B. Source Code for Stochastic Simulations 75 //this array will contain the average total //population size at each generation i n t number_of_fails = 0; //variable for the number of extinctions f o r ( i n t n=0;n<maximum_generations;n++) //initializes the population array populat ion[n]=0.0; i n i t i a l i z e _ m a t r i x ( d e n s i t y _ w a v e _ m a t r i x , maximum_generations); //initializes the density-wave-matrix array f o r ( i n t n=0; n<maximum_generations; n++) //initializes the average distance array maxdist_matrix[n] = 0.0; density_wave_stream.open(density_wave_fi le_name.c_str()) ; popsize_s t ream.open(popsize_f i le_name.c_s t r ( ) ) ; //opens the density wave and population size files f o r ( i n t n=0; n<number_iterations; n++) //this loop runs the stochastic simulation the //desired number of iterations (because it is //stochastic a large number of runs are conducted //and averaged { d o u b l e inv_speed_this_run = 0.0; //variable to track how quickly each individual //run expands maxdistance_stream.open("maxdistance.dat") ; //opens the file for the location of the //farthest individual i n i t i a l i z e _ m a t r i x ( l o c a t i o n , 6); //zeros the location matrix ca lcu la te_x_coords ( loca t ion) ; //determines the corresponding x-coordinate for //each grid point in the location array i n i t i a l _ f i l l ( l o c a t i o n ) ; //inputs the hard-coded initial population sizes Appendix B. Source Code for Stochastic Simulations 76 f o r ( i n t k=0;k<maximum_generations;k++) //generations loop. k = current generation { f o r ( i n t p=0; p< number_of_bins; p++) //this loop tallies up the total population //for this generation (note, the population //array is not zeroed each iteration, and //is divided by the number of non-extinct //populations at the end of the sigma-comp //loop for each sigma-comp value to give the //average population size in each generation { popula t ion[k] = popula t ion[k] + l o c a t i o n [ 1 ] [ p ] ; } c a l c u l a t e _ N e f f ( l o c a t i o n , sigma_comp); //calculates the effective population //density at all grid points determine_surv ivors ( loca t ion) ; //calculates the number of survivors (ie. //applies mortality) for all grid points c a l c u l a t e _ o f f s p r i n g ( l o c a t i o n ) ; //calculates the number of offspring //produced at each grid point i n i t i a l i z e _ m a t r i x ( t e m p l o c ) ; //zeroes the temploc array move_offspr ing( loca t ion , temploc); //determines the new spatial coordinates of //all offspring and places each at the //appropriate grid point in the temploc //array maxdist_matrix[k] = maxdis tance( locat ion) + maxdis t_mat r ix[k] ; //determines the location of the farthest //individual from the origin and adds this Appendix B. Source Code for Stochastic Simulations 77 //value to the maxdist array value for this //generation (note, the maxdist array is //averaged over non-extinct populations //after all iterations are complete for this //sigma-comp value) maxdistance_stream << k << " " << maxdis tance( loca t ion) ; //sends the generation number and the //distance from the origin of the farthest //individual to the "maxdistance.dat" file copy_mat r ix ( loca t ion , temploc, 1) ; //overwrites the individuals in the //location array with the offspring in the //temploc array (parents all die after //reproduction) f o r ( i n t p=0; p<number_of_bins; p++) //adds the current population size at each //grid point to the density wave matrix { density_wave_matrix[k][p] = density_wave_matrix[k][p] + l o c a t i o n [ 1 ] [ p ] ; ' } } //end of generations loop inv_speed_this_run = average_wave_movement(); //determines the average wave speed this //iteration i f ( i n v _ s p e e d _ t h i s _ r u n < 0.01) //if there is virtually no movement in the //non-discarded generations the population is //assumed to have gone extinct numbe r_o f_fa i1s++; maxdistance_stream.close (); //closes the maxdistance.dat file } Appendix B. Source Code for Stochastic Simulations 78 maxdistance_stream.open("maxdistance.dat"); //reopens the maxdistance.dat file (thereby //clearing it) f o r ( i n t n=0; n<maximum_generations; n++) //for each generation (n=current generation) this //loop averages the location of the farthest //individual, the population size, and the //population density at each grid point { maxdist_matrix[n] = maxdis t_matr ix[n] / double (number_ i t e ra t ions -number_of_ fa i l s ) ; maxdistance_stream << n << " " << maxdist_matrix[n] << end l ; //averages the location of farthest individual //and outputs the generation and location to //the maxdistance.dat file populat ion[n] = popu l a t i on [n ] / double (number_ i t e ra t ions -number_of_ fa i l s ) ; popsize_stream << n << " " << populat ion[n] << endl ; //averages the population size in each //generation and outputs it to the average //population size file f o r ( i n t p=0; p<number_of_bins; p++) //calculates the average population density in //this generation at each grid point and outputs //to the density wave file { density_wave_matrix[n][p] = densi ty_wave_matr ix[n][p] / (number_i terat ions-number_of_fai ls) ; density_wave_stream << loca t ion [0 ] [p ] << " " << density_wave_matrix[n][p] << endl ; Appendix B. Source Code for Stochastic Simulations 79 } density_wave_stream << endl << end l ; //just for formatting in density wave file invasion_speed_avg_stream << sigma_comp << " " << average_wave_movement() << end l ; //outputs the average distance between the average //location of farthest individuals in consecutive //generations to the invasion speed file no_fa i l s_s t ream << sigma_comp << " " << number_of_fails << end l ; //outputs the total number of extinctions for this //sigma-comp value sigma_comp=sigma_comp+sigma_comp_increment; //increments sigma-comp by the desired increment maxdistance_stream.close() ; //closes the files that are sigma-comp specific density_wave_stream.close() ; pops ize_s t ream.c lose ( ) ; } invasion_speed_avg_stream.close() ; no_fa i l s_s t ream.c lose ( ) ; //closes non-sigma-comp specific output files r e t u r n 0; } v o i d i n i t i a l i z e _ m a t r i x ( d o u b l e matr ix [] [number_of_bins], i n t rows) { f o r ( i n t i=0; i<rows; i++) //runs over all entries in the array "matrix" and sets //their values equal to zero { Appendix B. Source Code for Stochastic Simulations 80 f o r ( i n t j=0;j<number_of_bins;j++) { matrix [ i ] [j] = 0.0; } } } v o i d i n i t i a l i z e _ m a t r i x ( d o u b l e matrix[number_of_bins]) { f o r ( i n t i=0; i<number_of_bins; i++) //runs over all entries in the array "matrix" and sets //their values equal to zero m a t r i x [ i ] = 0.0; } v o i d c a l c u l a t e _ x _ c o o r d s ( d o u b l e matrix[][number_of_bins]) { f o r ( i n t i=0;i<number_of_bins;i++) //determines the x-coordinate corresponding to each //grid point in the array "matrix" by calculating how //far it must be from the left end of the grid { m a t r i x [ 0 ] [ i ] = -max_dist + mesh_size*i ; } } v o i d ou tpu t_ loca t i on_sc r een ( d o u b l e matr ix[][number_of_bins] , i n t rows) { f o r ( i n t i=0; i<number_of_bins; i++) //runs over all entries in the array "matrix" and //outputs them to the screen. output is in standard //math matrix format { f o r ( i n t j=0; j<rows; j++) { cout << m a t r i x [ j ] [ i ] << " "; } cout << endl ; } Appendix B. Source Code for Stochastic Simulations 81 } v o i d i n i t i a l _ f i l l ( d o u b l e matr ix [] [number_of_bins]) { matrix[1][number_of_bins/2]=1; //places one individual at the grid point corresponding / / t o x=0 } v o i d ca lcu la te_Nef f ( d o u b l e matr ix [] [number_of_bins], d o u b l e sigma) { d o u b l e e f f e c t ; f o r ( i n t i=0; i<number_of_biris; i + + ) / / t h i s loop calculates the effective population density //at each grid point { i f ( m a t r i x [ 1 ] [ i ] > 0.0) //if there i s no one alive at t h i s grid point the //effective population density i s set to 0.0 { effect = 0.0; f o r ( i n t j=0; j<number_of_bins; j++) / / t h i s loop sums up the weighted effect //(weighted by the quasi-local interaction //kernel) of all individuals on the grid //essentially t h i s i s performing a riemann sum //approximation of the integral form { effect=effeet + mat r ix [1] [ j ]*exp( - (pow(mat r ix [0 ] [ j ] -ma t r i x [0 ] [ i ] , 2 . 0 ) / (2*s igma*s igma) ) ) ; } ef fec t = mesh_s ize*(ef fec t ) / ( sqr t (2*pi )*s igma) ; //multiplies by the width of each bin to //approximate riemann sum and normalizes the //effective density as in the integral form m a t r i x [ 2 ] [ i ] = e f fec t ; Appendix B. Source Code for Stochastic Simulations 82 //places the value of the effective population //density in the third row of the array "matrix" } } } v o i d copy_matrix ( d o u b l e m a t r i x l [] [number_of_bins], d o u b l e matrix2[number_of_bins], i n t row_to_copy) { f o r ( i n t i=0; i<number_of_bins; i++) //simply copies row-to-copy of matrix 2 into the same //row of matrix 1 { matr ix l [ row_to_copy] [ i ] = m a t r i x 2 [ i ] ; } } v o i d d e t e r m i n e _ s u r v i v o r s ( d o u b l e matrix[][number_of_bins]) { f o r ( i n t i=0; i<number_of_bins; i++) / / t h i s loop determines the number of survivors at each //grid point { m a t r i x [ 4 ] [ i ] = 0; i f ( m a t r i x [ 1 ] [ i ] > 0.0) //if there are individuals at t h i s grid point, //calculate their survival probability. otherwise //set i t equal to 0.0 { matrix [3] [ i ] = 1 / ( l+a*mat r ix [2 ] [ i ] ) ; f o r ( i n t j=0; j < m a t r i x [ l ] [ i ] ; j++) { i f ( m a t r i x [ 3 ] [ i ] > UniformDev()) matr ix[4][ i ]++; } } } } Appendix B. Source Code for Stochastic Simulations 83 v o i d c a l c u l a t e _ o f f s p r i n g ( d o u b l e matrix[][number_of_bins]) { f o r ( i n t i=0; i<number_of_bins; i++) / / t h i s loop determines the number of offspring to be //produced at each grid point { m a t r i x [ 5 ] [ i ] = 0; i f ( m a t r i x [ l ] [ i ] > 0.0) //if there are adults at t h i s grid point, calculate //the number of offspring. Otherwise the number of //offspring i s set to 0 { f o r ( i n t j=0; j < m a t r i x [ 4 ] [ i ] ; j++) / / t h i s loop draws a Poisson(lambda) number of //offspring for each adult at t h i s grid point //and places the total number of offspring in //the sixth row { m a t r i x [ 5 ] [ i ] = m a t r i x [ 5 ] [ i ] + PoissonDev(lambda); //note: PoissonDev(lambda) draws a random //number from a Poisson distribution with //mean lambda } } } } v o i d move_off sp r ing ( d o u b l e m a t r i x l [] [number_of_bins], d o u b l e matrix2[number_of_bins]) { i n t d i r e c t i o n = 0; //if direction = 1 movement i s to the right, if //direction = -1 movement i s to the left d o u b l e d is tance = 0.0; i n t change__in_bin_number = 0; / / t h i s will be the number of grid points the individual //will be shifted from i t s current location Appendix B. Source Code for Stochastic Simulations 84 f o r ( i n t i=0; i<number_of_bins;i++) //this loop moves the offspring for each grid point { f o r f i n t j=0; j < m a t r i x l [ 5 ] [ i ] ; j++) //this loop is iterated once for each offspring at //this grid point { //first the direction of movement is determined //by pulling a uniform random (0,1) deviate i f (UniformDevO < 0.5) //if the deviate is less than 0.5 the individual //moves left d i r e c t i o n = - 1 ; e l s e //otherwise it moves right d i r e c t i o n = 1; dis tance = GaussDev(0, sigma_disp); //this determines the distance the individual //will move in that direction by pulling a //random deviate from a Normal(0, sigma-d) //distribution change_in_bin_number = d i r e c t i o n * i n t ( f l o o r ( d i s t a n c e / m e s h _ s i z e + 0 .5 ) ) ; //the number of grid points the individual will //be moved left or right is calculated here if(i+change_in_bin_number < number_of_bins && i+change_in_bin_number >= 0) matrix2[i+change_in_bin_number]++; //provided the new grid point is not outside of //the fixed grid established, the number of //individuals at the new grid point is //incremented by one e l s e cout « "OUT OF BOUNDS, LOST ONE" « endl ; //if the new grid point is out of range, an //error message is output to the screen } Appendix B. Source Code for Stochastic Simulations 85 } } d o u b l e GaussDev ( d o u b l e mean, d o u b l e std) { d o u b l e W = 1.0; d o u b l e wl = 0.0; d o u b l e w2 = 0.0; w h i l e (W>=1) //this is the standard method for drawing N(mean,std) //random deviates. See Press (2002) for details { wl = 2*UniformDev()- l ; w2 = 2*UniformDev()-1; W = wl*wl+w2*w2; } r e t u r n (mean + sqr t ( -2*log(W)/W)*wl*std) ; d o u b l e maxdistance ( d o u b l e matr ix [] [number_of_bins ]) { d o u b l e f a r t h e s t _ l e f t = 0.0; //the location of the farthest individual to the left d o u b l e f a r thes t_ r igh t = 0.0; //the loc. of the farthest indiv. to the right i n t coun te r l = 0; //counter that starts at left end of the grid/mesh i n t counter2 = number_of_bins - 1 ; //counter that starts at right end of the grid/mesh w h i l e ( f a r t h e s t _ l e f t = = 0 . 0 && coun te r l < number_of_bins) //this loop runs from the left end of the mesh and //checks each grid point until an individual is found { i f ( m a t r i x [ 1 ] [ c o u n t e r l ] != 0) f a r t h e s t _ l e f t = f a b s ( m a t r i x [ 0 ] [ c o u n t e r l ] ) ; counterl++; } Appendix B. Source Code for Stochastic Simulations 86 w h i l e (farthest_right==0.0 && counter2 >= 0) / / t h i s loop runs from the right end of the mesh until //an individual i s found { i f ( m a t r i x [ 1 ] [ c o u n t e r 2 ] != 0) f a r thes t_ r igh t = fabs(matr ix[0] [counter2]) ; counter2—; } i f ( f a r thes t_ le f t >= fa r thes t_ r igh t ) //if the leftmost individual i s farthest from the //origin return i t s distance r e t u r n f a r t h e s t _ l e f t ; r e t u r n f a r t he s t_ r i gh t ; //otherwise return the distance of the farthest right //individual (Note, if no one i s alive, 0 i s returned d o u b l e average_wave_movement() { i n t number_of_points = 1000; //the maximum number of data points expected //must be greater or equal to the number of generations d o u b l e xcoords[number_of_points]; //the array of the x-values d o u b l e ycoords[number_of_points]; //the array of the corresponding y-values i f s t r eam in_stream; //the input f i l e stream f o r ( i n t j=0; j<number_of_points; j++) //runs over the arrays xcoords and ycoords and //zeroes them { xcoords[ j ] = 0.0; ycoords[ j ] = 0.0; } Appendix B. Source Code for Stochastic Simulations 87 in_stream.open("maxdistance.dat") ; //assigns the input stream to and opens the input file i n t counter l = 0; //tracks the number of data points that have been read //in from the input file w h i l e (! in_stream. eof ()) //this loop places the input data from the file into //the xcoord and ycoord arrays { in_stream >> xcoords [coun te r l ] ; //xcoords[counterl] = next; in_stream >> ycoords [coun te r l ] ; //ycoords[counterl] = next; counterl++; //note that this increments counterl one too many //times, which is why I use counterl-1 below, //note also that there are two blank lines at the //end of the input file, and so this loop proceeds //for one iteration too many } d o u b l e sum_y = 0.0; f o r ( i n t j=discard ; j<counter l - l ; j++) //This loop runs over the data values for generations //that are not desired to be discarded and sums the //differences in the location of the invasion front { sum_y = sum_y + ycoords[ j ] - y c o o r d s [ j - 1 ] ; } i n_s t r eam.c lose ( ) ; r e t u r n ( sum_y/double(counter l - l -d iscard)) ; //returns the average distance between successive //wavefront locations for the non-discarded generations
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Analysis of an integrodifference model for biological...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Analysis of an integrodifference model for biological invasions with a quasi-local interaction Merchant, Sandra M. 2003
pdf
Page Metadata
Item Metadata
Title | Analysis of an integrodifference model for biological invasions with a quasi-local interaction |
Creator |
Merchant, Sandra M. |
Date Issued | 2003 |
Description | The behaviour of a new model for the spatial spread of biological invasions with non-overlapping synchronous generations and well-defined dispersal and sedentary stages is examined. In this integrodifference model, competition between conspecifics takes the form of a quasi-local interaction, where the strength of competition between two individuals depends on their physical distance from each other. Both the deterministic model and a stochastic analogue are examined by numerically simulating the spread of a localized initial population over several generations. By modelling intraspecific competition with a quasi-local interaction, the shape of the travelling waves changed significantly from that of the classical model with only local competition, creating more variable and complex wavefront shapes than are possible with the classical model. The addition of quasi-local competition was also found to alter several aspects of the initial behaviour of this model, including the invasion speed and spatial structure, although in the deterministic case the asymptotic invasion speed and population density behind the front of the wave agreed with those of the classical model. In the stochastic analogue, however, the rate of spread of the invasion was found to be considerably lower than that of the classical model, both initially and asymptotically. Furthermore, the speed achieved by the stochastic invasions was found to depend on the parameters of the quasi-local interaction kernel. |
Extent | 3795157 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-10-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079411 |
URI | http://hdl.handle.net/2429/14518 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2003-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2003-0564.pdf [ 3.62MB ]
- Metadata
- JSON: 831-1.0079411.json
- JSON-LD: 831-1.0079411-ld.json
- RDF/XML (Pretty): 831-1.0079411-rdf.xml
- RDF/JSON: 831-1.0079411-rdf.json
- Turtle: 831-1.0079411-turtle.txt
- N-Triples: 831-1.0079411-rdf-ntriples.txt
- Original Record: 831-1.0079411-source.json
- Full Text
- 831-1.0079411-fulltext.txt
- Citation
- 831-1.0079411.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}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079411/manifest