Space matters: Evolution and ecology in structured populations by Joshua William Anthony Zukewich B.Sc., University of British Columbia, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Master of Science in THE FACULTY OF GRADUATE STUDIES (Mathematics) The University of British Columbia October 2012 c Joshua William Anthony Zukewich, 2012 Abstract The inclusion of spatial structure in biological models has revealed important phenomenon not observed in “well-mixed” populations. In particular, cooperation may evolve in a network-structured population whereas it cannot in a well-mixed population. However, the success of cooperators is very sensitive to small details of the model architecture. In Chapter 1 I investigate two popular biologically-motivated models of evolution in finite populations: Death-Birth (DB) and Birth-Death (BD) processes. Under DB cooperation may be favoured, while under BD it never is. In both cases reproduction is proportional to fitness and death is random; the only difference is the order of the two events at each time step. Whether structure can promote the evolution of cooperation should not hinge on a somewhat artificial ordering of birth and death. I propose a mixed rule where in each time step DB (BD) is used with probability δ (1 − δ). I then derive the conditions for selection favouring cooperation under the mixed rule for all social dilemmas. The only qualitatively different outcome occurs when using just BD (δ = 0). This case admits a natural interpretation in terms of kin competition counterbalancing the effect of kin selection. Finally I show that, for any mixed BD-DB update and under weak selection, cooperation is never inhibited by population structure for any social dilemma. Chapter 2 addresses the Competitive Exclusion Principle: the maximum number of species that can coexist is the number of habitat types (Hardin, 1960). This idea was borne out in island models, where each island represents a different well-mixed niche, with migration between islands. A specialist dominates each niche. However, these models assumed equal migration between each pair of islands, and their results are not robust to changing that assumption. Débarre and Lenormand (2011) numerically studied a two- niche model with local migration. At the boundary between niches, generalists may stably persist. The number of coexisting species may be much greater than the number of habitat types. Here, I derive the conditions for invasion of a generalist using an asymptotic approach. The prediction performs well (compared with numerical results) even for not asymptotically small parameter values (i.e. ≈ 1). ii All parts of this thesis were carried out under the supervision of Dr. Christoph Hauert and Dr. Michael Doebeli. Chapter 1 - Venu Kurella first investigated the mixed Birth-Death and Death- Birth update rule for structured populations in a graduate course. He derived the condition for cooperation being beneficial following the methods in Ohtsuki et al. (2006) for a one- parameter Prisoner’s Dilemma and a one-parameter Snowdrift Game. I used a different method that does not rely on the continuous-population assumption and the subsequent use of ODE’s. My method also modularizes the calculation of fixation probabilities such that the first half of my results (Section 1.1.1) may be used for different updates and different games. The modularization also made it easy to derive results for a general symmetric 2-player game (Appendix B.1-B.5) which I applied to a two-parameter class of games that encompasses all social dilemmas (Section 1.4). I produced all Figures using Matlab and OmniGraffle. I ran all the simulations using the EvoLudo package produced by Christoph Hauert (in Java), with an extension I wrote to accommodate the mixed update rule used here. The results were submitted to the Journal of Theoretical Biology, but it was rejected. Thanks also to helpful comments from the anonymous reviewers at JTB. A manuscript containing the results of Chapter 1 has been submitted to PLoS ONE. Chapter 2 - Dr. Michael Ward provided reading materials and great suggestions, without which I would have made very little headway. I carried out all calculations and used Matlab for all figures and numerical solutions. iii Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Abbreviations and Symbols . . . . . . . . . . . . . . . . . . . . . . . . . x Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Ordering of Birth and Death Updates in Structured Populations . . . 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 The Moran Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Structured Populations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2.1 Birth-Death (BD) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.2.2 Death-Birth (DB) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.3 Mixed DB-BD Update . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3 Well-Mixed Populations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.1 Birth-Death (BD) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3.2 Death-Birth (DB) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.4.1 Cooperation Games . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 iv 1.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2 Invasion of a Generalist at Habitat Boundaries . . . . . . . . . . . . . . . 19 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.2 Two Species Steady State . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Invasion of a Generalist . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4 Matrix Approximation, Numerical Solution . . . . . . . . . . . . . . . . . . 22 2.5 Asymptotics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.6 Exactly Solvable Potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Appendix A Evolutionary Games on Graphs and Pair Approximation . . 33 Appendix B Ratio of Transition Probabilities for Weak Selection . . . . . 36 Appendix C Stability of the Two Species Steady State . . . . . . . . . . . . 40 v List of Tables 1.1 Conditions for C or D being beneficial (columns 2 and 4, respectively) and the condition for C being favoured over D (column 3). Note that BD and DB yield the same conditions in well-mixed populations. k̂ ≡ 1 + 1/k and δ̂ ≡ 1 + δ/k are used for convenience. Note: 1 ≤ δ̂ ≤ k̂. . . . . . . . . . . . . 11 vi List of Figures 1.1 Favoured and beneficial strategies in social dilemmas for well- mixed populations. Parameter space of social dilemmas in well-mixed populations with the cost-to-benefit ratio u as the y-axis and the syn- ergy/discounting parameter v as the x-axis (see Eq. (1.28)). The dashed lines divide the plane into five regions, which correspond to the Prisoner’s Dilemma (PD), Snowdrift Game (SD), Stag-Hunt Game (SH), Deadlock Defection (DD) and Byproduct Mutualism (BM). The three solid lines are predictions for (i) ρD = 1/N - above this line defection is beneficial; (ii) ρC = ρD - below this line cooperation is favoured; and (iii) ρC = 1/N - be- low this line cooperation is beneficial. The three lines intersect at v = 1; for v < 1, cooperation and defection may be simultaneously beneficial, while for v > 1 cooperation and defection may both be detrimental. Each data point represents simulation results for 107 invasion attempts by a single cooperator and 107 invasion attempts by a single defector. Parameters are: selection strength w = 0.05, population size N = 100. . . . . . . . . . . . . 13 vii 1.2 Favoured and beneficial strategies in social dilemmas for struc- tured populations. Parameter space of social dilemmas in structured populations for BD updates (Panel A), DB updates (Panel B), and mixed BD-DB updates (Panel C ). The three solid lines indicate asymptotic pre- dictions based on Pair Approximation for (i) ρD = 1/N , (ii) ρC = ρD, and (iii) ρC = 1/N . In Panel C , DB and BD updates are chosen with equal chances (δ = 0.5). Parameter space organized as in Figure 1. Popula- tion structure is a lattice with connectivity k = 4. The simulation results (for w = 0.05 and N = 100), show good agreement with the analytical predictions (for w 1 and N 1). . . . . . . . . . . . . . . . . . . . . . . 14 1.3 Comparison of update mechanisms in social dilemmas for struc- tured populations. Parameter regions for cooperation in social dilem- mas: predictions for well-mixed versus structured populations. The solid blue lines indicate predictions for ρC = ρD under different updates: (a) DB, (b) mixed BD-DB for δ = 0.5, and (c) BD in structured populations. In well-mixed populations the condition is the same as (c). Population struc- ture can extend the parameter region where cooperation is favored. The shaded area marks the extended parameter region, which has no biological interpretation in this framework. . . . . . . . . . . . . . . . . . . . . . . . . 15 2.1 The potential, V(x), is plotted along with some example E = 0.4. V (x) = E at “turning points.” Here the turning points are at x = ±x∗ where x∗ = ln [(E + 1)/2E]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.2 The vertical axis is the largest eigenvalue of the discretized version of the operator: J(x) + ∆. The original operator was defined on the real line, whereas we define the discretized version on [-10,10], divided into 1000 subintervals. The largest eigenvalue is plotted as a function of α, which is the fitness of the generalist. The critical value is: αc ≈ 0.84. For α > αc the steady state with the generalist extinct is unstable. . . . . . . . . . . . . 23 viii 2.3 The LHS of Eq. (2.21) plotted against E. . . . . . . . . . . . . . . . . . . . 26 2.4 The largest eigenvalue λ that satisfies Eq. (2.22) as a function of α from the WKB approximation. Note that the approximation was made for α 1 though the results are not bad even for α not asymptotically large. The dashed red lines indicate the possible range for λ; the blue line is the solution from the matrix approximation. . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.5 The original potential is given by the black curve, along with various curves of the form Asech2(x/B). Each of the curves is matched to the original potential for two out of the following factors: height at 0, total area, and behaviour for large |x|. These and other exactly solvable potentials could be helpful in generating approximate solutions. . . . . . . . . . . . . . . . . 28 ix List of Abbreviations and Symbols Chapter 1 BD Birth-Death [update or rule]. DB Death-Birth [update or rule]. δ Probability of using DB in the mixed DB/BD update rule. N Population size. k Every individual has exactly k neighbours in a structured population. w Selection strength. Weak selection means w 1. C Arbitrary label before Applications. Specifically Cooperator after. D Arbitrary label before Applications. Specifically Defector after. i Number of C’s T+i The probability that i increases given currently i C players. T−i The probability that i decreases given currently i C players. Π Total payoff to an individual. πCD Payoff a C gets from each D with which it interacts. f Fitness of a C. g Fitness of a D. θi Under weak selection, we define T + i /T − i ≈ 1 + wθi +O(w2). pC , pD Frequency of Cooperators (≡ i/N) or Defectors (≡ 1− pC). qC|D Probability an individual is a C given that it has one D neighbour. ρC , ρD Probability that a single C (or D) replaces the resident population. x The variable number of C neighbours of the focal individual. φ Total fitness of all individuals. Used to normalize the fitnesses. α = πCC − πCD − πDC + πDD, a compound parameter. β = N(πDD − πCD) + πCC + πDD, a compound parameter. BM The Byproduct Mutualism Game (or Deadlock Cooperation). DD The Deadlock Defection Game. SH The Stag-Hunt Game. SD The Snowdrift Game. PD The Prisoner’s Dilemma Game. c Cost a C pays for cooperating. b Benefit that a C player gives, b > c > 0. c0 Cost a C pays for cooperating from the Game used in Ohtsuki et al. (2006). b0 Benefit that a C player gives from the Game used in Ohtsuki et al. (2006). v Synergy/Discount parameter. Each additional benefit received is weighted by v. u = 2c/b, the rescaled cost-benefit ratio, a compound parameter. O Big O Notation. LHS Lefthand side. RHS Righthand side. Chapter 2 i Species types: 1, 2, and 3. x Spatial variable. pi(x, t) Density of species i, a function of space and time. ri(x) Fitness of species i, a function of space. α Fitness of the generalist (species 3) everywhere. r3(x) = α. xi p̂i(x, t) Density of species i at the generalist-extinct steady state. i.e. p̂3 = 0 φi(x) Perturbations away from the generalist-extinct steady state. J Jacobian evaluated at the generalist extinct steady state. λ Growth rate of the perturbation (eigenvalue). V ≡ (2e|x| − 1)−1 analog of the potential well in the quantum tunneling problem. E ≡ (λ+ 1− α)/α analog of the energy in the quantum tunneling problem. WKB Asymptotic method named after Wentzel-Kramers-Brillouin. ≡ 1/α used as an asymptotic parameter: 1 Sn(z) The spatial coefficients of a series used in WKB theory. x∗ Turning points where V = E and solutions change from oscillatory to exponential. LHS Lefthand side. RHS Righthand side. xii Acknowledgements Thanks to my research supervisors, Christoph Hauert and Michael Doebeli, for their support, and to all the delightful staff of the Mathematics Department. Special thanks to Eric Cytrynbaum for introducing me to the field of mathematics by intuition and ensuring I didn’t drop out of academia each year afterward. Much love and thanks to Kelly Paton for introducing me to living like a grown up and ensuring that I didn’t drop out each year afterward. Joshua William Anthony Zukewich The University of British Columbia October 2012 xiii Chapter 1 Ordering of Birth and Death Updates in Structured Populations 1.1 Introduction Evolutionary game theory was developed to model frequency-dependent selection (May- nard Smith, 1982). Arguably the most captivating system that exhibits frequency depen- dent selection is the evolution of cooperation under social dilemmas, a problem that has puzzled researchers across disciplines for decades (for a review, see Doebeli and Hauert (2005)). In a social dilemma, cooperators provide a benefit to a group at some cost to self, while defectors pay no cost and contribute nothing. Groups of cooperators “do better” than groups of defectors, yet in any mixed group defectors “do best” (Dawes, 1980). The tension in social dilemmas is that defection maximizes a given individual’s payoff while cooperation maximizes the total payoff to the group. Network reciprocity is one of many approaches taken to explain the evolution of cooperation (Nowak, 2006). The network describes a spatially or socially structured population (see Nowak and May (1992) for the first such example and Szabó and Fath (2007) for a comprehensive review), in contrast to traditional evolutionary game theory, which assumes that populations are ‘well-mixed’ (where each individual interacts with 1 every other individual with equal likelihood) (Hofbauer and Sigmund, 1988). Evolution on networks is often modeled as a discrete-time birth-death process. Two popular updating rules (see e.g. Ohtsuki et al. (2006); Ohtsuki and Nowak (2006); Ohtsuki et al. (2007); Ohtsuki and Nowak (2008); Fu et al. (2009); Morita (2008); Nowak et al. (2010); Tarnita et al. (2009, 2011); Taylor et al. (2007)) are based on the frequency dependent Moran process (Nowak et al., 2004): • Birth-Death Update (BD): At each update, an individual is chosen for reproduction with a probability proportional to its fitness; its offspring then replaces a randomly selected neighbour. • Death-Birth Update (DB): At each update, an individual is chosen randomly for death; the vacant site is then filled by the offspring of one of its neighbours, selected with a probability proportional to fitness. After many updates, eventually the finite population will be composed of only one type (in the absence of mutation), and we say that this type has reached fixation (further explanation in Section 1.1.1). Note that these updates can also be used for well-mixed populations, which is the special case where every individual is a neighbour of every other individual (the graph is ‘fully connected’). For a structured population, the offspring of a parent is located close to the parent (called ‘limited dispersal’), which plays a crucial role in evolution in structured populations. Ohtsuki et al. (2006) found that spatial structure can promote the evolution of cooperation under DB, but not under the superficially similar rule, BD. In a complemen- tary inclusive fitness approach, Taylor et al. (2007) expressed the same disparity in results between BD and DB. In both rules reproduction is proportional to fitness and death is random, and yet they yield qualitatively different dynamics. This disparity is the focus of the present treatment. The motivation is that whether cooperation is fostered due to population structure should not hinge on an externally imposed ordering of birth and death events. 2 In order to investigate the differences between BD and DB with limited dispersal, we introduce a parameter δ ∈ [0, 1] that allows a smooth transition between the two updating rules. Specifically, in each updating step we use DB with probability δ and BD with probability 1 − δ. For δ = 1 this results in pure DB dynamics, whereas for δ = 0 this results in pure BD dynamics. By varying δ we can then identify qualitative changes in the evolutionary dynamics. 1.1.1 The Moran Process To model evolution in finite populations we use the Moran process (Moran, 1962; Nowak et al., 2004). In each discrete time step one birth and one death occur, so the population size, N , is constant. This assumes that ecological dynamics have come to a steady state in order to focus on evolutionary dynamics. We consider the evolution of a population with two strategies, C and D. The state of the population is the number of C players, i (and N − i D-players). The Moran process is defined by the transition probabilities to go from i → i + 1 C-players (T+i ) and from i→ i− 1 C-players (T−i ). These probabilities depend on how likely it is to interact with either type (using the variable i and information about the population structure) and the fitness effects that result from these interactions (using the parameters from the payoff matrix and the selection strength parameter, w). We define fitnesses as a baseline fitness of 1 − w plus the payoff an individual receives weighted by w. If w = 0, interactions have no effect on fitness and evolution is ‘neutral.’ Under neutral selection T−i /T + i = 1 (the transition probabilities no longer depend on the payoffs). If interactions only have a small effect on fitness, selection is said to be ‘weak’ (w 1), and T−i /T+i ≈ 1 +wθi, where we have introduced the coefficient θi that captures the effects of population structure and update rules. The quantity of interest in finite populations is the probability that a single C eventually replaces a resident D population (or the converse). This is termed the fixation probability of C, ρC (or conversely, ρD). The fixation probabilities can be calculated as 3 (Nowak et al., 2004): ρC = 1 1 + N−1 j=1 j i=1 T−i T+i , (1.1) ρD = 1 1 + N−1 j=1 j i=1 T+N−i T−N−i . (1.2) Under weak selection (w 1), by Taylor expanding Eq. (1.1)-(1.2) in w and ignoring higher order terms we find: ρC ≈ 1N − w N2 N−1 j=1 j i=1 θi, (1.3) ρD ≈ 1N + w N2 N−1 j=1 j i=1 θN−i. (1.4) In the neutral case, in which T+i = T − i , we have θi = 0, and hence the fixation probability is 1/N . We say that C is a beneficial mutation (or simply beneficial) if ρC > 1/N , and C is detrimental if ρC < 1/N . However, for some payoff matrices both ρC and ρD can simultaneously be less than 1/N or greater than 1/N . In these cases we say that C is favoured over D (or simply favoured) if ρC > ρD. For example, if ρC > ρD > 1/N , C and D are both beneficial but C is favoured; or if ρC < ρD < 1/N , C and D are both detrimental but D is favoured. Using Eq. (1.3) and Eq. (1.4), the conditions of interest are: 1/N < ρC : N−1 j=1 j i=1 θi < 0, (1.5) ρD < ρC : N−1 i=1 θi < 0, (1.6) 1/N < ρD : N−1 j=1 j i=1 θN−i > 0. (1.7) In the following we use these conditions for discussing effects of different population struc- tures and update rules. 4 1.2 Structured Populations In a structured population an individual’s fitness depends on who its neighbours are. The transition probabilities depend not only on the number of C’s in the population but also on the detailed configuration of the population – an unwieldy amount of data. To simplify, we first use Pair Approximation (see Appendix A). The local density of C’s around a C is qC|C and indicates the conditional probability that the neighbour of a C is another C. Similarly, qD|C is the conditional probability of finding a D next to a C, and so forth. Finally, pC ≡ i/N denotes the global frequency of C’s. The exact type of network is not of central importance in this treatment. We do assume all individuals have the same number of neighbours (k is a constant). The Pair Approximation method used in our analysis is exact for infinite trees with no loops or leaves. Nevertheless, our numerical results agree well with our analytical approximations, despite the simulated networks being small and containing many loops. In structured populations the weak selection limit (w 1) leads to a separation of timescales (see Appendix A). In the initial phase of the population dynamics the local densities (e.g. qC|C) change quickly while the global densities (e.g. pC) remain approx- imately constant. In a relatively short period, local densities reach a quasi-steady state where individuals are more likely to be surrounded by others of the same type. We can solve for this quasi-steady state by assuming that pC is constant and finding solutions to q̇C|C = 0; i.e. where the local densities are no longer changing (see Appendix A, Eq. A.5). This leads to the quasi-steady state solution: qC|C − qC|D = 1/(k − 1), (1.8) which means that, on average, a C has k/(k − 1) more C’s in its neighborhood than a D has in its neighborhood. Once the quasi-steady state is reached, the next phase of the population dynamics proceeds much slower. Gradually, the global densities change while local densities approximately satisfy the quasi-steady state solution, until eventually one type is lost completely. 5 1.2.1 Birth-Death (BD) For BD updating, an individual is selected from the population to reproduce with a prob- ability proportional to fitness. Its offspring randomly replaces one of its neighbours. Since the fitness of a focal individual depends on its neighbourhood, the transition probabili- ties are averages over all possible neighbourhoods. If we let the focal individual have C-neighbours and k − D-neighbours, then the transition probabilities are: T+i = k =0 pC k qC|Cq k− D|C · f φ · k − k , (1.9) T−i = k =0 pD k qC|Dq k− D|D · g φ · k , (1.10) where φ denotes the total fitness of all individuals in the population, which normalizes the fitnesses in order to use them probabilistically. The first factor in each term of Eq. (1.9) [or Eq. (1.10)] is the probability of finding a focal C [D] with a given neighbourhood; the second factor is its resulting relative fitness; and the third factor is the probability that the offspring replaces a D [C] since death occurs uniformly randomly. f and g are the fitness of a focal C (or D) with C-neighbours and k − D-neighbours (see Appendix B.3). Recall that the ratio of the above transition probabilities (T−i /T + i ) is the crucial quantity to determine the fixation probability of C’s and D’s (see Eq. (1.1) and Eq. (1.2)). In the limit of weak selection, T−i /T + i ≈ 1 + wθi and it is this θi which tells us when cooperation is favoured, according to conditions (1.5)-(1.7). Here, we solve for θBDi (see Appendix B.3), as: θBDi = −α+ k(πDD − πCD)− (k − 2)α i N . (1.11) The πij are the payoffs a type-i gets from each interaction with a type-j. To get Eq. (1.11) we have used the quasi-steady state condition [Eq. (1.8)] and introduced α = πCC−πCD− πDC + πDD for convenience. 6 1.2.2 Death-Birth (DB) For DB updating, one focal individual is randomly selected to die and its neighbours compete to fill the vacant site. The transition probabilities are again averages over the neighbourhoods of the focal individual: T+i = kC pD k qC|Dq k− D|D · kC f̃D kC f̃D + kDg̃D , (1.12) T−i = kC pC k qC|Cq k− D|C · k − g̃C f̃C + (k − )g̃C . (1.13) The first factor in each term of Eq. (1.12) [or Eq. (1.13)] is the probability of finding the focal player (which is selected for death) in each arrangement and the second is the probability that it gets replaced by the opposite type. The f̃j and g̃j are the fitness of C and D neighbours of a focal j individual (see Appendix B.4). As for Birth-Death, we must find θDBi in order to determine when cooperation is favoured. We take the ratio of T−i and T + i and expand as T − i /T + i ≈ 1 + wθDBi , where θDBi = 1 k k2(πDD − πCD) + k(πCD − πCC)− α− α(k − 2)(k + 1) iN . (1.14) Again, we have used the quasi-steady state condition [Eq. (1.8)]. 1.2.3 Mixed DB-BD Update Under the mixed DB-BD rule (DB is used with probability δ and BD with 1− δ) we find that the structural coefficient, θδi , is a weighted average of θ BD i and θ DB i (see Appendix B.5): θδi = θ DB i δ + θ BD i (1− δ). (1.15) This result holds in the limit of weak selection (w 1). In this limit, any mixed update rule behaves the same to zero-th order in w; it is only the first order term where differ- ences arise. Using θδi in conditions (1.5)-(1.7) determines whether C or D are beneficial mutations. 7 1.3 Well-Mixed Populations For contrast and comparison we include the analyses of well-mixed populations under BD and DB. In a population with i C-players, the fitness of a C-player (fi) and a D-player (gi) are: fi =1− w + w N − 1[πCC(i− 1) + πCD(N − i)] (1.16) gi =1− w + w N − 1[πDCi+ πDD(N − i− 1)]. (1.17) 1.3.1 Birth-Death (BD) For BD updating the ratio of the transition probabilities simplifies to: T−i /T + i = gi/fi (Nowak et al., 2004) (and see Appendix B.1), which is approximately 1+wθBDi for w 1, where: θBDi = 1 N − 1 (β − αi) , (1.18) with α = πCC − πCD − πDC + πDD as before and β = N(πDD − πCD) + πCC − πDD. To find when C and D are favoured and beneficial, we insert Eq. (1.18) into con- ditions (1.5)-(1.7) and get (to leading order in 1/N and w): 1/N < ρC : 2(πDD − πCD) < πCC − πDC , (1.19) ρD < ρC : πDD − πCD < πCC − πDC , (1.20) 1/N < ρD : 2(πCC − πDC) < πDD − πCD. (1.21) 1.3.2 Death-Birth (DB) Under DB the individual chosen for death cannot reproduce, and hence the ratio of the transition probabilities is slightly different than under BD updating. We find (see Ap- pendix B.2): θDBi = 1 N (β − αi). (1.22) Note that θBDi = N N−1θ DB i . 8 Substituting into conditions (1.5)-(1.7) we find the same result as for BD. The only difference is the deviation from neutral selection, |ρC − 1/N |, which is larger when using BD rather than DB. In both cases birth is affected by fitness whereas death is uniformly random. In each BD time step an individual with high fitness may be selected for reproduction before it risks being killed. Under DB, an individual with high fitness could be selected to die before it ever has a chance of being selected for reproduction. Because the random step occurs before the step affected by selection, the noise in DB exceeds that in BD - and hence the fixation probabilities under DB are closer to those of an entirely random process, i.e. to 1/N . 1.4 Applications 1.4.1 Cooperation Games Up to this point the types C and D have been arbitrary labels, but the vast majority of game theory has been developed for social dilemmas between cooperators (C) and defectors (D). Social dilemmas are characterized by: (i) two C’s do better than two D’s (πCC > πDD); (ii) interacting with a C is always preferable to interacting with a D (πCC > πCD and πDC > πDD); and finally (iii) a D does better than a C when they interact (πDC > πCD). These restrictions leave four possible orderings of the payoffs (Hauert et al., 2006): πCC > πDC > πCD > πDD (Byproduct Mutualism, BM) (1.23) πCC > πDC > πDD > πCD (Stag Hunt Game, SH) (1.24) πDC > πCC > πCD > πDD (Snowdrift Game, SD) (1.25) πDC > πCC > πDD > πCD (Prisoner’s Dilemma, PD) (1.26) with popular names of the games in parentheses. Note that there is no ‘dilemma’ in BM games since cooperation is trivially favoured. 9 Here, we use a two-player version of the game in Hauert et al. (2006): Cooperators pay a cost c > 0 to provide a benefit b > c to a common pool, which will be equally split between the two players regardless of their strategies. Defectors pay no cost and contribute no benefit. In addition, we let the benefits be non-additive: the first contribution has weight one and the second has weight v. For v > 1 accumulated benefits are synergistically enhanced, whereas for v < 1 the benefit from the additional cooperator is discounted. The payoff matrix for the row player is: C D C b2(1 + v)− c b2 − c D b2 0 (1.27) We normalize the payoff matrix by adding c, then dividing by b/2. After rescaling, the payoff matrix is: C D C 1 + v 1 D 1 + u u (1.28) where u = 2c/b is the adjusted cost-benefit ratio. This payoff matrix encompasses all four social dilemmas: (i) 1 + v > u > 1, u > v Prisoner’s Dilemma; (ii) v > u > 1 Stag-Hunt Game; (iii) 1 > u > v Snowdrift Game; and (iv) u < 1, u < v Byproduct Mutualism. Note, however, that if 1 + v < u then we no longer have that two cooperators do better than two defectors (and so the game is not a “social dilemma”). The resulting game is Deadlock Defection, so called because there is no reason to ever cooperate. Byproduct Mutualism could similarly be termed Deadlock Cooperation. 1.4.2 Results The conditions for selection favouring C or D and for C or D mutations being beneficial are summarized in Table 1.1 and Figures 1.1-1.3 for well-mixed and structured populations under BD, DB and mixed update rules. Table 1.1 illustrates that taking δ → 0 or δ → 1 for the mixed rule recovers the results for the BD and DB updates. The limit k 1 10 Condition: ρC > 1/N ρC > ρD ρD > 1/N Well-Mixed 3u < [v + 2] 2u < [v + 1] 3u > [2v + 1] BD 3u < k̂v + 3− k̂ 2u < [v + 1] 3u > (3− k̂)v + k̂ DB 3u < k̂ k̂v + 3− k̂ 2u < k̂[v + 1] 3u > k̂ (3− k̂)v + k̂ Mixed Rule 3u < δ̂ k̂v + 3− k̂ 2u < δ̂ [v + 1] 3u > δ̂ (3− k̂)v + k̂ Table 1.1: Conditions for C or D being beneficial (columns 2 and 4, respectively) and the condition for C being favoured over D (column 3). Note that BD and DB yield the same conditions in well-mixed populations. k̂ ≡ 1 + 1/k and δ̂ ≡ 1 + δ/k are used for convenience. Note: 1 ≤ δ̂ ≤ k̂. recovers the results for well-mixed populations - a highly connected population behaves like a well-mixed one. In the prisoner’s dilemma (PD) cooperators are favoured (ρC > ρD) on a triangle in the vu-plane with vertices: (1/y∗, 1), (1, 1), and (y∗, y∗), where y∗ = (k + δ)/(k − δ) (bold black bordered in Fig. 1.2 Panel B). This triangle has area 2δ/(k2−δ2) which grows with increasing δ (see Fig. 1.2 Panel B for δ = 1) or decreasing k. Hence, cooperation is enhanced for smaller neighbourhood sizes and when increasing the proportion of DB up- dates. In highly connected populations (k 1) the triangle gets asymptotically small and disappears for well-mixed populations (see Fig. 1.1). The same happens when increasing the proportion of BD updates, and the triangle disappears for δ = 0 (see Fig. 1.2 Panel A). Whenever cooperators are favoured in well-mixed populations (ρC > ρD), they are also favoured in structured populations (so long as v > −1, but the synergy/discounting factor v needs to be positive to be biologically meaningful). Thus, in the limit of weak selection, structure never inhibits cooperation in social dilemmas. A comparison of ρC > ρD in the vu-plane is displayed in Fig. 1.3 for different scenarios. The critical lines (ρC = ρD) for structured and well-mixed populations intersect at (v, u) = (−1, 0) for any value of δ. Within the range of biologically meaningful u, v, structure always promotes the 11 evolution of cooperation. More generally, the conditions for ρC > ρD are: Well-Mixed: πCC + πCD − πDC − πDD > 0, (1.29) Structured: δ(πCC − πCD + πDC − πDD) + k(πCC + πCD − πDC − πDD) > 0. (1.30) Social dilemmas require πCC > πCD and πDC > πDD, and hence if condition (1.29) is satisfied then so is condition (1.30), but the converse is not true in general. The condition ρC > ρD for structured populations under BD is the same as for well- mixed populations, reaffirming that δ = 0 is a critical value for the evolutionary dynamics. Increasing the proportion of DB updates (increasing δ) always makes cooperation more likely to evolve, at least in the limit of weak selection. 12 0. 2 1 2 Co st- Be ne fit Pa ra m et er Synergy/Discount Parameter 0 1 2 Homework 3 J. Zukewich 1 Rock-Paper-Scissors Consider the Rock-Paper-Scissors game with payoff matrix ρC > ρD > 1/N ρC > 1/N > ρD 1/N > ρC > ρD 1/N > ρD > ρC ρD > 1/N > ρC ρD > ρC > 1/N The average payoffs are πR = −sxP + xS ,πP = xR − sxS and πS = −sxR + xP . 1. Show that the average payoff of the population is given by (1− s)(xRxP + xPxS + xSxR). π̄ = xRπR + xPπP + xSπS = xR(−sxP + xS) + xP (xR − sxS) + xS(−sxR + xP ) = −sxRxP + xRxS + xPxR − sxPxS − sxSxR + xSxP = −s(xRxP + xPxS + xSxR) + xRxS + xPxR + xSxP (1− s)(xRxP + xPxS + xSxR) At p = (1/3, 1/3, 1/3), π̄ = (1− s)/3. Replicator dynamics as regularly defined. Due to symmetry, there is a fixed point p = (1/3, 1/3, 1/3). The replicator equation lives on 1 (i) (ii) (iii) Figure 1.1: Favoured and beneficial strategies in social dilemmas for well-mixed populations. Parameter space of social dilemmas in well-mixed populations with the cost-to-benefit ratio u as the y-axis and the synergy/discounting parameter v as the x- axis (see Eq. (1.28)). The dashed lines divide the plane into five regions, which correspond to the Prisoner’s Dilemma (PD), Snowdrift Game (SD), Stag-Hunt Game (SH), Deadlock Defection (DD) and Byproduct Mutualism (BM). The three solid lines are predictions for (i) ρD = 1/N - above this line defection is beneficial; (ii) ρC = ρD - below this line cooperation is favoured; and (iii) ρC = 1/N - below this line cooperation is beneficial. The three lines intersect at v = 1; for v < 1, cooperation and defection may be simultaneously beneficial, while for v > 1 cooperation and defection may both be detrimental. Each data point represents simulation results for 107 invasion attempts by a single cooperator and 107 invasion attempts by a single defector. Parameters are: selection strength w = 0.05, population size N = 100. 13 !"""""""""""""""""""""""""""""""""#"""""""""""""""""""""""""""""""""$" !"""""""""""""""""""""""""""""""""#"""""""""""""""""""""""""""""""""$" !"""""""""""""""""""""""""""""""""#"""""""""""""""""""""""""""""""""$" %&' %&&' %&&&' (( )( *( +, *- *./012.3(&4567/8")919:0801 ;6 48< +0 /0 !8" )9 19 : 08 01 != $" """ """ """ """ """ """ ""# """ """ """ """ """ """ """ """ """ "$ " !" # Homework 3 J. Zukewich 1 Rock-Paper-Scissors Consider the Rock-Paper-Scissors game with payoff matrix ρC > ρD > 1/N ρC > 1/N > ρD 1/N > ρC > ρD 1/N > ρD > ρC ρD > 1/N > ρC ρD > ρC > 1/N The average payoffs are πR = −sxP + xS ,πP = xR − sxS and πS = −sxR + xP . 1. Show that the average payoff of the population is given by (1− s)(xRxP + xPxS + xSxR). π̄ = xRπR + xPπP + xSπS = xR(−sxP + xS) + xP (xR − sxS) + xS(−sxR + xP ) = −sxRxP + xRxS + xPxR − sxPxS − sxSxR + xSxP = −s(xRxP + xPxS + xSxR) + xRxS + xPxR + xSxP (1− s)(xRxP + xPxS + xSxR) At p = (1/3, 1/3, 1/3), π̄ = (1− s)/3. Replicator dynamics as regularly defined. Due to symmetry, there is a fixed point p = (1/3, 1/3, 1/3). The replicator equation lives on 1 Homework 3 J. Zukewich 1 Rock-Paper-Scissors Consider the Rock-Paper-Scissors game with payoff matrix ρC > ρD > 1/N ρC > 1/N > ρD 1/N > ρC > ρD 1/N > ρD > ρC ρD > 1/N > ρC ρD > ρC > 1/N The average payoffs are πR = −sxP + xS ,πP = xR − sxS and πS = −sxR + xP . 1. Show that the average payoff of the population is given by (1− s)(xRxP + xPxS + xSxR). π̄ = xRπR + xPπP + xSπS = xR(−sxP + xS) + xP (xR − sxS) + xS(−sxR + xP ) = −sxRxP + xRxS + xPxR − sxPxS − sxSxR + xSxP = −s(xRxP + xPxS + xSxR) + xRxS + xPxR + xSxP (1− s)(xRxP + xPxS + xSxR) At p = (1/3, 1/3, 1/3), π̄ = (1− s)/3. Replicator dynamics as regularly defined. Due to symmetry, there is a fixed point p = (1/3, 1/3, 1/3). The replicator equation lives on 1 Figure 1.2: Favoured and beneficial strategies in social dilemmas for structured populations. Parameter space of social dilemmas in structured populations for BD updates (Panel A), DB updates (Panel B), and mixed BD-DB updates (Panel C ). The three solid lines indicate asymptotic predictions based on Pair Approximation for (i) ρD = 1/N , (ii) ρC = ρD, and (iii) ρC = 1/N . In Panel C , DB and BD updates are chosen with equal chances (δ = 0.5). Parameter space organized as in Figure 1. opulation structure is a lattice with connectivity k = 4. The simulation results (for w = 0.05 and N = 100), show good agreement with the analytical predictions for w 1 and N 1). 14 PD (-1,0) (a) biologically relevant: u, v > 0 DD SD SH BM -1 0 1 2 0 1 2 Synergy/Discount Parameter Co st- Be ne fit Pa ra m et er (b) (c) Figure 1.3: Comparison of update mechanisms in social dilemmas for struc- tured populations. Parameter regions for cooperation in social dilemmas: predictions for well-mixed versus structured populations. The solid blue lines indicate predictions for ρC = ρD under different updates: (a) DB, (b) mixed BD-DB for δ = 0.5, and (c) BD in structured populations. In well-mixed populations the condition is the same as (c). Population structure can extend the parameter region where cooperation is favored. The shaded area marks the extended parameter region, which has no biological interpretation in this framework. 15 1.5 Discussion and Conclusion We set out to resolve the disparity between Birth-Death and Death-Birth updates for structured populations demonstrated by Ohtsuki et al. (2006). They considered a simple Prisoner’s Dilemma where cooperators pay a cost c0 > 0 to donate a benefit b0 > c0 (sub- scripts used to distinguish from the b, c used in this paper) to their interaction partner and defectors neither provide benefits nor suffer costs. Cooperation is favoured and beneficial under the DB update if b0/c0 > k, where k indicates the average number of neighbours. However, cooperation is never favoured or beneficial under BD. Here, we introduced a mixed rule where at each time step DB is used with prob- ability δ and BD with probability 1− δ. This allowed us to investigate the cause for the qualitative change in the evolutionary dynamics. To compare with Ohtsuki’s work we can substitute their payoff matrix into equations (1.11), (1.14), (1.15) and (1.5)-(1.7). For our mixed BD-DB updating the conditions for cooperators being beneficial and favoured all simplify to: b0 c0 > k δ . (1.31) The condition for DB is recovered for δ = 1 whereas δ → 0 recovers the result for BD. The only qualitatively different outcome on the continuum between BD and DB occurs when using exclusively BD (δ = 0). This suggests that in general, results based on BD updating may not be robust to small changes in the updating procedure. For any δ > 0, there is a critical cost-to-benefit ratio above which cooperation is favoured. This shows that the success of cooperators does not hinge on the sequence of events particular to DB, but is a more general phenomenon. Spatial models capture the effect of limited dispersal, one consequence of which is that individuals are more likely to interact more with others of the same type (called positive assortment) than they would be in well-mixed populations. Positive assortment has a two-fold effect on populations facing social dilemmas: (1) cooperators may achieve a higher fitness through their interactions with other cooperators, but (2) this increased fitness may be for naught if cooperators just replace other cooperators. (1) has often 16 been called ‘kin selection’, while (2) has been termed ‘kin competition,’ a distinction introduced by Hamilton (Hamilton, 1964). Later work by Taylor (1992) showed that (1) and (2) always balance in patch structured populations and hence altruism cannot evolve. However, this balancing does not necessarily occur in network-structured populations (e.g. Nowak and May (1992); Ohtsuki et al. (2006)). Kin selection and competition provide an intuitive framework to understand the differences between BD and DB. Under both rules the fitness of all individuals is calculated before births and deaths and so the effect of kin selection does not depend on whether BD or DB is used. The difference then must lie in kin competition. Under DB, each individual has a 1/N chance to die in each time step and so there is no effect of kin competition. Under BD, the likelihood of an individual dying depends on the fitnesses of its neighbours - the more fit its neighbours are, the more likely the focal individual will be replaced. This means that cooperators, who provide benefits to their neighbours, are actively increasing their own mortality. Under BD the two effects of kin selection and competition exactly counterbalance. This balance represents a critical point in the evolutionary dynamics. Results derived at such a critical point are not robust. The mixed BD-DB update allows variable strength of kin competition: increasing δ decreases kin competition. However, having kin competition outweigh kin selection would require further extension of the model - for instance by letting the interaction and replacement networks be different (Ohtsuki et al., 2007). Our results fit cleanly within the framework introduced by Tarnita et al. (2009). They studied evolution in structured populations by defining a parameter σ that captures the effect of structure and determines which type is favoured in the limit of weak selection. Under our assumptions, we find for the mixed BD-DB update that σ = (k + δ)/(k − δ). To get the σ values for BD and DB in Tarnita et al. (2009) we simply set δ = 0 or δ = 1, respectively. Again, it is apparent that δ = 0 is a critical value as this implies σ = 1: where structure has no effect on strategy selection (Tarnita et al., 2009). The differences between BD and DB stem from a different balance of two biological 17 effects of spatial structure. The mixed BD-DB update provides a transition between the two extremes and highlights that conclusions based on the BD update may not be robust. Finally, we have shown that for the mixed BD-DB update, structure never inhibits the evolution of cooperation. For other update rules, or with strong selection, spatial structure may be detrimental to cooperation, e.g. in the Snowdrift Game (Hauert and Doebeli, 2004; Hauert, 2006; Fu et al., 2010). This highlights the fact that results for structured populations should be explored for robustness to changes in model architecture. 18 Chapter 2 Invasion of a Generalist at Habitat Boundaries 2.1 Introduction The competitive exclusion principle claims that if n habitat types exist in an environment then at most n species can coexist (Hardin, 1960). Any other species that tries to invade one of the habitats must compete with the specialist resident there and will either drive the resident to extinction or be driven extinct itself. However, Débarre and Lenormand (2011) have shown that the competitive exclu- sion principle need not hold; more than n species may exist in n habitat types provided dispersal is limited. The mechanism is termed “habitat-boundary polymorphism.” Near habitat boundaries, specialists in one habitat type may migrate into the other habitat, where they are maladapted. A generalist, that is able to do well in either habitat, but is a specialist for neither, is less prone to this local maladaptation. As a result, generalists can stably persist close to habitat boundaries. In this chapter I provide analysis to supplement the numerical investigation by Débarre and Lenormand (2011). They provide a necessary condition for the generalist to be able to invade: (1) the generalist has to be more fit than either of the specialists, 19 averaged over the whole environment. I derive the fitness where a generalist may invade using an asymptotic approach, and it proves to be significantly higher than (1). I then confirm the results numerically. 2.2 Two Species Steady State For simplicity, we consider only one infinite spatial dimension: the two species live on a line (−∞ < y < ∞). The density of species i is pi(y, t). The densities of the two species evolve according to Eq. (2.1), which describes births/deaths (reaction term) and movement (diffusion) of species i (Débarre and Lenormand, 2011): dpi dt = βpi ri r̄ − 1 +D d2pi dy2 . (2.1) ri(y) is the fitness of species i at location y. The average fitness across the environment is: r̄(y, t) = i pi(y, t)ri(y). β has units of per time. D has units of space squared per time. We can non-dimensionalize time by the substitution τ = βt, and non-dimensionalize space as x = β/Dy. Equation (2.1) in non-dimensional form is: ṗi = pi ri r̄ − 1 +∆pi, (2.2) where the dot represents a derivative with respect to τ and ∆ represents two derivatives with respect to x. Consider the special case when species 1 is a specialist for habitat 1 (x < 0) and species 2 is a specialist for habitat 2 (x > 0). Their fitness is 1 where they are a specialist and 0 otherwise. That is: r1(x) = 1 if x < 00 if x > 0 and: r2(x) = 0 if x < 01 if x > 0 (2.3) We can solve explicitly for the steady states of (2.2) (where ṗi = 0) as: p̂1(x) = 1− 12ex if x < 01 2e −x if x > 0 and p̂2(x) = 12ex if x < 01− 12e−x if x > 0 , (2.4) by assuming that p̂i must be bounded as well as continuous and differentiable at x = 0. This solution is stable (see Appendix C). 20 2.3 Invasion of a Generalist We now introduce a generalist that has a fitness α in both habitats, where 0.5 ≤ α ≤ 1. For some α numerical integration of the PDE’s has shown that the generalist may persist near the habitat boundary (Débarre and Lenormand, 2011). In this section we consider the steady state derived for two species and add the generalist at zero density everywhere. This is a three-species steady state with the gen- eralist extinct (p̂1 and p̂2 from Eq. (2.4), and p̂3 = 0). We derive a stability condition by introducing a small perturbation away from this steady state. The resulting eigenvalue problem is a second order ODE with variable coefficients of the same form as a Schrodinger equation. In general these equations are not solvable, but we can investigate the problem numerically and solve approximate equations exactly. To begin, rewrite the densities as: p1 = p̂1 + e λtφ1, (2.5) p2 = p̂2 + e λtφ2, (2.6) p3 = p̂3 + e λtφ3, (2.7) where the φi are small perturbations and the λ is the growth rate of the perturbation. If λ > 0 for some φ then that perturbation will grow exponentially - ie. the steady state is unstable. Inserting these definitions into equation (2.2) we find: λφ = [J(x) +∆]φ+O(φ2), where (2.8) J(x) = −1 0 −α(1−H(x)) 1−0.5e−|x| 0 −1 −αH(x) 1−0.5e−|x| 0 0 α 1−0.5e−|x| − 1 , (2.9) where H(x) is the Heaviside function and φ = [φ1 φ2 φ3]T . Our problem is reduced to finding the eigenvalues of the operator L = J(x) +∆. The linearization for φ3 only depends on φ3, so all we need to solve is: 1 α d2φ3 dx2 + [V (x)− E]φ3 = 0 (2.10) 21 where V (x) = 1 2e|x|−1 and E = (λ+1−α)/α. We also require that: limx→±∞φ3 = 0 since φ3 is supposed to be a small perturbation. This is a Schrodinger equation, where the potential is V (x) and the energy level is E (Figure 2.2). These equations are exactly solvable for certain V (x), but not in general. The next sections detail different approaches for this particular V (x). !" # $%!& '!" Figure 2.1: The potential, V(x), is plotted along with some example E = 0.4. V (x) = E at “turning points.” Here the turning points are at x = ±x∗ where x∗ = ln [(E + 1)/2E]. 2.4 Matrix Approximation, Numerical Solution The simplest approach is to consider a large but finite domain and discretize the space into n bins. This turns the operator-eigenvalue problem into a matrix-eigenvalue problem which we can solve numerically. We choose the interval [-10,10] and 1000 subintervals (the 22 result is stable to changes in interval size and step size). We get the entire spectrum of eigenvalues, but are only concerned with the largest. We do this for a range of α values and plot the resulting largest eigenvalue against α, shown in Figure 2.2. 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 ï0.5 ï0.4 ï0.3 ï0.2 ï0.1 0 0.1 0.2 0.3 _ lar ge st eï va lue Figure 2.2: The vertical axis is the largest eigenvalue of the discretized version of the operator: J(x)+∆. The original operator was defined on the real line, whereas we define the discretized version on [-10,10], divided into 1000 subintervals. The largest eigenvalue is plotted as a function of α, which is the fitness of the generalist. The critical value is: αc ≈ 0.84. For α > αc the steady state with the generalist extinct is unstable. 23 2.5 Asymptotics What else can we do with equation (2.10)? WKB theory yields good asymptotic solutions, even when the asymptotic parameter is not particularly small. We can use α 1 (even though we are mostly interested in when 0.5 < α < 1). Define = 1/ √ α. We can then write Eq. (2.10) as: 2 d2φ3 dx2 + [V (x)− E]φ3 = 0, (2.11) which admits a WKB-type solution. There are turning points at ±x∗ such that V (±x∗) = E. At turning points the solutions switch from exponential growth/decay to oscillatory behaviour. Without turning points the solutions are not well-behaved, so we require 0 < E < 1 (or α− 1 < λ < 2α− 1). Using WKB theory, we pose an asymptotic solution of the form: φ3 ≈ e 1 ∞ n=0 nSn(x) . (2.12) Substituting (2.12) into (2.11), the leading order equation is [S0(x)]2 = E − V (x), which has solution: S0(x) = ± x x0 E − V (t)dt, (2.13) where the Sn(x) are functions to be determined. We subdivide the domain into regions where E − V (t) > 0 and E − V (t) < 0 separated by the two turning points: ±x∗, where x∗ = −ln[0.5(1 + E−1)] (See Figure 2.1). Around the turning point at −x∗, the leading order WKB solution is: φ3(x) = c1[E − V (x)]−1/4 · exp −1 −x∗ x E − V (t)dt x+ x∗ O(2/3) connecting Airy function x+ x∗ = O(2/3) 2c1[V (x)− E]−1/4 · sin 1 x −x∗ V (t)− Edt+ π4 x+ x∗ O(2/3) (2.14) 24 Around the turning point at x∗, the leading order WKB solution is: φ3(x) = 2c2[V (x)− E]−1/4 · sin 1 x∗ x V (t)− Edt+ π4 x− x∗ O(2/3) connecting Airy function x− x∗ = O(2/3) c2[E − V (x)]−1/4 · exp −1 x x∗ E − V (t)dt x− x∗ O(2/3) (2.15) We require that φ3 is bounded (it must decay at ±∞), so we throw away the terms that are unbounded as |x|→∞. As written, all the radicals are real. This solution behaves as decaying exponentials for |x| > x∗ and for −x∗ < x < x∗ the solution is oscillatory. Near the turning points these solutions break down, and by rescaling space close to the turning points we get Airy functions to connect the oscillatory and exponential behaviours, and hence the same constants in the solutions on either side of the turning points. The region: −x∗ < x < x∗ (away from the endpoints) is described by two equations, one from each turning point analysis. These two must match, that is: 2c1[V (x)−E]−1/4·sin 1 x −x∗ V (t)− Edt+ π 4 = 2c2[V (x)−E]−1/4·sin 1 x∗ x V (t)− Edt+ π 4 (2.16) Or, −c1 sin −1 x −x∗ V (t)− Edt− π 4 = c2 sin 1 x∗ −x∗ V (t)− Edt+ π 2 − 1 x −x∗ V (t)− Edt− π 4 , (2.17) where we used sin(−θ) = −sin(θ) on the LHS, and broke up one integral into two on the RHS. Two terms in the argument of sin on the LHS show up in the arguments of sin on the RHS, therefore (2.17) is satisfied so long as: 1 x∗ −x∗ V (t)− Edt = π 2 + nπ, (2.18) c1 = (−1)nc2, (2.19) 25 where n ∈ I. This restriction gives the values of E which correspond to eigenfunctions that work. Equation (2.18) is: 2 x∗ 0 1 2et − 1 − Edt = π 2 + nπ, (2.20) 2 π x∗ 0 1 2et − 1 − Edt = 1√ α 1 2 + n . (2.21) When E = 0, the LHS of Eq. (2.21) is 1, and when E = 1, the LHS is 0. The LHS is a decreasing, concave up function of E for 0 < E < 1 (See Figure 2.3). If n is any integer other than 0 then the RHS of Eq. (2.21) does not overlap [0,1] for 0.5 ≤ α ≤ 1. n = 0 is the only discrete eigenvalue (NB: as α gets bigger there will be more discrete eigenvalues). 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 LH S of E q. 2 .2 1 E Figure 2.3: The LHS of Eq. (2.21) plotted against E. 26 Our solution is a curve of E − α pairs which satisfy: π 4 x∗ 0 1 2et−1 − Edt 2 = α, (2.22) This integral is a bit unwieldy. We can solve it numerically for some value of E. This is done for a range of E values to get the corresponding α values. Recall that E ≡ (λ+α− 1)/α. We plot λ vs. α in Figure 2.4, along with the possible range of λ values (between the red dashed curves) and the solution from numerically solving the matrix eigenvalue problem (i.e. the curve in Figure 2.2). 0.5 1 1.5 2 2.5 3 −1 −0.5 0 0.5 1 1.5 2 2.5 3 3.5 4 λ α Figure 2.4: The largest eigenvalue λ that satisfies Eq. (2.22) as a function of α from the WKB approximation. Note that the approximation was made for α 1 though the results are not bad even for α not asymptotically large. The dashed red lines indicate the possible range for λ; the blue line is the solution from the matrix approximation. 27 2.6 Exactly Solvable Potentials There are a vast number of potentials for which the time-independent Schrodinger’s equa- tion is exactly solvable. Another method of finding eigenvalues could be to approximate our potential by an exactly solvable one. We would be looking for an even potential. For instance, V (x) = Asech2(x) is a possible approximation, since we have an exact solution to: Ψ+(λ+U0sech2x)Ψ = 0 (Drazin & Johnson, 1989). This or other functions may be of further help, though it may not provide any additional insight that the WKB solution doesn’t already provide. ï5 ï4 ï3 ï2 ï1 0 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1 x V(x) E large |x|, height height, area area, large |x| Figure 2.5: The original potential is given by the black curve, along with various curves of the form Asech2(x/B). Each of the curves is matched to the original potential for two out of the following factors: height at 0, total area, and behaviour for large |x|. These and other exactly solvable potentials could be helpful in generating approximate solutions. 28 2.7 Conclusion The n-habitat n-species concept was based on distinct habitats within which dynamics are well-mixed (Hardin, 1960). In reality, the edges of a habitat are perceived quite dif- ferently than the interior. Even when there are sharp boundaries between habitat types, local migration smooths the transition between them such that a great diversity of per- ceived habitats exist near the boundary. A specialist near the boundary of its habitat is no longer in fact a specialist since it can migrate to spaces where it is extremely mal- adapted. The most fit organism here is a generalist which does sufficiently well whichever habitat it migrates to and may stably coexist with the specialists; this is habitat boundary polymorphism. Here, we have shown that the generalist must be substantially more fit, averaged over the whole domain, than the specialists; much more fit than the sufficient condition given by Débarre and Lenormand (2011). This suggests that while habitat boundaries may be many, the number of generalists fit enough to capitalize on them may be few. Lastly, this problem is also of interest as an example of a widely-studied class of Schrodinger’s equations, which have received much interest over the last forty years. An exact solution to this particular potential may be of interest in that field as well. 29 Bibliography Dawes, R., 1980. Social dilemmas. Annu. Rev. Psychol. 31, 169–193. Débarre, F., Lenormand, T., 2011. Distance-limited dispersal promotes coexistence at habitat boundaries: reconsidering the competitive exclusion principle. Ecology letters 14, 260–266. Doebeli, M., Hauert, C., 2005. Models of cooperation based on the prisoner’s dilemma and the snowdrift game. Ecol. Lett. 8, 748–766. Fu, F., Nowak, M., Hauert, C., 2010. Invasion and expansion of cooperators in lattice populations: Prisoner’s dilemma vs. snowdrift games. J. Theor. Biol. 266, 358–366. Fu, F., Wang, L., Nowak, M., Hauert, C., 2009. Evolutionary dynamics on graphs: Effi- cient method for weak selection. Phys. Rev. E 79, 046707. Hamilton, W., 1964. The genetical evolution of social behaviour. J. Theor. Biol. 7, 1–16. Hardin, G., 1960. The competitive exclusion principle. Science 131, 1292–1297. Hauert, C., 2006. Spatial effects in social dilemmas. J. Theor. Biol. 240, 627–636. Hauert, C., Doebeli, M., 2004. Spatial structure often inhibits the evolution of cooperation in the snowdrift game. Nature 428, 643–646. Hauert, C., Michor, F., Nowak, M., Doebeli, M., 2006. Synergy and discounting of coop- eration in social dilemmas. J. Theor. Biol. 239, 195–202. 30 Hofbauer, J., Sigmund, K., 1988. The theory of evolution and dynamical systems. Cam- bridge University Press, New York, NY, USA. Matsuda, H., Ogita, N., Sasaki, A., Sato, K., 1992. Statistical mechanics of population. Prog. Theor. Phys 88, 1035–1049. Maynard Smith, J., 1982. Evolution and the Theory of Games. Cambridge University Press, New York, NY, USA. Moran, P., 1962. The statistical processes of evolutionary theory. Clarendon Press, Oxford, UK. Morita, S., 2008. Extended pair approximation of evolutionary game on complex networks. Prog. Theor. Phys 119, 29–38. Nowak, M., 2006. Five rules for the evolution of cooperation. Science 314, 1560. Nowak, M., May, R., 1992. Evolutionary games and spatial chaos. Nature 359, 826–829. Nowak, M., Sasaki, A., Taylor, C., Fudenberg, D., 2004. Emergence of cooperation and evolutionary stability in finite populations. Nature 428, 646–650. Nowak, M., Tarnita, C., Antal, T., 2010. Evolutionary dynamics in structured populations. Philos. Trans. R. Soc. Lond. B 365, 19. Ohtsuki, H., Hauert, C., Lieberman, E., Nowak, M., 2006. A simple rule for the evolution of cooperation on graphs. Nature 441, 502–505. Ohtsuki, H., Nowak, M., 2006. The replicator equation on graphs. J. Theor. Biol. 243, 86–97. Ohtsuki, H., Nowak, M., 2008. Evolutionary stability on graphs. J. Theor. Biol. 251, 698–707. Ohtsuki, H., Nowak, M., Pacheco, J., 2007. Breaking the symmetry between interaction and replacement in evolutionary dynamics on graphs. Phys. Rev. Lett. 98, 108106. 31 Szabó, G., Antal, T., Szabó, P., Droz, M., 2000. Spatial evolutionary prisoner’s dilemma game with three strategies and external constraints. Phys. Rev. E 62, 1095–1103. Szabó, G., Fath, G., 2007. Evolutionary games on graphs. Phys. Rep. 446, 97–216. Szabó, G., Tőke, C., 1998. Evolutionary prisoner’s dilemma game on a square lattice. Phys. Rev. E 58, 69–73. Tarnita, C., Ohtsuki, H., Antal, T., Fu, F., Nowak, M., 2009. Strategy selection in structured populations. J. Theor. Biol. 259, 570–581. Tarnita, C., Wage, N., Nowak, M., 2011. Multiple strategies in structured populations. Proc. Natl. Acad. Sci. 108, 2334. Taylor, P., 1992. Altruism in viscous populations - an inclusive fitness model. Evol. Ecol. 6, 352–356. Taylor, P., Day, T., Wild, G., 2007. Evolution of cooperation in a finite homogeneous graph. Nature 447, 469. Van Baalen, M., 2000. Pair approximations for different spatial geometries, in: Dieckmann, U., Law, R., Metz, J. (Eds.), The Geometry of Ecological Interactions: Simplifying Spatial Complexity. Cambridge University Press, New York, NY, USA. 32 Appendix A Evolutionary Games on Graphs and Pair Approximation We model population structure by assigning individuals to the nodes of a graph where edges indicate their interaction partners. Each individual is connected to exactly k neigh- bours, forming a k-regular graph. The fitness of an individual depends on its neighbours, and so the transition probabilities (T+i and T − i ) depend on the configuration of the entire population. In order to simplify the transition probabilities we use Pair Approximation (Mat- suda et al., 1992; Van Baalen, 2000), which reduces the complex information about struc- ture to information about pairs of individuals. Formally, let qm|n be the conditional probability of finding an m-type in the neighbourhood of a focal n-type. Pair Approxi- mation then states that qm|n = qm|n: the conditional probability that a neighbor of the focal n individual is of type m does not depend on any other type individual connected to the n. Pair Approximation is the most common approach, but see Morita (2008) and Van Baalen (2000) for use of triplets or Szabó and Tőke (1998); Szabó et al. (2000) for n-point approximations (n ≤ 6), and Fu et al. (2009) for a numerical approach to improve Pair Approximation. Pair Approximation involves two probabilities (pC and pD) and four conditional 33 probabilities (qC|C , qD|C , qC|D, qD|D) plus four conserved quantities: pC + pD = 1, qC|C + qD|C = 1, qD|D + qC|D = 1 and qD|CpC = qC|DpD. The last equation follows from the fact that the number of C −D edges must be the same as the number of D−C edges. Hence, we have two free variables. We choose the global and local frequency of cooperators pC and qC|C . The expected change per time step in pC is: E(∆pC) ∆t = T+i N − T − i N ≈ −wT + i θi N , (A.1) where we have assumed w 1 such that T−i ≈ T+i (1 + wθi). Note that the T±i depend not only on i (or, more precisely NpC) but on qC|C as well. To find the expected change in qC|C per time step, we determine the number of C−C links that are created or destroyed when a C replaces a D or vice versa. The total number of C −C links is qC|CpC · (Nk/2) (note: there are Nk/2 total links in the population). If a C replaces a D, then the number of C−C links increases by 1+(k−1)qC|D. If a D replaces a C, then the number of C−C links decreases by (k − 1)qC|C . The expected change in the number of C − C links is: Nk 2 E(∆(qC|C · pC)) ∆t = T+i 1 + (k − 1)qC|D − T−i [(k − 1)qC|C ]. (A.2) Using T−i ≈ T+i (1 + wθi), we have: E(∆(qC|C · pC)) ∆t ≈ 2T + i Nk 1− (k − 1)(qC|C − qC|D) +O(w). (A.3) Now we can expand ∆(qC|C · pC) = qC|C∆pC + pC∆qC|C +∆pC∆qC|C = pC∆qC|C +O(w) because ∆pC is O(w) (Eq. (A.1)). Then we solve for ∆qC|C as: E(∆qC|C) ∆t ≈ 2T + i NkpC 1− (k − 1)(qC|C − qC|D) +O(w). (A.4) qC|C changes fast (O(1), see Eq. (A.4)) relative to pC (O(w), see Eq. (A.1)). Hence, we let the fast variable, qC|C , go to its quasi-steady-state while the slow variable, pC , stays approximately constant Ohtsuki et al. (2006). The quasi-steady state is satisfied when: qC|C − qC|D = 1k − 1 . (A.5) 34 Eq. (A.5) (Eq. (1.8) in the main text) states that the updating mechanism leads to more C’s being around a focal C than there are around a focal D. Population structure and limited dispersal provide positive assortment of types. We take Eq. (A.5) to always be satisfied and use it to simplify the transition probabilities. Once we choose a particular update mechanism, the transition probabilities become a function of just the number of C-players (as in the well-mixed case) and we can then use Eq. (1.5)-Eq. (1.7) to find when C and D are beneficial and when C is favoured over D. 35 Appendix B Ratio of Transition Probabilities for Weak Selection For each update rule and structure, the transition probabilities T+i , T − i are different. In each case we determine T−i /T + i in the limit of weak selection (w 1). In this limit T−i /T + i ≈ 1+wθi, where the coefficient θi captures the effect of population structure and update rules. B.1 Well-Mixed Population (BD) Under the BD rule the number of C’s increases by one when a C reproduces and a D then dies: T+i = ifi ifi + (N − i)gi N − i N (B.1) Similarly, the number of C’s decreases if a D reproduces and a C dies: T−i = (N − i)gi ifi + (N − i)gi i N (B.2) Hence T−i /T + i = gi/fi, where fi = 1 − w + w[(i − 1)πCC) + (N − i)πCD] and gi = 1 − w + w[iπDC + (N − i − 1)πDD] are the fitness of C and D given that there are i C-individuals and N − i D-individuals (Nowak et al., 2004). By Taylor expanding and 36 neglecting higher order terms we get: T−i /T + i ≈ 1− w(fi − gi), (B.3) and hence θBDi = −αi+N(πDD − πCD) + πCC − πDD. (B.4) where α = πCC − πCD − πDC + πDD is used for convenience. B.2 Well-Mixed Population (DB) Under the DB rule the number of C’s increases by one when a D dies and a C then reproduces: T+i = N − i N ifi ifi + (N − i− 1)gi . (B.5) Similarly, the number of Cs decreases if a C dies and a D reproduces: T−i = i N (N − i)gi (i− 1)fi + (N − i)gi . (B.6) Hence T−i /T + i = gi fi (i− 1)fi + (N − i)gi ifi + (N − i− 1)gi , (B.7) or, up to first order, T−i /T + i = 1− w(fi − gi) N N − 1 (B.8) and θDBi = N N − 1θ BD i . (B.9) B.3 Structured Population (BD) The transition probabilities under BD for structured populations are given in Eq. (1.9)- Eq. (1.10) based on the fitness of a focal cooperator (fkC ) and a focal defector (gkC ) with kC C-neighbours: fkC = 1− w + w(kCπCC + (k − kC)πCD), (B.10) gkC = 1− w + w(kCπDC + (k − kC)πDD). (B.11) 37 In the limit of weak selection the separation of time scales results in the quasi- steady state condition (A.5) that can be used to simplify Eq. (1.9)-Eq. (1.10): φT+i pCD = 1 + w[πCC + πCD(k − 1)− 1 + (πCC − πCD)(k − 2)pC ], (B.12) φT−i pCD = 1 + w[πDD + πDC(k − 1)− 1 + (πDD − πDC)(k − 2)pD], (B.13) where φ indicates the total fitness of all individuals in the population. Using T−i /T + i ≈ 1 + wθi, we find: θBDi = −α+ k(πDD − πCD)− (k − 2)α i N , (B.14) with α = πCC − πCD − πDC + πDD. B.4 Structured Population (DB) Under DB the transition probabilities are given in Eq. (1.12)-Eq. (1.13) where f̃j and g̃j denote the fitness of C and D neighbours of a focal j individual: f̃j = 1− w + w{[δjC + (k − 1)qC|C ]πCC + [δjD + (k − 1)qD|C ]πCD}, (B.15) g̃j = 1− w + w{[δjC + (k − 1)qC|D]πDC + [δjD + (k − 1)qD|D]πDD}, (B.16) where δj = 1 if j = and δj = 0 otherwise. Using the quasi steady-state condition (A.5), Eq. (1.12)-Eq. (1.13) simplify to: T+i = pCD 1− wξCD k [k − 1− (k − 2)pC ] , (B.17) T−i = pCD 1− wξDC k [k − 1− (k − 2)pD] , (B.18) with ξij = (k − 2)pi(πji − πjj − πii + πij) + k(πjj − πij)− πii + πij . (B.19) Using T−i /T + i ≈ 1 + wθi, we find: θDBi = 1 k k2(πDD − πCD) + k(πCD − πCC)− α− α(k − 2)(k + 1) iN . (B.20) 38 B.5 Structured Population (Mixed Update) Under the mixed update, DB is used with probability δ and BD with probability 1 − δ. The probabilities to increase or decrease the number of C players by one C are: T+δi =T +BD i (1− δ) + T+DBi δ, (B.21) T−δi =T −BD i (1− δ) + T−DBi δ. (B.22) Therefore, θδi = θ BD i (1− δ) + θDBi δ. 39 Appendix C Stability of the Two Species Steady State The two species steady state is: p̂1(x) = 1− 12ex if x < 01 2e −x if x > 0 and p̂2(x) = 12ex if x < 01− 12e−x if x > 0 , (C.1) And ˆ̄r = p̂1r1+ p̂2r2 = 1− 12e−|x|. Define p1 ≡ p̂1+ z1 and p2 ≡ p̂2+ z2 close to the steady state (i.e. |zi| 1). Note that r̄ = r1p1 + r2p2 = ˆ̄r + i ziri We insert these expressions into: ṗi = pi ri r̄ − 1 +∆pi, (C.2) To get: ∂ ∂τ (p̂i + zi) = (p̂i + zi) ri ˆ̄r + i ziri − 1 +∆(p̂i + zi), (C.3) Which simplifies to: żi = p̂i ri ˆ̄r − 1 +∆p̂i + zi ri ˆ̄r − 1 − rip̂i ˆ̄r i ziri ˆ̄r +∆zi +O(z2) (C.4) The term in the square brackets is zero (steady state solution). In the habitat where i is a specialist (i.e. ri = 1 and ˆ̄r = p̂i), we have: żi = zi 1 p̂i − 1 − p̂i p̂i zi p̂i +∆zi +O(z2) = −zi +∆zi +O(z2) (C.5) 40 Alternatively, in the habitat where i is not a specialist (i.e. ri = 0) żi = −zi +∆zi +O(z2) (C.6) Take the perturbation to be of the form: zi = ui(t)eikx, then: u̇i = −ui(1 + k2) +O(), (C.7) The right hand side is always negative, and so the steady state is stable. 41
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Space matters : evolution and ecology in structured...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Space matters : evolution and ecology in structured populations Zukewich, Joshua William Anthony 2012
pdf
Page Metadata
Item Metadata
Title | Space matters : evolution and ecology in structured populations |
Creator |
Zukewich, Joshua William Anthony |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | The inclusion of spatial structure in biological models has revealed important phenomenon not observed in “well-mixed” populations. In particular, cooperation may evolve in a network-structured population whereas it cannot in a well-mixed population. However, the success of cooperators is very sensitive to small details of the model architecture. In Chapter 1 I investigate two popular biologically-motivated models of evolution in finite populations: Death-Birth (DB) and Birth-Death (BD) processes. Under DB cooperation may be favoured, while under BD it never is. In both cases reproduction is proportional to fitness and death is random; the only difference is the order of the two events at each time step. Whether structure can promote the evolution of cooperation should not hinge on a somewhat artificial ordering of birth and death. I propose a mixed rule where in each time step DB (BD) is used with probability δ (1 − δ). I then derive the conditions for selection favouring cooperation under the mixed rule for all social dilemmas. The only qualitatively different outcome occurs when using just BD (δ = 0). This case admits a natural interpretation in terms of kin competition counterbalancing the effect of kin selection. Finally I show that, for any mixed BD-DB update and under weak selection, cooperation is never inhibited by population structure for any social dilemma. Chapter 2 addresses the Competitive Exclusion Principle: the maximum number of species that can coexist is the number of habitat types (Hardin, 1960). This idea was borne out in island models, where each island represents a different well-mixed niche, with migration between islands. A specialist dominates each niche. However, these models assumed equal migration between each pair of islands, and their results are not robust to changing that assumption. Débarre and Lenormand (2011) numerically studied a two-niche model with local migration. At the boundary between niches, generalists may stably persist. The number of coexisting species may be much greater than the number of habitat types. Here, I derive the conditions for invasion of a generalist using an asymptotic approach. The prediction performs well (compared with numerical results) even for not asymptotically small parameter values (i.e. epsilon ≈ 1). |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-10-19 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-ShareAlike 3.0 Unported |
DOI | 10.14288/1.0073348 |
URI | http://hdl.handle.net/2429/43512 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-sa/3.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_fall_zukewich_joshua.pdf [ 966.74kB ]
- Metadata
- JSON: 24-1.0073348.json
- JSON-LD: 24-1.0073348-ld.json
- RDF/XML (Pretty): 24-1.0073348-rdf.xml
- RDF/JSON: 24-1.0073348-rdf.json
- Turtle: 24-1.0073348-turtle.txt
- N-Triples: 24-1.0073348-rdf-ntriples.txt
- Original Record: 24-1.0073348-source.json
- Full Text
- 24-1.0073348-fulltext.txt
- Citation
- 24-1.0073348.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:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0073348/manifest