BROKEN SYMMETRY AND CRITICAL PHENOMENA IN POPULATION GENETICS: THE STEPPING-STONE MODEL. By Timothy Lee Duty B.S., Virginia Polytechnic Institute and State University, 1990 M.Sc, University of British Columbia, 1993 A T H E S I S S U B M I T T E D I N P A R T I A L F U L F I L L M E N T O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F P H Y S I C S A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F B R I T I S H C O L U M B I A September 2000 © Timothy Lee Duty, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I 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 Physics and Astronomy The University of British Columbia 6224 Agricultural Road Vancouver, B.C., Canada V6T 1Z1 zooo Abstract In this thesis, I study the behaviour of the Stepping-Stone model: a stochastic model from theoretical population genetics. It was introduced by Kimura and Weiss [1][2] as a simple model to investigate the interplay of the evolutionary processes of random genetic drift, mutation, migration and selection. In particular, they were interested in the behaviour of spatially-structured populations when these processes were mediated by local interactions. From the point of view of statistical physics, the Stepping-Stone model can be viewed as a 2-species, non-equilibrium reaction-diffusion model with spatial degrees of freedom and unique fluctuations that arise from irreversible processes at the microscale. My thesis begins with a brief overview of the theory of stochastic processes that includes both the classical treatment using master equations, the Fokker-Planck equation and the Langevin equation, and modern formalisms that map these equations to operators and functional integrals. This is followed by a discussion of the steady-state and critical phenomena in the Wright-Fisher and Moran models, the non-spatial predecessors of the Stepping-Stone model. Critical phenomena in these models is shown to be associated with the breaking of a discrete symmetry that results in a discontinuous order parameter. Next, the Stepping-Stone model is introduced and reformulated using operators and path-integrals. Kimura and Weiss were able to solve for the steady-state of the Stepping-Stone model, under certain conditions, but the dynamics of the model remained elusive. This model is related however, to the "Voter" model, which has been well-studied and is known to have a critical spatial dimension of 2, in that only for d < 2 does the system asymptotically reach one of the 2 degenerate absorbing states. I extend the results of Kimura and Weiss and obtain exact results for the dynamics ii of the Stepping-Stone model. Two regimes of the model are analyzed: 1) the neutral regime, where selection has vanished and mutation is viewed as the control parameter. Here a steady state exists, characterized by a correlation length that diverges as the mutation rate, /j, becomes zero; and 2) the broken symmetry or "fixed" state, where selection is small but finite, and mutations so rare that the relevant description concerns the dynamics of "avalanches" or "cascades" of new alleles induced by the initial mutation and perturbing the system from the absorbing state. A unique kind of dynamical critical phenomenon occurs when selective advantage and mutation become negligible. Like the Voter model, it is qualitatively different both above and below a critical dimension dc = 2. For spatial dimension d < 2 , the critical behaviour is associated with the breaking of a discrete symmetry and corresponds to fixation of one of the two alleles (genotypes). Symmetry breaking in the Stepping-Stone model is shown to be a consequence of the asymptotic return probability of a random walk—identically one for d < 2 and strictly less than one for d > 2. In addition to the correlation length, the steady-state of the neutral regime is further characterized by a measure of variance defined as the amplitude of the two-point correlation function evaluated at vanishing separation. In genetics this measure of variance is known as FST, the fixation index. Exact results are derived for both the steady-state value of FST, and its asymptotic time-dependence at the critical point, which approaches 1 in both d = 1 and d = 2. At the critical point, 1 — F S T ~ t"1/2 for d = 1. For d = 2, a much slower decay is found, 1 — FST ~ l / ln(£). The d = 3 critical behaviour of FST is that it approaches a constant C < 1 as (C — FST) ~ t"1/2. The constant C is non-universal and related to the return probability of a random walk. It approaches 1 for very large values of the dimensionless constant K = 27r2^nr i where A - 1 is the spatial scale of the interaction, rg is the generation time, D is the diffusion constant and n is the population density. iii A related infinite-alleles model was studied by both Sawyer [3] and Nagalaki[4]. These authors found similar asymptotics for the probability that two randomly choses indi-viduals are genetically identical under certain assumptions. Sawyer[3] also proved that fixation in the infinite alleles model occurs iff the random walk followed by the individuals is recurrent. The broken symmetry or "fixed" regime of the Stepping-Stone model has been ex-plored from the point of view of survival of rare mutant alleles, here parameterized by a coefficient of selection s. The exponent u, governing the divergence of the relaxation time as s —> 0 is calculated and found to be v = 2 for d — 1, while for d > 2 it is given by the mean field value of v = 1. Two other exponents for critical spreading processes are determined and scaling arguments are presented and used to find the decay exponents that characterize the time-dependent survival probability and its asymptotic value for very long times. These results also establish the upper critical dimension of the model, dc = 2. Finally, these results are rederived and supplemented by a dynamical renormalization group (RG) analysis. The critical behaviour in d = 1 is found in both regimes to be controlled by non-trivial fixed points. The critical behaviour of the broken symmetry regime for d > 2 and its RG fixed point is that of a critical branching process. For the neutral regime, however, the renormalization group flow is qualitatively different in rf = 2 and 3, reflecting the existence of broken symmetry for rf = 2. The RG flow for the rf = 3 neutral regime contains a line of fixed points with the effective description at large scales given by a Gaussian version of the time-dependent Landau-Ginsburg model. Finally, the renormalization results are found to be valid to all orders of perturbation theory. iv Table of Contents Table of Contents v List of Figures x 1 Introduction 1 2 Stochastic Processes. 6 2.1 The Master Equation 6 2.2 The Fokker-Planck Equation 10 2.2.1 The forward equation 10 2.2.2 Solution of the Fokker-Planck equation 12 2.2.3 The backward equation 13 2.2.4 System size expansions: from master equation to Fokker-Planck equation 15 2.3 Langevin Equations 16 2.3.1 A single variable Langevin equation 17 2.3.2 Continuum Langevin equations 18 2.4 Operators and Path Integrals 20 2.4.1 Response functional formalism 21 2.4.2 The master equation in occupation number representation 25 2.4.3 The Fokker-Planck equation using operators 38 2.4.4 System size expansions, again! 44 v 3 Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 49 3.1 The Wright-Fisher and Moran Models 50 3.2 The Fokker-Planck Equation for Moran's Model 51 3.2.1 Expansion of the master equation 53 3.2.2 Operator forms of Moran's model 55 3.3 Heterogeneous diffusion in the Wright-Fisher and Moran models 57 3.4 Mean-field theory 60 3.5 The Effective Potential 61 3.6 The Broken Symmetry State: Fate of new mutants 65 3.7 Discussion 72 4 The Stepping-Stone Model Revisited 74 4.1 Introduction 74 4.1.1 The Exchange Coupled Stepping-Stone Model (ECSS) 77 4.1.2 The Diffusively Coupled Stepping-Stone Model (DCSS) 77 4.2 The Neutral Regime 78 4.2.1 Fourier transformation to wave-vector space 79 4.2.2 Equation of motion for the two-point correlation function. . . . . 81 4.2.3 The steady-state 83 4.3 Broken Symmetry or d < 2 86 4.4 Critical dynamics of the two-point correlation function 86 4.4.1 The decay of heterozygosity 86 4.4.2 The infinite alleles model 91 4.5 Alternative Formulations of the Stepping-Stone Model 93 4.5.1 The Fokker-Planck equation and operators 95 4.5.2 The path-integral 95 vi 4.5.3 The response functional method 96 4.5.4 The neutral regime 97 4.6 The Broken-Symmetry State: Fate of a single mutant gene 98 4.6.1 Diagrammatic Methods in Field Theory. 102 4.6.2 Perturbation theory for the relaxation rate E 103 4.6.3 The effective relaxation rate and "screening" of selection in low dimensions 107 4.6.4 Survival probabilities: scaling arguments 110 5 Renormalization Group Analysis 115 5.1 Introduction 115 5.2 A Simple Example of Renormalization: The biased random walk 116 5.3 Perturbations from the Random Walk Fixed Point 119 5.3.1 n-particle to m-particle processes 119 5.3.2 Relevance of independent decay and branching processes 121 5.3.3 Behavior of coagulation, or annihilation processes 121 5.3.4 A paradox 121 5.4 The Critical Branching Fixed Point 122 5.4.1 Perturbations from the Critical Branching fixed point 124 5.5 The Gaussian Fixed Point 125 5.5.1 Perturbations from the Gaussian Fixed Point 127 5.6 Summary of non-interacting fixed points 128 5.7 Renormalization of the Stepping-Stone Model 130 5.7.1 The Neutral Regime 130 5.7.2 The Broken-Symmetry State 135 5.8 Determining critical exponents from RG eigenvalues 138 vii 6 Summary of Results and Conclusions 141 6.1 The Wright-Fisher and Moran Models 141 6.2 The Stepping-Stone Model 142 6.2.1 The neutral regime 142 6.2.2 The broken symmetry state. 144 6.3 Renormalization Group For Reaction-Diffusion and Related Stochastic Processes 145 6.3.1 Renormalization to all orders 146 6.3.2 Modern Theory of Stochastic Processes 146 Bibliography 148 A Random Walks on the Linear, Square and Cubic Lattices. 152 A. l Random Walk Green's Functions 152 A.2 Relation of the first-passage time to the Green's function 155 A.3 The Critical Dimension of the Random Walk 157 A.4 Asymptotic Behaviour of the Random Walk Green's Function Ga (0,£). . 157 A. 5 The Laplace Transform of Ga (0,t) 158 A.5.1 The Linear Lattice 159 A.5.2 The Square Lattice 159 A.5.3 The Cubic Lattice 160 B Computer Simulations 162 B. l Stochastic Differential Equations 162 B.2 Continuous-time Reaction Diffusion Processes 162 C The Branching Process 165 C l The model 165 vm C.2 Survival probability 166 C.2.1 The method of characteristics 167 C.2.2 The backward equation 169 C.3 Size distribution 171 ix List of Figures 3.1 Schematic diagram of the Wright-Fisher model 52 3.2 Stochastic Processes in the Moran Model 58 3.3 The mean-field potential V (x) 62 3.4 Comparison of the stochastic "effective potential" U (x) with the mean-field potential V (x) 66 3.5 Simulations of the Wright-Fisher model using stochastic differential equa-tions 67 3.6 Discontinuous behaviour of the order parameter at the critical point of the Wright-Fisher model 68 4.1 Computer simulations of the d = 2 Stepping Stone model steady state. . 87 4.2 Numerical solution for the heterozygosity h (i) — 1 — FST (t) 92 4.3 Vertices representing the action for propagation of new mutant alleles. . . 101 4.4 Feynman diagram representing the l-loop correction to the two-point prop-agator 104 4.5 Perturbation expansion for the two-point propagator 105 4.6 Perturbation expansion and Dyson equation for E 106 4.7 Divergence of the relaxation time for d = 1, 2 and 3 109 4.1 Simulation results for the time-dependent survival probability for d = 1. . 114 5.2 Vertices representing the Lagrangian of the Neutral regime 131 5.3 Bubble diagrams contributing to the renormalization of v and u 132 x 5.4 Diagram representing the 1-loop integral that contributes to the differen-tial RG recursion relation 133 5.5 Renormalization group flow for the neutral theory 136 5.6 RG flow for the broken symmetry state 139 C.l A typical tree from a branching process 171 xi Chapter 1 Introduction "Evolution is a stochastic process of change in gene frequencies in natural populations." —Motoo Kimura [5] Population geneticists have long recognized the importance of random or stochastic processes in producing and maintaining genetic variation but have always disagreed on their significance to the mechanisms of evolution. The work of both Seawell Wright [6] and R. A. Fisher [7] marked the beginning of a quantitative theory of the role of genetic drift in the evolutionary process that was outlined by Darwin much earlier in his Origin of Species. What has become known as the Wright-Fisher model now occupies a central role of this theory and forms a starting point for the application of molecular population genetics to natural populations. The Wright-Fisher model has been further generalized and aspects of it find their way into many areas of the biological sciences. Surprisingly, Wright argued that random processes were essential to every aspect of evolution while Fisher took the opposite view that evolution worked by mass selection alone, with random processes playing the minor role of generating rare new variations. Fisher was responsible for "The Fundamental Theorem of Natural Selection", a theory of evolution based upon deterministic equations that has influenced generations of evolutionary biologists. The controversy surrounding the importance of random genetic drift reached a cul-mination with the "Neutral Theory of Molecular Evolution" outlined by Kimura [8] [9]. He proposed that most if not all the observed genetic variation at the molecular level was 1 Chapter 1. Introduction 2 due to recent mutations of selectively neutral alleles balanced by fixation due to random genetic drift. In this thesis I show that from the standpoint of modern statistical physics, his Neutral Theory in the limit of rare mutations can be considered a theory of critical phenomena in population genetics. The significance of this finding is not merely academic as observed mutation rates are in fact small ( ~ 1CT9 - 10"10 per nucleotide per gener-ation), and selective differences, at least at the molecular level are tiny or non-existent [8]. Kimura also investigated the influence of spatial structure on genetic drift by for-mulating the so-called Stepping-Stone model, following Wright's work on what became known as "isolation by distance" or equivalently, correlation of nearby individuals. In the Wright-Fisher model, a predecessor to the Stepping-Stone model, each individual of the population is equally likely to interact with any other. This gives rise to questions concerning the implications of purely local interactions on the evolutionary process (in genetics this is often referred to as spatial population structure ). Wright, argued for the necessity of such spatial subdivision of populations in order to facilitate the "crossing of adaptive valleys" in the fitness landscape. Physicists tend to think of this as "hopping over a potential barrier". His theory, which he christened the "Shifting Balance Theory", was first proposed in 1931 and has met with a stern array of opponents including (of course) Fisher, who down-played the role of small population size and genetic drift in favour of his mass selection. Debate over the relevance of the shifting balance the-ory for the evolutionary processes continues to rage and despite some recent criticisms of Wright's theory, the question is far from being resolved [10] [11][12]. Interestingly enough, Kimura's neutral theory seems to sidestep Wright's Shifting Balance Theory by providing a very different picture of the evolutionary process at the molecular level. However, in both the theories of Kimura and Wright, random processes are fundamental and therefore unfavoured by the opposing camp of determinists. Chapter 1. Introduction 3 Some of the questions concerning spatial structure were answered by Kimura and Weiss [1][2] with their analysis of the Stepping-Stone model. Others, including those related to the shifting balance theory were left unresolved. The work carried out in this thesis describes an attempt to extend their results using some of the ideas and techniques from modern statistical physics and as such, does not directly address the debate over the shifting balance theory1. On the other hand, several new results have been obtained in this study that unambiguously highlight the importance of local interactions and spatial dimensionality on large scale gene flow in distributed populations. These new results are even more relevant for the neutral theory. Hopefully, the extension of the theory to more complicated systems will succeed in addressing ideas such as Wright's shifting balance theory and further elaborating Kimura's neutral theory. With regard to understanding complex many-body populations, it is hoped that the techniques for treating stochastic processes described in this thesis take us a few steps further along our own fitness landscape that represents our understanding of many-body problems in stochastic population models. Kimura realized the need for new mathematical techniques in 1957, "In his theories of evolution Wright put forward an important concept of 'balance,' especially of balance between directional factors such as selection, mutation, and migration and non-directional or stochastic factors such as random sampling of gametes and random fluctuation of environmental con-ditions. It appears that new methods of stochastic processes will be needed for a satisfactory treatment of Wright's theory of evolution.—Motoo Kimura [13] The research contained in this thesis employs several alternative, yet equivalent the-oretical frameworks for treating stochastic processes. These are briefly reviewed in 1That would require a multi-locus, perhaps many allele model in contrast to the single-locus, two-allele model studied in the following chapters. Chapter 1. Introduction 4 chapter 2. Each theoretical approach has its particular advantages and disadvantages and therefore contributes in complementary ways to the overall understanding that has emerged from their synthesis. In this thesis, we will employ all of these to varying degrees in order to complete our understanding of some very fundamental models of population genetics. Perhaps the most intuitive, easiest to understand and oldest of these frame-works is the Fokker-Planck equation which describes the explicit time evolution of the probability density. Unfortunately, its generalization to many degrees of freedom has limited usefulness. The operator formalism for either discrete or continuous variables is much less intuitive and quite unwieldy for describing probability densities. On the other hand, its utility is found in deriving and solving equations for moments and corre-lation functions. Finally, a functional integral can be obtained from either the response functional method or the operator approach and is analogous to the partition function in equilibrium statistical mechanics or the path-integral in quantum field theory. The path-integral descrition of stochastic processes is invaluable for setting up a consistent perturbation theory and renormalization group analysis. In chapter 3 I review the Wright-Fisher and Moran models as a warm-up for under-standing the behaviour of the Stepping-Stone model. The next two chapters contain the main results of this thesis. In chapter 4 exact results for the Stepping-Stone model are derived and discussed. In chapter 5 , I re-examine the Stepping-Stone model from the point of view of the Renormalization Group (RG). The results obtained from the RG are shown to correspond qualitatively and quantitatively to the exact results of chapter 4 . We are fortunate that the Stepping-Stone model admits a few exact, results and is amenable to the different approaches mentioned above. This situation creates a much firmer foundation for the application of these methods to the many, more dif-ficult problems in stochastic population dynamics, reaction-diffusion theory and other non-equilibrium models. Indeed, the Stepping Stone model may be somewhat unique in Chapter 1. Introduction 5 that complementary to the exact results, a non-trivial renormalization group analysis to all orders is possible. The benefits afforded by this observation form a two-way street. Not only can the RG help us understand the Stepping Stone model, the Stepping-Stone model provides a deeper understanding of the application of renormalization group meth-ods to non-equilibrium processes in general, and sheds light on the conceptual foundation of the RG, itself. Chapter 2 Stochastic Processes. This chapter is intended as an overview of stochastic processes and serves to introduce the reader to the terminology and notation used throughout this thesis. First I discuss the classic approaches to working with stochastic processes: the master equation, the Fokker-Planck and Langevin equations. I also point out how these are related and highlight the various equivalences between formalisms. Next we describe several more modern methods of stochastic processes that map the classic approaches to operators, functional integrals or both. These latter techniques, while more abstract and somewhat harder to grasp, are much more valuable for the many-body problems encountered in the analysis of lattice and continuum models. Readers with experience in the description of stochastic processes may only need to skim this chapter to familiarize themselves with the notation. Newcomers to stochastic processes should first read the entire chapter and may need to consult some of the references, even though the presentation is meant to be self-contained. Those whose are familiar with the more traditional methods, the master equation for discrete variables and Fokker-Planck and Langevin equations for continuous variables, should review section 2.4 on the mapping of these equations to operators descriptions and path-integrals. 2.1 The Master Equation. The master equation offers the most flexible description for stochastic processes. Con-sider a stochastic model where the state or configuration space of the system is the set 6 Chapter 2. Stochastic Processes. 7 {C} . The generality of the master equation is that C can represent almost any kind of variable: a configuration of spins on a lattice, a m-tuple (na,rib, ...nm), for example describing the numbers of a m-species system, a set of vertices and edges, a set of binary numbers and so on. The continuous time master equation is a linear equation for the time derivative of the probability Pc (t)—the probability that the system is in configu-ration C at time t. It is written in terms of the transition rates w (C'/C) that give the probability per unit time of jumping from C to C, ^pc (*) = E \-w iClC') Pc (*) - w iC'lC) Po (*)], (2-1) c and simply expresses the conservation of probability: The first term represents the net flow of probability into configuration C and the second, the net flow out of C. We will be concerned with stochastic processes that have the Markov property. This means that the transition rates w(C'/C) depend only on the current configuration as opposed to earlier ones and furthermore are not explicitly time-dependent. As a simple example, we consider the master equation for a nearest-neighbor random walk on the integer line or linear lattice. Here the state space is Z (the integers) and we label a particular configuration as n, the position of our walker. The master equation reads -Pn (t) = a {Pn_! (t) + Pn+1 (t) - 2Pn (t)} , (2.2) where a is the jump or migration rate. The special solution of this equation for the initial condition Pn (0) = 6n;Tn is denoted the random walk Green's function or propagator Ga (n, TO, t). Note that this quantity is really a conditional probability which is commonly denoted P (n, t\m, 0) in probability textbooks. I, however, will use the physics convention here of calling it G. The translational invariance of Eq. 2.2 tells us that the dependence of G on TO and n can only be through their difference, that is Ga (n, TO, t) = Ga {n — m,t). If Qm is the probability that the walker is initially at TO, then the general solution of Eq. Chapter 2. Stochastic Processes. 8 (2.2) is given by pn {t) = ^2 Ga (n - m, t) Qm. m Another example of a master equation would be that describing radioactive decay. Here the configuration space is the positive integers N, and the master equation given by ^Pn (t) = A {(n + 1) P n + 1 (t) - nPn (t)} , (2.3) having a single parameter A, the decay rate. Radioactive decay is a special case of one-step or birth and death processes. These stochastic processes can be modeled as taking place one event at a time so that the transition rates have the form w (n'/n) = A (n) 8n',n-i + 7 W *„',„+!, where A (n) is the death rate and 7 (n) the birth rate. If n is restricted to the positive integers, one usually has what is called a "natural boundary" at n = 0, A (n) = 0. In this case the steady state is readily found, for if dtPn (t) — 0 for all n, one must have A(n)P n = 7 ( n - l ) P n _ i . One can therefore express Pn in terms of Po through p = 7(w - l)7(w - 1) • • - 7(0)p A (n) A ( n - l ) - - . A (1) °' and determine Po from normalization n Birth and death processes are usually solved (when possible) by making use of the generating function 00 F ( ^ ) = 5 > „ ( £ K - (2.4) n=0 Chapter 2. Stochastic Processes. 9 From the generating function, we can find the moments by taking derivatives with respect to z and then setting z — 1 . For example, the expectation of n is dF(z,t) <»> = 2=1 dz Often the moments dmF{z,t) cf>m = (n(n - l)(n - 2) • • • (n - m + 1)) dzr' (2.5) lz=l are referred to as the factorial moments. From the factorial moments, one can always get back the probabilities 1 oo p - • = -> £ (-1) n! ^—' TO! m TO! To illustrate the use of the generating function, we take our example of radioactive decay above and multiply both sides of the master equation in Eq. (2.3) by zn. Next we sum over n to produce the equation !LF{Zit) = -\{z-l)^-zF{z,t). (2.6) Equation (2.6) can be explicitly solved for arbitrary initial conditions using the method of characteristics (see appendix C, section 2.1 for an example of this method). For the present purposes, however, we'll simply look at one particular solution. Suppose we start with a Poisson distribution of particles so that r>n( ~nn nl and hence F(z, 0) = exp [no (z - 1)]. In this case the solution of Eq. (2.6) is easily verified F{z,t) =exp [n0e-xt (z - 1)] . Chapter 2. Stochastic Processes. 10 From this we see that (n (t)) = n0e~xt, and a = y V (t))-(n(t))2 = nQe'xt. Another simple, yet important birth-death process, the branching process, is solved in detail in appendix C and further illustrates the use of generating functions. 2.2 The Fokker-Planck Equation. The Fokker-Planck equation is a generalized "diffusion" partial differential equation (p.d.e) that comes in two slightly different forms. The most common form describes the time evolution of the probability distribution of continuous random variable, which we shall denote as W (x,t). This equation was first used by Fokker [14] and Planck [15] to describe the Brownian motion of particles. A very thorough reference on applications and methods for solving the Fokker-Planck equation is the book by Risken [16]. Risken also shows how the Fokker-Planck equation can be regarded as a special form of the master equation for continuous random variables. 2.2.1 The forward equation. The more general form of the Fokker-Planck equation is also known as the forward Kolmogorov equation and is written not for W (x, t), but in terms of the conditional probability density P (x,t\x0,t0) , which is the probability density for x at time t > t0, given the initial condition W(x,t0) — 6(X — XQ). The distinction between these two forms is necessary at this point because in the discussion that follows we will introduce the backward equation—another equation for the conditional probability P (x, t\xo, to). Chapter 2. Stochastic Processes. 11 Following Risken, we write the forward equation as d d Id2 —P (x, t\x0,t0) = -—A (x) P (x, t\x0, t0) + - — B (x) P (x, t\x0, t0) , (2.7) or defining the linear differential operator / W ( x ) — | A ( x ) + i ^ B ( « ) , (2.8) we can write —P (x, t\x0, t0) = CFP (x) P (x, t\x0, t0), (2.9) We will often refer to the operator Cpp as the Liouville operator of the Fokker-Planck equation. Multiplying Eq. (2.9) by W (xo,to), integrating over Xo, and recognizing that W(x,t) = jP(x,t\x0,t0)-W (x0,t0)dx0, produces the more commonly known Fokker-Planck equation for W (x, t) , -W(x,t) = CFP(x)W(x,t), (2.10) The first term on the right hand side of Eq. (2.8), containing A {x) , is referred to as the "drift" term. The second term, containing B (x), is called the "diffusion" term. In the most simple diffusion processes B (x) is simply a constant. Loosely speaking, the two terms have a relatively simple interpretation. Suppose we start with the initial state above, W (x,to) = 8 (x — x0). A(xo) describes the instan-taneous rate of change of the expectation of x — XQ (i.e. (x) — Xo), while B (x) gives the instantaneous rate of change of the variance of x (i.e. (x2) — XQ). A steady state describing fluctuations around Xo is achieved if A (x0) is such that it focuses trajectories back to XQ, counteracting the spreading behaviour of trajectories caused by the diffusion term. Chapter 2. Stochastic Processes. 12 The Fokker-Planck equation can be generalized to M variables X\...XM — {x} and has the form M a M f) a t=i «,j=i J W({x},t). (2.11) In general, the drift vector Aj ({#}) and diffusion tensor Bij ({x}) depend on all the variables {x}. 2.2.2 Solution of the Fokker-Planck equation. The Fokker-Planck equation is obviously linear in W ({x} ,t), however, two classes of Fokker-Planck equations are often discussed the literature. These go by the names "linear" and "non-linear". This terminology can be somewhat confusing since the two categories are not mutually inclusive and refer to different aspects of the Fokker-Planck equation. This distinction is important to make so it is included in this brief overview. A linear Fokker-Planck equation of the form of Eq. (2.11) is one where the drift vector Ai ({x}) is linear in the variables X\...XM = {x} and the diffusion tensor Bij ({x}) is constant (i.e. independent of {x}, but still a tensor). Linear Fokker-Planck equations have Gaussian distributions for the stationary as well as non-stationary solutions and this is true even if Aj ({x}) is explicitly time-dependent. A non-linear Fokker-Planck equation has a diffusion tensor ({x}) that is explicitly dependent on at least one of the variables of {x}. A more appropriate term for this type of stochastic process is heterogeneous diffusion. Fokker-Planck equations describing heterogeneous diffusion are difficult to solve ex-cept for the single variable case, where a change of variables produces a constant diffu-sion term and a modified drift term. Conceptually one can think of the heterogeneous Chapter 2. Stochastic Processes. 13 diffusion as giving rise to "fluctuation generated forces". For the single-variable Fokker-Planck equation of Eq. (2.10), the steady state distribution Ws (x) can always be found, (2.12) where C is some normalization constant. This allows us to define a potential U (x) , U(x)=\nB(x)-2 f^jjdy so that Ws(x) = Cexp(-U(x)). 2.2.3 The backward equation. The backward Kolmogorov equation, or adjoint equation, as it is sometimes called, is also an equation for the conditional probability P (x, t\xo, to). However, the partial derivatives are with respect to the initial time to and initial "position" XQ. We write it as Q —P(x,t\x0,t0) = £ p P (xo)P(x,t\xo,to), (2-13) Oto and the adjoint operator CpP (xo) corresponding to £pp (x) defined in Eq. (2.8) is given by £ F P (x0,to) = A(x0) ^ + \ B (XO) (2.14) The adjoint operator is obtained by: (i) changing the sign of the drift term, (ii) reversing the order of differentiation and multiplication by A (xo) and B (xo), (iii) changing the partial derivatives with respect to x to those with respect to x0. We should point out that the backward equation is not peculiar to stochastic processes described by a Fokker-Planck equation. The master equation 2.1 is a forward equation. One can also construct a backward equation. However, the relation between forward and backward Ws (x) = C B(x) exp J My) B(y) dy Chapter 2. Stochastic Processes. 14v operators for the master equation is slightly more convoluted than the relation of £ F" p (XQ) to £pp (%) for the Fokker-Planck equation given above. The backward equation is very useful for first-passage time problems. The first-passage time can be defined in a couple of different ways. The first is as the time at which the stochastic variable x (t) first leaves a given domain [17] [18]. For example, suppose the domain is the interval x G [x\, X2]. For realizations of x (t) that start at t = 0, the first-passage time is the time when x (i) reaches either X\ or x2 for the first time. The second definition is more relevant for systems with absorbing states. Absorbing states are present if both A (x) and B (x) vanish for some value of x. In this case, one usually defines the first-passage time as the time at which x (t) hits a particular absorbing state. As an example of the latter, suppose X\ and x2 are absorbing states. The the probability u (XQ, t) of hitting one of the absorbing states in a time interval t satisfies an equation similar to Eq. (2.13), d —u(x0,t)==C£p(x0)u(x0,t). (2.15) The probability for hitting X\ is obtained by solving Eq. (2.15) with the boundary conditions u(x1,t) = 1, u(x2,t) = 0. Conversely, the probability of hitting x2 is obtained from the boundary conditions u(xi,t) = 0, u(x2,t) = 1. Chapter 2. Stochastic Processes. 15 2.2.4 System size expansions: from master equation to Fokker-Planck equa-tion. I have yet to describe how one arrives at a Fokker-Planck description for a given stochastic process. In the following section I discuss one method that relies on the equivalence of the Langevin equation to the Fokker-Planck equation. This approach, however, does not go beyond a merely phenomenological description of the fluctuations. Ideally, one should start with a master equation, written in terms of microscopic configurations and transition probabilities, and systematically derive a Fokker-Planck equation for some sort of coarse-grained variables. Along these lines, van Kampen [19] [20] [21] has formulated a general technique for treating the birth-death and one-step processes described in the previous section which involves an expansion of the master equation in inverse powers of fi, which is a measure of the system size (f2 is usually the volume). The fi-expansion of van Kampen begins with the ansatz that fluctuations in a mi-croscopic variable n should be oc VT2. This ansatz is quite reminiscent of the Central* Limit Theorem for random variables. This situation for a single variable is described by writing n = fl(f){t) + VrLx, (2.16) where 4> (t) is the non-fluctuating "macroscopic" variable and x, a fluctuating continuous variable with probability distribution W (x, t). Van Kampen shows that to O ( f i - 1 ) , 4> (t) obeys a deterministic ordinary differential equation that is precisely the one obeyed by the expectation (n(t)} /O under the stochastic process. The fluctuating part described by W (x, t) satisfies a linear Fokker-Planck equation. Furthermore, the expansion is systematic in that its validity is readily apparent. The appropriateness of the expansion can be briefly summarized by saying that stability of the differential equation for <\> (t) in the standard dynamical systems sense implies the validity of the expansion for large Chapter 2. Stochastic Processes. 16 enough system size. The reader should consult the references above to determine just what "large enough" means [19] [20] [21]. The important point to make here is that van Kampen's method fails near the instabilities, or critical points as they are called in statistical physics. This is unfortunate since these points of the phase diagram are precisely those we would like to study. So the question remains, does there exist a Fokker-Planck representation of the mas-ter equation which is asymptotic in system size but retains an accurate description of the fluctuations, even near critical points? By accurate we mean both qualitatively in that the phase diagram is preserved, and quantitatively in that the critical exponents charac-terizing the instabilities are unchanged. The answer is emphatically yes, as was pointed out by Horsthemke and Brenig [22] [23]. Their work was based on the mathematically rigorous results of Kurtz [24] [25] who proved similar results for the asymptotic represen-tation of stochastic processes by stochastic differential equations (s.d.e.'s). Horstemke and Brenig's results are simple to state, in contrast to the change of variables of Eq. (2.16), one makes the change of variables n = fix. In other words, we transform to the intensive variable x = n/fl. Expansion of the master equation now produces a non-linear Fokker-Planck equation that preserves the character of the fluctuations near unstable critical points. We will illustrate this type of expansion in section 3.2.1 when we derive a non-linear diffusion approximation for Moran's model. 2.3 Langevin Equations. Many stochastic processes are formulated as stochastic differential equations (s.d.e.'s). • Such processes are most easily visualized in terms of ordinary differential equations (o.d.e.'s) subject to random perturbations or "noise". In the physics literature first Chapter 2. Stochastic Processes. 17 order s.d.e.'s are often referred to as Langevin equations due to the work of Langevin [26], who also investigated Brownian motion. 2.3.1 A single variable Langevin equation. A single variable Langevin equation is usually written in the physics literature as *£& = A(x) + y/W(x)V(t), (2.17) where A (x) is the deterministic force and r\(t) , & Gaussian distributed random impulse with zero mean and auto-correlation (n(t)n(t')) = 6(t-t'). (2.18) The main advantage of this formulation is its conceptual simplicity. It is easy to simulate a Langevin equation by first integrating the deterministic part for a time interval At, followed by adding the random "kick" of \JB (x)r\ (t) At. The disadvantages of using Langevin equations are twofold. First of all, Eq. (2.17) as written, is ambiguoiis. This point is discussed in detail by van Kampen [27]. The problem is that n (t) is a stochastic process with discrete jumps. When B (x) is not a constant, one must choose what value of x is used to compute the noise. Should it be the value before or after the jump? Or possibly an average of the two? The Ito interpretation [28] is to use the value before the jump, /t+6t n(t')dt', and is typically appropriate for internal fluctuations, like those of birth-death processes. One can show that this interpretation of Eq. (2.17) is equivalent to the Fokker-Planck equation i t w {x>t] = ~TxA { x ) W { x ' t ) + \ & B { x ) w ( x ' t ] • Chapter 2. Stochastic Processes. 18 Alternatively, one can use the Stratonovich interpretation [29] and take the average x{t + 6t) + x(t)^ r t + s t x(t + 6t) = x(t) + A{x)6t + \lB[ ^ - r ^ - r - - ^ ^ ^ ^ ^ The latter interpretation is appropriate for external noise, for example, a system cou-pled to an environment where the Gaussian noise is never really "white" (delta-function correlated), but has a non-zero relaxation time. One can show that this interpretation leads to the Fokker-Planck equation i w ( *' t ) = ~kA { x ) W { x ' t ) + \ L ^^Vx ^ ^ ) w { x ' f ) • The second limitation of the Langevin equation formulation is that the form of the noise term B (x) must come from a more detailed knowledge of the stochastic process. As pointed out by van Kampen [27], both of these problems relating to specification of the noise are resolved if one starts with a microscopic description in terms of the master equation and proceeds using the system size expansion discussed in the previous section to arrive at a Fokker-Planck equation. One can then use the equivalence described above to get the appropriate Langevin equation and interpretation. This procedure is convenient for simulating a Fokker-Planck equation using the Langevin formulation. 2.3.2 Continuum Langevin equations. The Langevin equation defined by Eq. (2.17) describes a single variable, or degree of freedom, and is easily generalized to several degrees of freedom by using a system of o.d.e.'s supplemented with appropriate auto and cross correlators for the noise terms. This procedure can be further extended to continuous systems with infinitely many de-grees. For concreteness, a large class of continuum Langevin equations for a field cf) (r,t) in d spatial dimensions has the form Chapter 2. Stochastic Processes. 19 (r,i) = F[cj) (r,t)] + y/BWMV (r, t). (2.19) Here the force term F [(f) (r,t)] and the noise strength B [(f) (r,t)] are Junctionals of cf> (r,t). The force term F [(f) (r,t)} is the functional derivative of a potential 50 (r,t) and the Gaussian noise term n (r,t) is specified by (n (r,t) r, (r',t')> = « fa - ri) <5 (r2 - r'2) • • • 8 (rd - r'd) 8 (t - t'). (2.21) In this thesis, (f)T (t) will be synonymous with (f)(r,t). This choice of notation serves to emphasize that r is to be regarded as an index for the stochastic variables 0r (t). For convenience we will often write (f)T and the -^dependence will be implied. The same notation will be used when r is a vector of some d-dimensional lattice. In this case the functional derivative in Eq. (2.20) will simply be a partial derivative, i.e. For many stochastic models, especially those found in statistical physics, the mag-nitude of the noise is independent of (f>(r,t), i.e. B[(f>(r,t)] equals a constant which I denote here T. For these models, the steady-state (t-mdependent) distribution of<<^ P;,-j W [4>r] can be shown to be W[(f)r} = Z-1eW(-U[(f>r}/r), (2.23) and the normalization Z, in the continuum limit is given by the path-integral or functional integral obtained by integrating over all possible realizations of (f)r. In this thesis the notation J [d(f>] will be used to represent functional integrals. Z is also known as the partition function and is given by Z = j [#] exp /T). Chapter 2. Stochastic Processes. 20 As a concrete example we consider the so-called 04 model of equilibrium statistical physics. This model correctly describes the critical behaviour of Ising-like systems. In this case U is the Landau-Ginzburg free energy where t is a linear function of the temperature and u is a constant. It turns out that there is more that one way to extend dynamics to this model. Here we discuss "model A" from the literature of equilibrium critical phenomena [30] [31] (see also chapter 10 of [32]). The thermodynamic force is 2.4 Operators and Path Integrals. The previous sections of this chapter reviewed classical methods for the formulation and solution of stochastic processes. These methods are adequate for the treatment of stochastic models having only a few degrees of freedom but become extremely cumber-some for the analysis of stochastic processes defined on the lattice or in the continuum. Beginning in the 1970's, two different formalisms emerged for the treatment of these types of many-body problems. The response functional formalism is perhaps the most conceptually straightforward and takes as its starting point the lattice or continuum Langevin equation discussed in the preceding section. The Langevin equation is recast as a functional integral from which a consistent perturbation theory and other diagram-matic methods can be applied [33]. For the study of spatial versions of birth-death (2.24) F [<j}r] = D V V r - 2t<j)r - 4u<f>l, (2.25) and to be consistent with equilibrium statistical mechanics, we must set r = kbr. Chapter 2. Stochastic Processes. 21 or one-step processes introduced in section 2.1, a more exotic field-theoretic method is available that uses creation and annihilation operators similar to those of non-relativistic quantum field theory. We give a brief overview of these two approaches in this section. In addition, a third formulation is discussed that involves writing the Fokker-Planck equation in terms of operators that are similar to the creation and annihilation opera-tors for birth-death processes. In contrast to birth-death processes, the operators of this theory represent continuous stochastic variables. This method has not been extensively discussed in the literature and represents the investigation of the author. We show that this is an alternative to the response functional method for obtaining a functional integral but has the added advantage of retaining the operator character when needed. 2.4.1 Response functional formalism. The response functional formalism was introduced by Martin, Siggia and Rose [34] who applied it to the problem of fluid critical behaviour. It was further developed for renormalization group calculations of equilibrium critical dynamics by Janssen [35] and Bausch, Janssen and Wagner [36]. A nice introduction to this method can be found in chapter 10 of the book by Cardy [32]. An extremely thorough exposition of this method and the underlying mathematics is found in the book by Zinn-Justin [33]. Derivation of the functional integral. One starts with a Langevin equation like that given in Eq. (2.19) but now we add a source term hT (t) and for simplicity we will consider the case of constant noise. We have and use this to construct a functional integral. The construction of the functional integral can be viewed as an application of the transformation of random variables. Our (2.26) Chapter 2. Stochastic Processes. 22 aim is to calculate averages and correlation functions of quantities depending on (j)r (t). For example we may want to calculate G<p$, the two-point correlation function which is given by the expectation G^(r,t;r ,,t) = (^P(t)0P.(t')>. (2.27) This means we need to evaluate path-integrals over 0 r (t). Furthermore, we need to include only those values of cj)r (t) that satisfy Eq. (2.26). We move the terms on the right hand side of Eq. (2.26) to the left hand side and call the resulting expression C (r,t). This expression must vanish for all values of <pT (t) that satisfy the Langevin equation. We can accomplish this by writing C (r,i) as the argument of a delta function which results in the dynamic functional Z [hr] = j [dxj>] 6 [C (r,t)} = j [d<t>] [#] exp (- J 0 r (t) C (r,t) ddrd?j (2.28) This requires a functional integral because the delta function constraint has to be applied for each value of r and t. We have introduced the 0 r (t) field for the delta function and the 4>T (t) integrals are taken along the imaginary axis. Now we can average Z over the noise as this is a Gaussian integral for each point r,t. We multiply Z by exp(—n2/2T) and average over nr at each point where we have exp (r)T4)T - ^]dnroc exp (\T(j) drjr oc exp [if^2 and therefore we have produced the functional integral Z[hr] = J [dab] [#]exp ( - J 0 r (t) C (r,t) - ^Tft (t) ddrdt (2.29) The result is that expectations are evaluated with respect to an ensemble weighted by exp (—5 [0r,0r]) , where we have defined the action as S [<f>t, 4>T] = j (^>Tdt4>v - 4>TF [0r] - 4>rhr - ddvdt. (2.30) Chapter 2. Stochastic Processes. 23 (For simplicity we have dropped the explicit time-dependence of the fields). Another term that is often used instead of action is the Lagrangian density, or simply Lagrangian, which is the quantity in Eq. (2.30) that is integrated over space and time. Correlation and response functions. We will need to evaluate expectations of 0 and 0 fields for arbitrary times and positions. For example, the most informative correlation is the density-density correlation function, which we also called the two-point correlation function in Eq. (2.27). The term spatial auto-correlation is often found in the ecological literature and is synonymous with the two-point correlation function. Another fundamental quantity is the G^, the response function, which is defined by the expectation G^(r,t;r',t) = (4>r(t)^(t')), (2.31) and for many-body systems is also called the two-point Green's function, two-point propa-gator, or simply propagator. The response function or propagator of Eq. (2.31) describes the response of the density field (fiT (t) to a delta function perturbation at the point (r',t'). This gives us a nice interpretation for the cpr field introduced in the derivation of the func-tional integral; it is often called the response field. The expectation of any number of cf> or (f> fields can, at least formally if not practically, be found if we introduce another source term, (pThT,to the action in Eq. (2.30), s [<f>r, 4>r] = J ddrdt (trdtfa - 4>rF [<pr] - 4>ThT - <f>X - ^r^j . (2.32) Expectations can now be found by successive functional differentiation of Z [hr, /ir] , Gm,n = (0X1 (t,) 0X 2 (t2) ...0X m (tm) 0 y i (ti) 0 y 2 (t'2) ...0 y„ (t'n)) (2.33) 8 8 8 8 8 8 . _, 8hXl fa)67^ {tiY'-Sh^ (tm)8hyi (t^Shy, (t>2y~8hyn (t'n) L " ' J ' Chapter 2. Stochastic Processes. 24 where h, h = 0 means we have set all h and h fields to zero everywhere. Knowledge of the full set of correlation and response functions is equivalent to a full solution of the Langevin equation combined with the required boundary conditions. The full solution is a probability functional that even in the best case, that of an explicit closed form solution, is usually much to cumbersome for any practical calculation. On the other hand, the correlation and response functions, especially the two-point correlation function and propagator, are invaluable for understanding the behaviour of the stochastic model. Many-body perturbation theory. In general, the action S [</>r,0r] for a given Langevin equation is not sufficiently simple to put Z [/ir, hr\ in a form useful calculating correlation and response functions using Eq. (2.33). The usual method in field theory is to separate S [0r, 0r] into a bi-linear, or Gaussian part, So and the remaining or "interacting" part Si. For a particular correlation or response function, one obtains a perturbation expansion containing correlation and response functions that are evaluated with respect to So as prescribed by Wick's theorem. To simplify the book-keeping, each term of the expansion can be represented by so-called Feynman diagrams. I will not be able to describe this procedure in detail in this thesis. The reader is encouraged to find the details in one of the numerous standard texts covering statistical field theory or zero-temperature, non-relativistic quantum field theory. Some recommendations are the books by Binney, et. al. [37] and Goldenfeld [38]. The reader is cautioned, however, that these books concern the time-independent theory of equilibrium statistical mechanics—there are some subtle differences between this and the field theory described in this section. A much more appropriate treatment for this thesis is the book by Zinn-Justin [33] mentioned in the introduction of this section. A more entertaining exposition of Feynman diagrams can be found in the book by Mattuck [39]. I do illustrate the application of field-theoretic methods for the Stepping-Stone model Chapter 2. Stochastic Processes. 25 in more detail in section 6 of chapter 4. It is unfortunate that such heavy machinery is needed to develop the diagrammatic techniques that in the end are quite simple and almost fun to use. 2.4.2 The master equation in occupation number representation. Master equations for stochastic birth-death processes can be described in terms of an oc-cupation number representation that has been given several other names in the literature including "second quantization representation", "Fock-space methods" and "quantum field methods". These names emphasize the fact that the occupation number represen-tation borrows heavily from the mathematical framework developed for the quantum-mechanical many-body problem. In addition to stochastic population models, these techniques are used for the closely related stochastic models of reaction-diffusion pro-cesses. The occupation number representation was introduced by Doi [40] and used to develop the theory of diffusion-controlled reactions and in particular the reaction-diffusion model dubbed A + A —> B [41]. This model was also studied by Zel'dovich and Ovchin-nikov using similar techniques [42]. The formulation used by these authors was further elaborated by Grassberger and Scheunert who applied it to several birth-death models [43]. Their article presents a good review of the underlying mathematical structure of the theory. The theory has subsequently been applied to many reaction-diffusion models. A recent review of this theory can be found in the article by Mattis and Glasser [44] who also provide an extensive bibliography of many applications and related works. Here we give a brief review of the formalism. Although the formalism is more appropriate for lattice or spatial models having many degrees of freedom, for simplicity of presentation we initially discuss the formalism as applied to a single discrete variable representing the number of particles or individuals of a stochastic birth-death process and follow this with the extension to many degrees of freedom. Chapter 2. Stochastic Processes. 26 The single-variable master equation. In section 2.1 we introduced birth-death processes, where the state of the system is specified by the discrete probabilities, Pn, of having n-individuals. The main idea of the occupation number formalism is to represent the state of the system specified by the Pn's as a vector \ip) in an abstract space spanned by the basis states |n) that correspond to exactly n individuals, n The \n) basis states are normalized by defining the scalar product as (n\m) = n!<5n;7n, (2.35) so that Pn = (n\rj>). (2.36) n\ The choice of scalar product in Eq. (2.35) will be convenient for stochastic processes describing identical particles. Using the vector space we have constructed, the master equation dt takes the form ^ ) = C IV'), (2-38) where we have introduced the Liouville operator C which, using the occupation number representation, is expressed in terms of creation and annihilation operators, a) and a. The o) and a operators are defined by their action on the |ra) basis states, a\n) = n\n — 1) o) In) = In + 1). (2.39) Chapter 2. Stochastic Processes. 27 They obey the "canonical" commutation relations [a,af] = 1, where the commutator [A,B] of two operators, A and B, is defined as it is in quantum mechanics [A,B] — AB - BA. We see that in terms of these operators, \n) = (a))n |0) and the state of zero particles is defined by a |0) = 0. The normalization requirement, 2~2n Pn — lj takes the form (Sty) = 1, (2.40) where \S) is defined by 15) = exp(at)|0) (2.41) n In statistical physics and quantum mechanics, |5") is known as a"coherent state" or "Glauber state". The normalization requirement of Eq. (2.40) places a constraint on the form of C Since (S\if)) — 1 must be preserved by the dynamics, i.e. ^t(S^) = (S\C\^)=0 (2.42) for an arbitrary state l^) . Therefore, one has the requirement (S\ C = 0. (2.43) Chapter 2. Stochastic Processes. 28 This implications of this requirement can be found by considering a process that takes n particles into m particles. Obviously we need a term with w {o^)™ an, where w is the transition rate. The requirement of Eq. (2.43) can be satisfied if we add to this a term with same number of creation and annihilation operators so that the Liouville operator for this process becomes Cm,n = w [(at)m a" - (a f) n an] . (2.44) These two terms correspond to the two terms in the master equation, Eq. (2.37) that represent the flow of probability into and out of the state of n particles. The form of Eq. (2.44) can be verified to satisfy the requirement of Eq. (2.43) by noting that (S\ a+ = (S\. (2.45) Interestingly enough, using the above definitions of operators and states, Eq. (2.38) can be shown to simultaneously be a representation of the master equation, the forward equation and the backward equation. The generating function. We now make contact with the generating function introduced in Eq. (2.4). One can show that F(z) = (z\xP), (2.46) where the coherent state \z) is defined by \z) = exp {zo)) |0). (2.47) Coherent states are also eigenstates of a, a \z) = z \z). Chapter 2. Stochastic Processes. 29 The normalization condition, [S\ip) = 1 becomes r " ( 4 = i = EP» = 1 ' n Using coherent states and the associated generating function, the creation and annihi-lation operators have a simple representation « = | , (2.48) and oJ = z. (2.49) This is a very practical result, for it allows one to very quickly derive the equation of motion for the generating function once C is specified. For example, consider a branch-ing process. This stochastic process is of fundamental importance for understanding stochastic population models. It consists of particles that can independently branch into two particles or decay and otherwise be removed from the system. The Liouville operator for the branching process is defined by £ B P = A (a — a)a) + 7 (a1")2 a — o)a , (2.50) where A is the death rate and 7 is the branching rate. By mere inspection, the generating function is found to satisfy the equation I™ = A(i-4)^>+42«H (2-51) = ( . - l ) ( i . - ^ F l « ) . Expectations of operators. A few remarks concerning expectations of operators are necessary at this point. For any operator A, its expectation is given by a scalar product (S\A\tp). An important Chapter 2. Stochastic Processes. 30 operator is the number operator Af = a^ a. Its expectation is the average number of particles (n) = (S\a1a\ip). (2.52) = {S\a\ip), where the last line follows from Eq. (2.45). Similarly, the k-th moment, pu is given by the expectation Hk = (S\ak\ip). The Schrodinger picture. Just as in quantum mechanics, we can think of our state \ip) as evolving in time and denote the time dependent state as \ip)t. The formal solution of the master equation (2.37) is |tfOt = e r ( t- t o ) l^ >0 • (2-53) The operator Z(t)=ect is the time-evolution operator that takes the initial state to the state at time t. If one can find the eigenvalue spectrum of C, then the model is completely solved. Needless to say, this is only possible for the simplest of birth-death processes. Heisenberg equations. One can also view the time-evolution from the Heisenberg picture. Here the operators evolve in time and expectations are taken with respect to the initial state. The equation of motion for an operator A is given by ^-A(t) = [A,C], (2.54) Chapter 2. Stochastic Processes. 31 and has the solution A(t)=ectA(t)e-ct. (2.55) A n example of Heisenberg equations: the branching process. Consider, once again, the branching process defined above by C^p in Eq. (2.50). First we calculate the average number of particles jta{t) = [a,C]. (2.56) We evaluate the right-hand side of Eq. (2.56) by commuting the a operator all the way v2l to the right. Using a, (a+) d = 2a\ we find a (t) = 2^o) (t) a (t) - (7 + A) a (t). (2.57) dt Taking the expectation, we find | (n(t))=.( 7 -A)(n(t)>, (2.58) which is trivially solved. For more complicated models, the Heisenberg equations couple higher moments. This does not occur for the average number of particles in the branching process since the particles behave independently. Continuing with our example, the branching process, we would like to know Ps (t), the time-dependent survival probability. This is the probability, 1 — P 0 (t), that the system has not yet reached the state of zero particles. Recall that Po (t) is given by the product (0 | e £ B p t | tp} . Suppose we start with only one particle. Then P0(t) = ( 0 | e £ B P V | 0 ) , (2.59) = ( 0 | e £ B F V e ~ £ B P * | 0 ) . Chapter 2. Stochastic Processes. 32 Adding e £ b p * in the product of the last line of Eq. (2.59) is permitted since £ B P is normal-ordered (all annihilation operators are to the right) and therefore We recognize the operator in the middle of Eq. (2.59) as of (t), so P s(t) = l -Po(t ) = l - ( 0 | a + (t)|0). It is convenient to define the operator a = — 1. It has the same commutation relation with a as does a), [a, d] = 1. This gives us Ps(t) = -(0\a(t)\0). Rewriting £ B P in terms of a, £ B P = (7 — A) aa + ja2a. Using the Heisenberg equation for a, da/dt — [a, £ B P ] , we find jtPs (t) = (7 - A) P s (t) - 7 P | {t). The ultimate survival probability P^ is given by the asymptotic value of Pg (t) for long times, Pco = (2.60) 7 Branching processes have the interesting property that for positive growth rate r = 7 —A, there is a finite probability to go extinct, yet the average number grows exponentially! The time-dependence of Ps (t) at the critical point, 7 — A = 0, is also interesting, Mt) = ~ , (2.61) Chapter 2. Stochastic Processes. 33 which is a power law with exponent —1. Away from the critical point, we find * « = : p £ £ _ , (2.62) or an exponential decay to Ps on both sides of the critical point. Many degrees of freedom. The extension of the occupation number formalism to several variables and lattice models is quite straight-forward. All that is necessary is to give indices to the creation and annihilation operators introduced above. For lattice models, we will use the convention of indexing these operators by the vector r which is decomposed into multiples of primitive lattice vectors. For example, on hypercubic lattices (linear, square, cubic, etc.) d = ^2^, (2.63) r where ri; is an integer, and the primitive lattice vector e; is simply the lattice spacing I times the unit vector in the ith direction. Commutators of operators at different sites are zero, (2.64) The normalization state IS*) becomes |S)=ex P ( ^ > r ) 1°) and represents a coherent state at each site. As an example Liouville operator, consider a collection of random walkers that hop to nearest-neighbour sites with migration (jump) rate TO. (Usually in the physics literature, m is used for indices. Here, we bend to the population genetics convention and use TO to denote the migration or hop rate). The Chapter 2. Stochastic Processes. 34 Liouville operator for a random walk is given by £RW = m E E ar+ear + 4 - e a r - 2aJaP, (2.65) r e where the sum over e is a shorthand for the sum over primitive lattice vectors. This model is solved in detail in appendix A using Heisenberg equations to find the Green's function or two-point propagator. One should note that in contrast to the single particle random walk master equation (2.2) presented in section 2.1, the Liouville operator £ R W describes an arbitrary number of random walkers. (In the absence of other interactions, this number is constant of the motion). For lattice models with many sites, and especially in the thermodynamic (infinite lattice) limit, one is usually not interested in obtaining the explicit probability of an arbitrary configuration, nor is this possible for a general C. For the same reasons outlined in the discussion of the response functional formalism, one needs to find correlation functions and propagators. In the present context, correlation functions are expectations of the number of particles at various sites, e.g. the two-point correlation function, G 0 0(r,t;r',t') = (ar (t) ar, (t')) . (2.66) Propagators for birth-death processes describe the response of the system at one point caused by the addition of new particles somewhere else, e.g. the two-point propagator defined by G^(r,t;r',t') = (aT(t)al,(t')). (2.67) Here we are using the Heisenberg picture and the correspondence of 0 (0) with a (a^ ) will become clear shortly. For translationally invariant systems, the spatial dependence in Eq.(2.66) and Eq.(2.67) can only be through the difference r — r'. This is also true for the time-dependence once the steady-state has been reached. Chapter 2. Stochastic Processes. 35 The correlation functions and propagators for interacting systems cannot in general be found in closed form. This is the same situation discussed with regard to the response functional formalism. One way of dealing with this problem is to pass from operators to a path-integral description and follow the same procedure as discussed for the response functional. I discuss this in the following section. An equivalent alternative is to stick with the operators and use the formalism of time-ordered operators, which also leads to diagrammatic expansions. I will not discuss this method in this thesis, but refer the reader to one of the many texts on quantum-field theory. A perennial classic is the book by Abrikosov, Gorkov and Dzyaloshinski [45]. Their discussion of zero-temperature methods is most relevant to the operator formalism discussed here. Before we discuss the path-integral, we point how one can take the continuum limit using the operator formalism. In this limit the creation and annihilation operators obey the standard "Bose" commutation relations a r, at, = 6d(r-r'), (2.68) al,al, = 0, [ar,ar/] = and the normalization state is given by | 5 ) = e x p ^ d r f r a+) |0). From operators to path-integrals. We can construct a functional integral that resembles that of the response functional formalism by considering the time evolution operator using the coherent state represen-tation. We will illustrate this for a single degree of freedom as the extension to many degrees of freedom does not involve any new concepts. A more detailed exposition of this construction in the context of quantum many-body theory can be found in the book Chapter 2. Stochastic Processes. 36 by Negele and Orland [46]. First, we need to know some additional properties of coher-ent states. We recall from Eq. (2.47) that a coherent state \z) is an eigenstate of the annihilation operator with eigenvalue z, a\z) = z\z). (2.69) It is also the left eigenstate of the creation operator, (z\ a) = (z\ z, where z is the complex conjugate of z. Coherent states are orthogonal, but overcomplete. Their scalar product is given by (z'\z) = <*'', and therefore the resolution of the identity operator becomes f d2z 1 = -^e~z2\z) (z\, and the integral is taken over the complex plane. The state at time tf given the initial state at time U can be written W)tf = ^ ~ t o ) \1>)u • (2-7°) Using the coherent state representation, this relates the generating function at tf to that at ti by F (z>, tf) = j ^e-z*Z (z', tf; z, U) F (z, U), with the time evolution kernel Z (z',tf; z,U) = Iz' ,c{tf-ti) Chapter 2. Stochastic Processes. 37 If we divide the time interval tf — ti into M equal steps of size e = tf1Jr, we can write /M-l U dzhdzhe-^ ( z M | e e / > M - i ) < ^ M - i | e £ £ k M - 2 ) - (zi\eeC\z) , k=\ (2.71) where ZM = z and Af=irM. Since C is normal-ordered, i.e. all the creation operators are to the left, and the annihilation operators are to the right, we can therefore let them act on their respective eigenstates in each term of Eq. (2.71) with the result M - l Z(z',tf;z,ti) = Af-1 lim TT dztdzke-*"*"-1 (2.72) M->oo k=l M-l , x 1 ^2 zk I — — J - L (zk,zk) x exp rz(ti)=z' Z(z',tf-z,ti)= / [dz}[dz}e-z'*' exp(-S[z(t) ,z(t)}) , (2.73) Jz(ti)=Z k=l where L {zk) zk) is obtained from the normal-ordered C by replacing all creation operators by the field zk and annihilation operators by the field zk. Taking the limit we have the functional integral fZ(ti)=z' with the action S [z (t), z (t)] = dt {z (t) dtz (t) -L[z(t),z (t)]} . (2.74) This result combined with the initial state given by F {z, U) tells us more than we really need to know in most circumstances. We can imagine turning on the interactions at some finite time, or perturbing the system in some other way and let tf —»• oo, and ti —> — oo, at which time the system will be found in the vacuum state or some other convenient steady-state. This allows us to ignore the initial state at ti and the e~z z term at tf. The only complication is getting the correct "ground state", or steady-state as the case may be. I will not discuss this technical but important problem in this introductory chapter. Chapter 2. Stochastic Processes. 38 If we now add the extra source terms as was done for the response functional formalism, all correlation functions can by found from the functional integral Z [h, h] = j [dz] [dz] exp (-S [z (t) , z (t)]), (2.75) with S [z (t), z (t)] = j dt {z (t) dtz (t) -L[z(t),z (t)}} . As we mentioned earlier, it is not much more difficult to repeat this derivation when one needs to formulate the path-integral for many degrees of freedom. For these prob-lems, we use the notation 0 r (t) and 4>T (t) for the fields that correspond to the creation and annihilation operators of the Liouvillian. To summarize the results, we find that we can construct a path integral description of a stochastic process specified by a Liouville operator C, by introducing the action S [4>T (t), 0 r (t)] = J ddvdt {4>r (t) dt<f>r (t) - L [4>T (t), 0 r (t)] } , (2.76) where L is obtained from C by substituting respectively, </>r (t) and </>r (t) for the creation and annihilation operators aj and a r. 2.4.3 The Fokker-Planck equation using operators. In this section, we develop an operator formalism for the Fokker-Planck equation in a similar manner as was done for the master equation. The operators we introduce here represent continuous random variables, in contrast to the discrete ones of the occupation number formalism. Once again, for simplicity of illustration, we initially consider a single degree of freedom. The extension to many degrees of freedom is obvious and will be discussed after we illustrate use of the formalism with an example. The basic idea is to represent state of system described by the probability density W (x) as a vector \yj) in a suitably constructed Hilbert space. The probability density is Chapter 2. Stochastic Processes. 39 given by a scalar product, W(x) = (x\V), (2.77) so that \tP) = J dx \x) {x\ip) = J dxW (x) \x). (2.78) In this formalism, the basis states \x) correspond to -^functions, (x\y)=6{x-y). (2.79) Normalization requires = Jw{x)dx = 1, (2.80) where we have introduced the normalization state |1), which corresponds to a uniform distribution Q(x) = 1, i.e. |1) = j dx \x). The dynamics must preserve the projection of Eq. (2.80). Furthermore, in this operator formalism, the expectation of an operator A is given by the product (A) = (1\A\1>). We define operators x,p, such that in x-representation {x\x\ip) = xW(x), (2.81) (x\m = ~W(x). (2.82) These operators obey the commutation relations [x,p] = 1. This identification afforded to us by Eq.(2.81) and (2.82) allows one to easily write the operator form of £ F P - F ° r continuous random variables, one often uses the characteristic Chapter 2. Stochastic Processes. 40 function W(k) = J dx e~ikxW (x), (the Fourier transform of W). In the -^representation (k\x\^) = ^W(k), (2.83) (k\p\7p) = kW(k). (2.84) This gives us a quick method to write the partial differential equation for the character-istic function. Expectations of x operators are the moments: the j'-th moment of the distribution, fij is given by the expectation ^ = {i\{xy |^ An additional result that is useful for evaluating expectations is the identity (l|p = 0. Using our Hilbert space, the Fokker-Planck equation reads jt \iP) = CFP | V ) , (2.85) and has the formal solution exp (CFPt) \ip)0 . (2.86) The path integral A this point, we could derive a path-integral from the time-evolution operator. This derivation follows closely that of the occupation number formalism, so I only give the result. One adds the source terms h (t) x (t) and h (t) p (t) to the resulting action and arrive at the functional integral Z [h, h] = J [dp] [dx] exp {-S [p, x]), (2.87) Chapter 2. Stochastic Processes. 41 with action S\p,x]= Jdt{p(t)dtx{t)- L[p(t),x{t)]-h{t)x{t)-h(t)p(t)} , (2.88) where we have obtained L [p (t), x (t)] by substituting x (t) for x, and p (t) for p in £ F P . An example. We now illustrate the formalism with an extended example. We start with a continuous random walk. It is described by the Fokker-Planck equation dW (x) _ 7 d2 dt ~ 2 8X2 and therefore W(x) The action, excluding the source terms, of this process is simply S\p,x] = J dt [p (t) dtx it) - | p 2 (t)] . (2.89) We can carry out the path integral over p (t) since at each point in time it is a Gaussian. The path-integral becomes Z = J [dx]exp(-S'[x]), with the new action (dx(t)/dtf b FJ - — ' which one recognizes as the well-known "Wiener functional". It is typically not conve-nient to perform the integration over [dp]. The only reason here was to make contact with a known result. We can use Heisenberg equations to look at the first two moments. jt(x) = ([x,CFP]) Chapter 2. Stochastic Processes. 42 = T O - I P I V ' } = o, and-±{X>) = ([x2,CFP]) = 7 ( l | l + 2p£|V> = 7-As expected, the variance grows linearly with time. For the next example, we add a damping term. dW (x) . d „ . . , 7 d2 „. . , (2.90) This Fokker-Planck equation describes diffusion in a parabolic potential. Alternatively, if one thinks of x as representing the velocity, Eq. (2.90) is a simple model of Brownian motion. We know from Eq. (2.12) that the steady-state distribution is Gaussian, i.e. In terms of the operators, W ^ = \ l ^ e x p ( ~ ) - (2-91) CFP= - \xp+-p2. Using the Heisenberg equations, we find equations for the first two moments d(x) = -\{x), (2.92) dt d(x2 Solving for the steady state, Pi = (x) = 0, (2.93) Chapter 2. Stochastic Processes. 43 Equations for higher moments can also be found in a similar manner, d~§ = -J^i + T^^Wa- (2-94) Solving these for the steady-state, we have ^ = ^ 0 - 1 ) ^ - 2 . (2.95) Therefore, the odd moments vanish. If we use units such that ^ = 1, we find / L i 2 = 1, //4 = 3,/X6 = 5 • 3,...,fJ,2n = (2n — 1)!!,which are precisely the moments of a Gaussian distribution. The two-point correlation function (or auto-correlation), and two-point propagator can also be found for this model using Heisenberg equations, we find j t P = Ap, (2.96) and d —x = -\x + -yp. (2.97) These equations do not couple x and p to higher products so we can write the solution for the operators, -xt x(t) = * ( 0 ) e ^ - ^ ( 0 ) e — sinh (At) (2.98) p{t) = p(0)eM. (2.99) Using the commutation relations, the fact that (l|p(0) = 0, and taking expectations results in GH (t) = (1 \x (t)p (0)| V) = e~xt, (2.100) and G u (t) = (1 \x (t) x (0)| V) = fi2 (0) e~xt. (2.101) Chapter 2. Stochastic Processes. 44 In the steady-state, Eq. (2.101) becomes GW (*) = Yxe~XK These results can also be found using the functional integral since it can be transformed to the standard Gaussian form. To conclude this section, we note that his formalism is easily extended to continuum models. For example, a spatial version of the example of Eq. (2.90) would be described by the continuum Fokker-Planck equation d W [ 4 > l { t ) ] = j ^ i - D V ^ ^ + x^^wttAt)] dt (2.102) 2 <50R (t) 8cf>r (t) Identifying ^ w i t h 4>r produces the action W [4>r (*)] • S [4>r,4>r] = J ddrdt 4>rdt<PT + 7J>0rV20R - A0 R 0 R + -4>\ (2.103) We find that for any given stochastic process, the operator formalism for the Fokker-Planck equation produces the same functional integral that we could derive using the response functional method. One advantage of the operator formulation is that one is not restricted to only using path-integral techniques. We can alternatively work with operator equations including time-ordered expansions. In addition, as was mentioned previously, we must specify the noise term of the Langevin equation using some additional information. 2.4.4 System size expansions, again! In section 2.2.4, I introduced two system size expansions for birth-death master equations. These expansions provide systematic procedures for passing from a microscopic picture Chapter 2. Stochastic Processes. 45 in terms of discrete random variables to a macroscopic one employing continuous random variables. In other words, they transform stochastic variables representing absolute num-bers of things (extensive variables) into stochastic variables that are densities (intensive variables). The key ingredient for both of these expansions is an appropriate scaling of the microscopic variables by the system size and proper application of the rules for transformation of random variables. The choice of system size scaling, 0 - 1 / 2 vs. fl'1, is also what distinguishes the two expansions. Here we point out that both expansions can be carried out using the occupation number formalism for birth-death processes. This is a very practical result since one is usually concerned with the behaviour in the limit of large system size. We will also find an important connection between the occupation number formalism and the operator form of the Fokker-Planck equation that parallels the classic results of Horsthemke, et. al. [22] [23], and Kurtz [24] [25]. The system size expansion Horsthemke, et. al. produces a non-linear Fokker-Planck equation. The expansion can be carried out using creation and annihilation operators, for example a and a\ by the canonical transformation a fie, (2.104) a = o) — 1 —> c/Q. In the limit of large fl, the c operator now becomes what we denoted as the x operator of the Fokker-Planck equation and the c operator becomes the p operator. The need for the shift to a = a^— 1 can be seen from action of a) and x on their respective normalization states. Recall that for the occupation number representation, expectations are taken with respect to the normalization state \S) = eat |0), (2.105) Chapter 2. Stochastic Processes. 46 which has the property (SI a f = (SI . (2.106) On the other hand, the normalization state |1) for the Fokker-Planck Hilbert space rep-resents a uniform distribution (2.107) and therefore has the property (l\p = 0. (2.108) This suggests the identification above since (S|a = o. This method can be shown to duplicate the results obtained with the classical technique although a more rigorous argument showing how one passes from a finite vector space to the infinite dimensional Hilbert space is desirable. We should mention that Grassberger and Scheunert [43] have discussed many of the implications of using shifted operators like a and the choice of the scalar product of the vector space. Van Kampen's expansion of the master equation [19] [20] can also be implemented with operators. This system size expansion of Van Kampen produces a linear Fokker-Planck equation, and therefore Gaussian fluctuations around the mean-field solution. In terms of operators, it is carried out by the canonical transformation, a Vtic (2.109) a = a) — 1 —> c/VT2. The only complication is that when one makes the identification c —> x, and c —>• p, one must also shift x by an amount equal to (x), i.e. the mean-field solution in order to get the Fokker-Planck operator describing the fluctuations around this value. This is not a Chapter 2. Stochastic Processes. 47 problem, however, for in the limit of large fi, one easily solves the Heisenberg equation of motion for (x) , as we did in the example of 2.4.3. Although this expansion can be quite useful and informative, we will not discuss it further. As we pointed out in section 2.2.4, Van Kampens system size expansion breaks down near the instabilities or critical points of the stochastic process. This is particularly disastrous for models with spatial degrees of freedom. An example: density dependent mortality. Here I illustrate the system size expansion of Eq. (2.104) by extending the branching process discussed in our introduction to the occupation number formalism (section 2.4.2) to include density-dependent mortality. Density-dependent mortality is added to the branching model via the process A + A —->• A, with rate w/fl. The choice of Q in the rate of this process insures that the death rate due to other individuals scales with the density and not (unreasonable) the total number of particles. The Liouville operator has terms for birth, or branching at rate 7 , a finite lifetime A, and the density dependent term with rate w/Q, (2.110) +w (a^a2 — a W ) / a First, we perform the shift to a a' — 1. This produces ja2a + (7 — A) aa (2.111) Now we implement the scaling transformation a (2.112) a a t 1 6/T2, Chapter 2. Stochastic Processes. 48 to find £ = (7- \)cc-wcc2 (2.113) + (7c2 c - wc2c2) /Q,. The non-linear Fokker-Planck equation for this stochastic process is therefore 1 d2 and is an equally valid description of the stochastic process originally formulated via the master equation. In section 3.2.1 we implement this type of system size expansion for a population genetics model. The reader should compare the ease of the approach using operators with that of working directly with the master equation. This section completes our survey of stochastic processes and was included to highlight the relationship between the occupation number representation of birth-death processes and the operator form of the Fokker-Planck equation. I use these ideas again when we look at the equivalent representations of the Stepping-Stone model in chapter 4. Chapter 3 Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model R. A. Fisher [7] was the first to mathematically analyze the problem of genetic drift using a partial differential equation (p.d.e.) known as the Fokker-Planck, or forward Kolmogorov equation. In this chapter, I will primarily use a non-linear Fokker-Planck equation [16] to describe the time-evolution of the probability distribution for a variable x G [0,1]. The variable x represents the frequency of a given gene or allele in a fixed-size population. Fisher's original work employed an equation for y = cos a;, which transforms the non-linear FP equation to a linear one. (This is discussed briefly in section 2.5). In either approach, the Fokker-Planck equation is a continuous-time "diffusion" equa-tion that describes the large system size, long-time limit of what has become known as Wright-Fisher model. This model is of fundamental importance in theoretical population genetics and provides a basis for the quantitative tools of modern phylogenetic analysis. Even though this model has been quite well studied over the years, understanding its be-havior is an important prerequisite for further development of the Stepping-Stone model and stochastic population dynamics in general. Methods other than the FP equation can be used to study the Wright-Fisher model. However, the FP equation would be most familiar to physicists (it resembles the Schrodinger equation) and more closely related to developments in modern statistical mechanics. 49 Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 50 3.1 The Wright-Fisher and Moran Models. The Wright-Fisher model is discrete, both in time and state space. It consists of a finite population of N haploid individuals which are either of type 'A' or 'a'1. Using genetics terminology, this would be called a "one locus, two allele" model. The population evolves under non-overlapping generations with sequential stages of producing an essentially infinite number of gametes (replication), followed by sampling of these gametes to reform a population of size N. Thus the total population size is conserved. Genetic drift, or fluctuation in the allele frequencies, results from this random sampling which is easily shown to be of binomial form. A schematic of the Wright-Fisher model is shown in figure 3.1. In the limit of large population size, the state of the system is described by the gene frequency x=^ (where na is the number of "a" individuals and we equally well could have chosen A as the allele of interest. Both Fisher [48], and Wright [6] studied the effects of this type of binomial sampling on gene frequencies for finite populations within the context of the discrete model. Subsequently Wright [49] introduced the diffusion approximation which lead to a non-linear FP equation that has been well-studied in the theoretical population genetics literature. A simpler alternative to the Wright-Fisher model was presented by Moran [50], and consists of overlapping generations. A schematic of Moran's model is presented in figure 3.2. In this model, two haploid individuals are picked at random. One dies without reproducing while the other is replicated. The basic effect is the same as the binomial sampling in the Wright-Fisher model. That is, not all individuals contribute to the next generation. Moran's model can be formulated in term of a birth-death master equation xThe convention in population genetics is to specify the population size as 2N. The N refers to the number of diploid individuals, each formed by the union of 2 haploid gametes. Under random mating and after one generation, the diploid population will be in Hardy-Weinberg equilibrium, whereby the diploid population behaves as a haploid population with twice the size[47]. To simplify the calculations, we use N for the number of haploid individuals. One can convert the results to diploid populations under random mating by substituting 2N for N. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 51 of the kind discussed in chapter 2. In the next section, start from this description of Moran's model and derive a non-linear Fokker-Planck equation for Moran's model which will be identical to that of the Wright-Fisher model. 3.2 The Fokker-Planck Equation for Moran 's Model. In section 2.2.4, I presented an asymptotic expansion of the master equation that pro-duces a non-linear Fokker-Planck equation. I use this expansion here derive a non-linear diffusion equation for Moran's model. First I carry out the expansion with the classical technique. Although I use the system size expansion of Horsthemke and Brenig [22] [23], the notation and style is closer to that of Van Kampen [19] [20]. The dynamics of Moran's model consist of two types of random events that occur independently in time. This means the stochastic process is Markovian and the time intervals between similar events are exponentially distributed. The first event type cor-responds to random sampling of gametes with a selective advantage parameterized by s. Two individuals are picked at random. One produces an offspring identical to itself while the other is discarded without reproducing. This is illustrated by the 'replacement' processes 1 + s A + a -> 2a, rate , (3.1) 1 — s A + a -»• 2A, rate , (3.2) where Af is the total number of haploids, N = na + HA- The processes 2A —> 2A and 2a —»• 2a do not change the configuration and hence do not enter the stochastic description. The A -^dependence of the rates in Eqs. (3.1) and (3.2) is chosen so that the replacement rate per individual becomes independent of N. Hence the concept of a "generation time", r g is meaningful. The second event type describes mutation. These Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 52 Generation Size Figure 3.1: Schematic diagram of the Wright-Fisher model. Each generation of size 7Y replicates into an infinite pool of gametes from which N are sampled to form the next generation. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 53 events are described by the processes a —> A, A —> a. We assume mutations are reversible and both events occur at rate //. In the expansion below, N will be the system size parameter we called Cl in chapter 2. In general, Cl is not the total number of particles, but the volume. For this model the total number is a conserved quantity and therefore we can take Cl = N. 3.2.1 Expansion of the master equation. With the stochastic events described above, the master equation for this process is +/x { (£^6a -l)na + ii (£-x£A -l)nA}P (na,nA,t) (3.4) where £ a is a "shift" operator defined by Zaf {na,nA) = f [na + l,nA), and C 1 / (na,nA) = f{na- l,nA), and likewise for the £'As. Now we transform to intensive quantities _ria nA x a - jy ; xA N . The shift operators under this transformation become Sa = Y = V N-*. (3.5) Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 54 Using the rules for transformation of random variables, we find the probability distribu-tion W (xa, xA; t) = P (na, nA, t) / d(xa,xA) d{na,nA) (3.6) where is the Jacobian of the transformation d(na,nA) d(xa,xA) d{na,nA) dxa dxa dna dnA dxA dxA dna dnA = N~2. If we now write the equation W (xa, xA; t) and use the shift operators for the expansion, we find no terms of order i V 0 . The highest order terms are O (N"1) , keeping only these gives us the equation dW(xa,xA,t) p ( d d \ . a ( d d \ . . - Q ^ — — = ^ \ ~ - T r - ) X A W { x a , x A , t ) + ^[ — -—)xaW{xa,xA,t) 2 \dxa dxA 2 \dxA dxa ltd2 Ar d AT d 82 0 d2 , .... + sN- sN^- + - 2-—— } xaxAW {xa, X A , but since 2Nrg \dxl dxa dxA dx\ dxAdxa xa + X A — 1, dxa = -dxA, we can define x = xa, so that d/dxa = -d/dxA = d/dx, and the equation becomes d W ^ / ] = ^(2x-l)W(x,t) + sr-1^x(l-x)W(x,t) (3.7) 1 d2 2NTz dx2 x(l-x)W{x,t). which is also the diffusion equation for the Wright-Fisher model [47]. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 55 3.2.2 Operator forms of Moran's model. Using the operator formalism discussed in section 4.2 of chapter 2, we can write the master equation for Moran's model in terms of creation and annihilation operators. This approach will be generalized in the next chapter as a reformulation of the Stepping-Stone model. We use and a to describe creation and annihilation of the 'a' individuals, and &tand b for the 'A' individuals. The Liouville operator reads 1 f(l-s) (6 f) 2a6+(l + s) (a*)2 a& - 2at6ta&| (3.8) C = 2NTS +p (b]a - o)a) + fi (a+6 - b]b) . Due to the conservation law, one might be tempted to replace the b operators with numbers, b -» N-a, (3.9) 6+ -> 2 - a f . and write a Liouville operator in terms of the a's only. There is a problem with this, however. It does not represent a canonical transformation. The transformation of Eq. (3.9) is equivalent to operator transformations d = a + b, (3.10) d) = at + tf, followed by the replacement d —»• N, and d^ —»• 2. The first problem is that [d,^]=2, (3.11) which is easily remedied. A much more serious problem arises from the commutator [d,a^] =1^0. (3.12) Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 56 The correct way to handle the problem is illustrated in two successive transformations. We could combine the two, but i) it is more instructive to proceed in two steps, ii) we use the separate transformations in the next chapter when we generalize this model to spatially-distributed populations, and iii) the first step is carried out in a manner that facilitates the operator form of the system size expansion discussed in chapter 2. First, we make a canonical transformation to sum and difference operators, c=(a-b)/N, d = a + b, ct = f (at -tf), dt = I ( at + &t) The resulting Liouville operator is 1 [ - i V s d W + N-2 ( c t ) 2 d2 + AT W d t d 2 - C t 2 c 2 (3.13) £ = -2/ic+c + 2NTS (3.14) Now we can make the substitution d —• N, and d f -» 1. This leaves us with 1 £ = -2/xctc + 2NT„ -NsJc2 + (ct) + iVsct - c t 2 c 2 (3.15) At this point, we could take N large and identify c and ct with Fokker-Planck operators and write a non-linear Fokker-Planck equation for the variable z = 2x — 1. We will skip this for now and make the (canonical) transformation to a and a, c + 1 a = a = 2ct. (We are re-using the notation; a here is actually N times the original a). The Liouville operator now reads 1 £ = -jio) (2a - 1) + sr^da (1 - a) 2NT 6 da (1 - a) (3.16) Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 57 The system size expansion is obtained by taking a and a to represent Fokker-Planck operators (see section 2.4.3) rather than (scaled) occupation number operators, i.e. x, P-In what we have called the ^-representation (see section 2.4.3), the operator x is multi-plication by x and the operator — p is partial differentiation with respect to x. In this manner we arrive at the same Fokker-Planck equation derived by the classical method and presented in Eq. (3.7). 3.3 Heterogeneous diffusion in the Wright-Fisher and Moran models. We now review the behaviour of the Wright-Fisher or Moran model as described by the associated Fokker-Planck equation. First, we consider the model in the absence of selective advantage or mutation. The Fokker-Planck equation for the allele probability density W (x, t) , in the limit of pure genetic drift is dW(x,t)_ 1 & dt -2Ndx*X{1 X h { S - U ) were time is in units of the generation time r g . One notices the effect of random sampling of gametes is to produce a singular, het-erogeneous diffusion term—it vanishes at the boundaries x = 0 and x = 1. This results in a dynamical buildup of probability at the absorbing states x = 0 and x = 1, and correspondingly, an asymptotic decay of the probability density everywhere else. This arises not only because of the absorbing states, it is also due to the state-dependent or heterogeneous diffusion (see section 2.2 of chapter 2) given by the diffusion term B (x) = x (1 — x) /2N. This term generates an effective force that drives the system toward the boundaries. In other words, the behaviour for vanishing s and p is more Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 58 a) "replacement" with selection b)m a « » A P(At)= exp(-uAt )/u Figure 3.2: Stochastic Processes in the Moran Model. The arrows represent random events, a) Replacement. This event occurs at time intervals At given by the distribu-tion P (At) oc exp (At/r) with r defined in the figure. The upper path is taken with probability p = (1 + s) /2, while the lower path is taken with probability p = (1 — s) /2. b) Reversible mutation at rate /x. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 59 complicated than a random walk with "cliffs" at the boundary. We shall see that one can associate a potential with this force that is symmetric about x = 1/2. We can make this symmetry manifest by transforming to the variable z = 2x — 1, z G [— 1,1]. In the preceding section I used a operator version of this transformation and will use it again in the next chapter when I discuss the Stepping-Stone model. In terms of this new variable, the Fokker-Planck operator becomes 1 d2 which is obviously symmetric with respect to z —• —z. When mutation and selection vanish, this symmetry is "spontaneously broken", meaning that the asymptotic proba-bility distribution is W(Z) = 6(Z-1),OTW (z) =6{z + l). Reversible mutation A A a and selective advantage in the sampling of gametes, as discussed in the preceding section, give rise to additional "drift" terms in the Fokker-Planck equation. If /J, is non-zero, the absorbing states vanish. (Recall from section 2.2.2, that both the drift term and the diffusion term must vanish at the same point to produce an absorbing state.) A non-zero coefficient of selection, s, breaks the z —*- — z symmetry and makes one absorbing state preferable. The steady-state distribution and the dynamics of the expectation and variance of o the (discrete) Wright-Fisher model had been known since the work of Fisher and Wright. However, the full solution of the FP equation (??), in terms of eigenvalues and eigenfunc-tions for various regimes of selection and mutation was presented in a series of papers by Kimura [51] [5] [13]. Kimura gave a rather complete analysis, not only of Eq. (??), but also of the conjugate or backward equation [52]. It will not be necessary for us to review these details. Instead, we will look at this problem from another point of view, Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 60 that of the "effective potential". For the moment, however, we will make a temporary digression and completely ignore the fluctuations to see how this system behaves using a rather naive type of mean field theory. 3.4 Mean-field theory. Mathematical modeling of population systems is dominated by the dynamical systems approach [53] . This requires the underlying assumption that in the limit of large popu-lation size, fluctuations are irrelevant, and average out in such a way that the behaviour of the system is given by dynamical equations for the densities of the various species comprising the population (or sub-populations in spatial models). This amounts to a kind of 'dynamical mean-field' ansatz, and is the foundation of the "Fisherian" school of neo-Darwinian mass selection. Later we will see precisely how this assumption fails but first it will be instructive to proceed with this assumption. Using the Wright-Fisher or Moran model defined above, we write an equation for the expectation value of x, which we can simply and without confusion denote by ix{ty in this section. The dynamical equation reads where time is measured in units of the generation time Tg. Alternatively, we can define a potential (or Lyaponov function) dx (t) p. (2x (t) -l) + sx (t) {l-x (t)) (3.18) dt (3.19) and equate the time derivative of x, to the negative derivative of V, dx ~dt dv dx Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 61 The steady state value of x (t) will be given by the minimum of V. In fact, Eq. (3.18) is quadratic so its explicit time-dependent solution is easily found. Figure 3.3 shows the potential for several values of p/s. Note that the potential varies quadratically about its minimum. Setting the right side of Eq. (3.18) to zero, one finds the steady state value of z, where the positive (negative) sign is chosen if p/s > 0 (p/s < 0). The steady state arises from a balance between the forces of mutation and selection and is known as mutation-selection balance in the population genetics literature [47]. We would also like to know about the stability of system based upon Eq. (3.18). This is given by the sign of the partial derivative that results from acting on the right-hand side of Eq. (3.18) with d/dx, and evaluating the result for x — x. Its value is — \/(s2 + 4/Lt2). So we see that the system is stable as long as either s or p is non-zero. In much of what follows, mean-field theory will prove to be inadequate, yet this result will still hold, at least qualitatively in the thermodynamic limit, even in the presence of fluctuations—both in the Wright-Fisher model, as well as the Stepping-Stone model. This is somewhat surprising since it is generally not the case in stochastic population 3.5 The Effective Potential. To go beyond the simple mean field picture presented in the previous section, I make use of an "effective potential" for our Fokker-Planck equation. For multivariate systems, such a potential does not, in general, exist. However, for any single variable Fokker-Planck equation and in particular, the one of Eq. (3.7), one can construct an "effective potential" [16]. models. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 62 Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 63 The steady-state distribution, Ws (x) of the stochastic process defined by equations 3.7 is easily found since for any Fokker-Planck equation of this form, one has W s { x ) = B i h ) e x p B(y) where C is some normalization constant. This allows us to define a potential U (x), such that Ws(x) = Cexp{-U(x)). If B (x) were constant, which is the case for homogeneous diffusion, for example, in models where Gaussian noise is added to the mean field dynamics, this potential would be identical to that of the deterministic differential equation for the expectation value. In the present case, A and B are given by A(x) = -fi{2x-l) + sx(l-x) (3.20) B{x) = N^xil-x), and where I have taken the units of time to be r g . B (x) depends upon x, so the effective potential becomes U (x) = (1 - 2Nfi) In (x (1 - x)) - 2Nsx. (3.21) This gives us the steady-state distribution Ws (x) = Ce2Nsx [x (1 - x ) ] 2 ^ 1 . The normalization constant can be found from the intergral C " 1 = f1 e2Nsx[x(l-x)]2N^1dx (3.22) Jo r2 (2N(M) r (4Nfj.) •F{2Nn,4Nfj,-2Ns) Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 64 where F (a, b; x) is the hypergeometric function. For s = 0, Ws (x) is simply a symmet-ric Beta distribution, however, in the genetics literature, Ws (x) is known as "Wright's distribution". In contrast to the mean-field potential V (x), the effective potential U (x) in Eq. (3.21) is strongly anharmonic. More importantly, however, the effective potential changes sign when 2Np = 1, and therefore is qualitatively different than its mean-field counterpart. Figure (3.4) shows a comparison of U (x) and V (x) for p/s ~ 1. An important point to make here is that our new, effective potential tells us not only about the steady-state distribution, it also tells us something about the dynamics of the system. The point has been somewhat overlooked in the population genetics literature. In Figure (3.5) I show two simulations of the Wright-Fisher model for the two regimes, 2Np » 1 and 2Np « 1, with s = 0. These simulations were produced by integrating the stochastic differential equations (s.d.e.'s) associated with the Fokker-Planck Equation (3.7). For these simulations, the population size was N = 1000 and the mutation rate was adjusted so that 2Np = 20 (top) and 2Np = 10~4 (bottom). The method of integration is discussed in section 1 of appendix B. Let us consider what happens if we start with s — 0, and "tune" the mutation rate slowly to zero. First, as p becomes less that the effective potential flips. As p is further reduced, we observe the behaviour shown in the bottom plot of Figure (3.5). In the limit p 0, the symmetry of U (x) is broken as the system is trapped in one of the absorbing states. The manner in which this occurs, however, is quite different from what is common in equilibrium statistical mechanics. Figure (3.6) shows the effective potential for 2Np slightly larger than 1 (lower plot) and then slightly larger than 1 (upper plot). What one should notice, is that at this "critical point", the order parameter (here the expectation value of x) changes discontinuously. But unlike first order transitions in equilibrium statistical mechanics, where the order parameter also Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 65 changes discontinuously, there is also critical behaviour near the transition. In the next chapter, we will find that this unusual behavior has important ramifi-cations for the scaling laws that relate critical exponents. In particular, the standard "hyperscaling" relation will have to be discarded and replaced with one that takes this effect into consideration. 3.6 The Broken Symmetry State: Fate of new mutants. Let us suppose that the mutation rate is effectively zero and enough time has passed for the system to reach an absorbing state. We would like to know the response of the system to a small perturbation caused by a single new mutation to a completely new allele. One might imagine an individual of the population being hit by a cosmic ray, or subject to some other rare process that produces this new allele. Furthermore, the new allele could be selectively advantageous, disadvantageous, or just plain neutral. The ensuing dynamics of the system will be to produce a "cascade", or "avalanche" of the new types in the population as a response to the new mutation. This is truly the Darwinian picture of the basic process of evolution, so it is quite desirable to know what can happen in this situation. One finds that even for an advantageous new allele, there is still a large probability that the resulting avalanche will die out. When the selective advantage is small, the fixation probability remains small regardless of the population size. Kimura has remarked that this fact has been never been truly appreciated by the mass selectionists, who imagine that every newly arisen beneficial allele becomes incorporated into the population [8]. Conversely, it is quite possible that the new allele will take over the population and reach fixation. For a neutral allele, the fixation probability does depend on the population size. It is arguable that the most important contribution by Kimura to Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 66 2 p N » l I i I i L i I i I 0 0.2 0.4 0.6 0.8 x Figure 3.4: Comparison of the stochastic "effective potential" U (x) with the mean-field potential V (x) for different values of 2fiN. Top: 2/J.N 3> 1, the minimum is shifted by an amount ~ AT - 1 . Center: 2fiN. ^ 1. Bottom: 2fiN <C 1, the effective potential U (x) has turned over. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 67 Figure 3.5: Simulations of the Wright-Fisher model using stochastic differential equations. For 2pN 1 (top), the fluctuations are approximately Gaussian around x = 1/2. For 2pN <C 1 (bottom), the system stays remains near the absorbing states. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 68 2u.N>l U(x) L U(x) L i 1 1 i I i I i l 0 0.2 0.4 0.6 0.8 1 x Figure 3.6: Discontinuous behaviour of the order parameter at the critical point of the Wright-Fisher model. The effective potential U (x) changes sign at 2pN = 1. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 69 theoretical population genetics was to propose that most, if not all genetic variation observed in Nature is due to such neutral mutations [8][9]. This postulate, know simply as The Neutral Theory provides one of the few truly quantitative theories in Biology. It should not surprise us then, that as a predecessor to his Neutral Theory, Kimura has contributed the most to the simple problem of fixation probabilities of new muntant alleles discussed in this section [13] [52]. The Fokker-Planck description we have been using so far is not well-suited for the type of calculation needed here, for what we want to calculate is effectively a first-passage type of problem. Using the Fokker-Planck equation to this end would require knowledge of the eigenvalues and eigenfunctions of the equation, which is indeed possible to derive [5]. Proceeding in this manner is much too cumbersome. The more general and alos aesthetic solution is the use the backward, or adjoint equation introduced in section 2.2.3. —u(x0,t) = C~pP(x0)u(x0,t). (3.23) This equation describes the fixation probability at time t, given that the system was initially at XQ. This is precisely how Kimura proceeded to his solution [52]. He solved Eq. (3.23) when the left hand side is set equal to zero. This produces the asymptotic fixation or survival probability, — lim^oo u (XQ, t). Kimura found I _ p-4Nsxo P ~ = , e . m s • (3-24) I will not review that calculation here. Instead a more illuminating exercise will be to find a first order approximation, valid for very large N, that gives us the desired results. I will show that in the thermodynamic, or large system size limit, the fate of a new mutant allele is described by a simple branching process. I will also show that when the new allele is selectively neutral, the branching process is critical. This limiting result is important to understand since many stochastic process have rare event regimes that converge to Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 70 branching processes for large system sizes. Indeed, I show in chapter 4 that the critical branching process can be understood within the framework of the renormalization group as a fixed point that describes the critical behavior of many stochastic models above their upper critical dimension. To show how this regime of the model emulates a branching process, we start with the master equation for Moran's model given in Eq. (3.3) and set ji = 0. Then we use the conservation law, N = na + nA, to eliminate nA and simply use n for na. We are left with a master equation describing only the A's. For clarity, and since we are not interested in a system size expansion, we do not use the shift operators introduced in Eq. (3.25) here. The master equation reads dPn l + s dt 2N 1 - s [(N-n+ 1) {n - 1) P n _! — (N — n)nPn] [(N-n + 1) (n + 1) Pn+1 -(N-n) nPn). (3.25) 2N (For ease of notation, I have dropped the explicit time-dependence of Pn.) For very large N, we can ignore terms of O (N'1) in the master equation, which simplifies to dPn l + s dt 2 1 - s [{n - 1) Pn-! - nPn] (3.26) [{n+l)Pn+1-nPn}. We recognize Eq. (3.26) as a master equation that describes a branching process with branch rate 7 = - ^ - , (3-27) and decay rate giving us the growth rate A = — - . (3.28) A = 7 - A = s. Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 71 I should mention that this result (3.26) can be very easily derived using the operator form of Moran's model (3.16). Branching processes were discussed as an example in section 2.4.2 where we illustrated operator techniques. A n alternative derivation using classical methods is presented in appendix C. We can use either results, Eq. (C.8) or Eq. (2.62). The time dependence of the survival probability is given by 1 + s — (1 — s)e ^l Its asymptotic value for s > 0, the "fixation" probability is, A 2s P ^ = lim Ps (t) = - = — — , (3.30) t-*oo 7 1 + S which for small s tells us Poo * 2s. (3.31) This agrees with a very early calculation by Haldane [54], and the result due to Kimura, Eq.(3.24) with XQ = 1/2AT (since Kimura uses the genetics convention of taking the population size to be 2N), P ~ = 7 ^ ~ 7 - 2 S , (3.32) for small s and Af —* oo. If we repeat the above argument for a selectively disadvantageous new allele A, in other words, we let s —> —s, we find from Eq. (C.12) the size distribution, Qn describing the total number n of new alleles produced in the avalanche. For large n, Qn ~ J - n - 3 / V ^ (3.33) with the characteristic size Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 72 When s = 0, the process is critical, £ diverges and both Qn and Ps (t) follow power laws Pn ~ n" 3 / 2 , (3.35) and 3.7 Discussion. Taking into account the fluctuations, we have found that the effective potential is not the same as the mean-field potential or "fitness landscape" as it sometimes is called in the literature only concerned with a mean-field dynamical systems description Even when 2N/J, » 1, minimization of the effective potential reveals a shift of O (-/V-1) in the steady state value of x, For 2N/J, < 1, the situation is qualitatively different as the stability is no longer given by the mean-field result. In discussions of evolutionary models, there is frequent reference to "fitness" landscapes or "selective surfaces". These approaches typically ignore the full implications of the particular type of state-dependent noise that arises from a microscopic description, and as we have shown are only valid in the limit 2N/i » 1. In fact, most models involving selective surfaces ignore mutation completely, making the concept of a fitness landscape quite irrelevant. This observation is the basis of Kimura's neutral theory that postulates a model of evolution at the molecular level that is very different than the sort of mass selection envisioned by Fisher and other followers of the determinist school of thought. We have re-examined the Wright-Fisher/Moran model from the point of view of the 'effective potential' and found that in the more realistic limit of small mutation rate, it Chapter 3. Spontaneous Symmetry Breaking in the Wright-Fisher/Moran Model 73 is this effective potential that is the relevant landscape upon which the system evolves. More germane to this thesis, however, the question arises as to whether or not the broken symmetry and critical behavior exhibited in the Wright-Fisher/Moran model persists when this model is generalized to include spatial structure, or in physics terminology, spatial degrees of freedom. Chapter 4 The Stepping-Stone Model Revisited 4.1 Introduction. Wright was the first to quantitatively examine the effects of random genetic drift for geographically distributed populations with his 'Island' model. According to this model, in which time is discrete, symmetric migration takes place between an infinite number of subpopulations or colonies, which otherwise behave according to the Wright-Fisher model discussed in chapter 2. As a generalization of the Island model, Kimura and Weiss [1] introduced the so-called "Stepping-Stone" model, which in addition to long-range migration, incorporated nearest-neighbour migration between the colonies, now arranged on a lattice. By solving the steady-state correlation function of this model, they were able to derive a number of interesting results concerning genetic differentiation among sub-populations. Unfortunately, the approach described in the previous chapter of finding an effective potential does not help us much for the Stepping-Stone model. Such a potential can actually be derived, but it does not provide a convenient way to calculate the desired results. We must resort to other means, and will be able to infer the behavior of this potential with respect to symmetry breaking and critical properties. We will consider two slightly different models in this chapter, both of which can be thought of as lattice generalizations of Moran's model. We find that these models are equivalent to the original Stepping-Stone model of Kimura and Weiss [1] in the same sense 74 Chapter 4. The Stepping-Stone Model Revisited 75 as the equivalence of Moran's model and the Wright-Fisher model discussed in the pre-vious chapter. In this chapter, I formulate these models using the creation-annihilation operator formalism described in chapter 2. First, I reproduce results reminiscent of those from Kimura and Weiss [1] [2] and Malecot[55] for the steady-state two-point correlation function. I then extend their work by solving the dynamics of the two-point correla-tion function for the neutral theory and also for small, but finite selection in the broken symmetry state. An alternative method for finding this correlation function, which is closer to the original work by Kimura and Weiss, would be to use the response function formalism or the operator form of the Fokker-Planck equation introduced in chapter 2. In section 5 of this chapter, I will give these alternative formulations of the Stepping-Stone model and will find that they are indeed equivalent to the occupation number description. Selection was absent in the Stepping-Stone models considered by both Kimura and Weiss[l][2] and Malecot[55]. Several authors, however, have discussed the Stepping-Stone model with selection and provided theorems concerning the ergodic properties of the model for various cases of mutation and selection. Both Itatsu[56], and Shiga and Uchiyama[57] consider asymmetric mutation rates, that is, A —> a with rate // and a —»• A at rate v, and prove that a heterozygous steady-state exists when /J, > 0, v > 0 and s is arbitrary. If one of the mutation rates is identically 0 and the other strictly positive, then a heterozygous steady-state is found only for s greater than some critical value sc. While Shiga [57] does not rule out the possibility that sc = 0, Itatsu finds that j^- < sc < 1. For the mean-field dynamical system, it is easy to show that for this case, sc = p.1 In section 4.6, I consider the case where p = 0, v = 0 and s ^ 0. For this case both Itatsu, and Shiga and Uchiyama find that the only stable steady state is the fixed state of the selectively advantageous allele. This finding is confirmed by the analysis in section 4.6. *It should be noted that Itatsu parametizes selection in slightly different manner. His s is half of that used in this thesis. Chapter 4. The Stepping-Stone Model Revisited 76 In addition to this, the asymptotic dynamics of the approach to this steady-state. The two models considered in this thesis have the same on-site processes, but differ in the way neighboring sites are coupled. In the first model, which I will refer to as the 'Exchange Coupled Stepping-Stone (ECSS) Model', the number of individuals at each site n, is fixed and independent of position in the lattice. Nearest neighbor sites are coupled by a random exchange of a pair of individuals, one from each site. Since all processes, both on-site and nearest-neighbor, conserve the number of individuals, we have only one independent field. The second model on the other hand, I refer to as the 'Diffusively Coupled Stepping-Stone (DCSS) Model', uses random diffusion or dispersal of individuals to neighbouring sites, does not conserve the number of individuals at each site and therefore has two independent fields. In the DCSS model for a homogeneously distributed population, however, the average number at each site is constant. I will show that the field representing the total number is independent of the gene-frequency field and therefore described by diffusive fluctuations. First, I consider the on-site processes which are identical in both models. In the following, a-operators describe a-type individuals, while 6-operators will describe A-type individuals. The derivation of the Stepping-Stone model here is a generalization of the operator form of Moran's model presented in chapter 3. Symmetric mutation is the same as in Moran's model and is described by one-body operators Anut = fi £ {ath - 6j6P + blar - arar) . (4.1) r Genetic drift along with a selective advantage is given by two-body operators £Sei/drift = 7^ £ ((1 - s) alalarbT + (1 + s) &r&PaP6P - 2ar6raP6P) . (4.2) r Without coupling, we would have N (the number of sites) independent populations of size n, evolving under Moran's dynamics. Chapter 4. The Stepping-Stone Model Revisited 77 4.1.1 The Exchange Coupled Stepping-Stone Model (ECSS). Now the total Liouville operator is given by £ E C S S = 4x + Cmut + £ S e l / d r i f t , (4.3) with the exchange processes described by £ e x , given by ^ex = ^ E E ( & i a i + e « r & r + e ~ <4&t + ea r0 r+e) r e + ^ E E ( a i 6 i + e & r a r + e - 6 Pa r + 06 ra P + e) , r e where e is a lattice vector of magnitude I, and the lattices can be chosen to be e.g the linear, square or cubic lattice in d = 1, 2,3, respectively. Here m is the migration rate. Usually in the physics literature, m is used as an index. We break with that tradition here and bend to the convention used in population genetics. In the continuum limit, however, we denote by D, the diffusion constant, which is given by D = ml2. 4.1.2 The Diffusively Coupled Stepping-Stone Model (DCSS). For this model, we have A D C S S = £ d i f f + £ m u t + ^sel/drift, the only difference from the ECSS model is the random walk term Cdis, given by Aiiff = m E E ( f l r + e a r + 4 - e a r ~ 2a^ar j (4.4) r e r e which describes "hopping" to neighboring sites. Once again, in both models, D — ml2 appears as the diffusion constant in the continuum limit. Chapter 4. The Stepping-Stone Model Revisited 78 4.2 The Neutral Regime. Under vanishing selection, no allele is favoured and as was discussed in the previous chapter, if the mutation rate also vanishes, the system becomes critical. Moreover, even in the spatial model, the critical point remains at /x = 0. This fact will become clear from inspection of the total Liouville operator—none of the interacting terms in £sei/drift produce corrections to the single particle propagator (the Green's function). We will also be able to show this by examining the Heisenberg equations for the densities. These do not couple to higher correlations and are dynamically stable for non-zero ji. In this section we will obtain exact results for the two-point correlation function for both finite and vanishing mutation rate. We will also find that the critical behavior of the d-dimensional system for fl = 0 is intimately tied to return probability of a random walk in d dimensions. Results similar to those of this section were found by Sawyer[3] and Nagalaki[4] for an infinite alleles model and are discussed below. First, we make a canonical transformation to sum and difference operators of the form used for the Moran model in chapter 3, This is analogous to the transformation z = 2x — 1, of the Wright-Fisher/Moran Fokker-Planck equation mentioned in the last chapter, but carried out at every site. Once again we stress the necessity of using a canonical transformation so that we preserve the commutation relations, i.e. so that [c, c^ j = 1 and [rf,^] = 1. For the on-site terms, this results in c = (a — b) jn c t = n ( a t _ 6 t ) dt = I ( a t + 6 t ) d = a + b (4.5) (4.6) Chapter 4. The Stepping-Stone Model Revisited 79 and Adrift = (2n)_1 (n-2cl4drdr - cPcPcpCp) (4.7) r In the ECSS Model, the conservation law at each site allows one to replace the dr oper-ators with the real number n, and the d) operators with 1. Including the exchange and diffusion parts, we now find respectively, ^ECSS = m £ XI (Cr+eCr + 4-eCr ~ 2crCr) (4.8) r e r + (2n)_1 £ (44 - 44CrCr) r and £ D C s s = m5Z£ (4+eCr + cP_ecP - 24cr) (4.9) r e +mJ2 ( 4+edr + 4-edr ~ Hdv) r e -2^X]cPcP r +7^ £ (n"2c rc rd rd r - 44crcr) r If we now look at the Heisenberg equation for the expectation of the difference field, we find j t (cr) = -2// (cr), (4.10) which tells us that the critical point remains at /i = 0. 4.2.1 Fourier transformation to wave-vector space. We will assume periodic boundary conditions. For translationally invariant systems, further analysis is accomplished and greatly simplified by transforming to wave-vector space. Chapter 4. The Stepping-Stone Model Revisited We use Fourier transforms defined by 80 cr = iV" 1 / 2 E c k e i k r , 4 = N-1'2 J2c{e~ikr (4.11) k k k k so that ck = i V - 1 / 2 E c r e - i k r 4 = N-V2Y4eikr (4.12) r r dk = N-1'2 £d r e~ i k r d[ = N~l>2 J24e*,r r r In the DCSS model, we have diffusion terms for both fields CD = -m^w k 4c k - m^]^w k 44 , (4-13) k k 5 where w * = ( 2 _ e ' k e _ e " i k e ) (4-1 4) e On hypercubic lattices in d dimensions, a ; k = 4 2 j s i n 2 ^ - , (4.15) 3=1 where I is the lattice constant. For small kl, o>k ~ k 2£ 2 . In the ECSS model, only the first term in Eq. (4.13) is present. The mutation term is the same for both ECSS and DCSS models, Anut = -2/nEckck- ( 4 - 1 6 ) k The selection/drift term for the ECSS model is £d"ft = ^ E C k c t - k - ^ E C k j ^ C k l + q ^ - q , (4-17) k ki,k 2 ,q Chapter 4. The Stepping-Stone Model Revisited 81 and for the DCSS model C d r i & = 2 n N ^ ( ^ ^ k l ^ ^ i + q ^ - q - C ^ C ^ C ^ + q C ^ . q ) . (4.18) k i , k 2 , q 4.2.2 Equation of motion for the two-point correlation function. I now proceed to solve the equation of motion for the Fourier transform of the two-point correlation function of the difference field, which is defined by the expectation G u (k,t) = (ckc_k) = (S |ckc_k| V) • (4.19) Here I use the Heisenberg representation, so that ^ (ckc-k) = [ckc_k, £] Taking expectations, one has ^GH,(k,t) = ( [c k c_ k ,£]} (4.20) The right side of Eq. (4.20) is found by direct evaluation of the commutator, using the canonical commutation relations. Examples of this procedure for creation and annihila-tion operators were given in chapter 2. A vast simplification occurs upon taking the expectation. This is due to the identity S ip) = 0 (4.21) I do not present the somewhat tedious algebra here. The resulting equation of motion for the ECSS model is ~GU (k,t) = -2 (mwk + 2(i) GH (k,t) +n~1-~Yl G4>4> (q,*) • (4-22) Chapter 4. The Stepping-Stone Model Revisited 82 For the diffusively-coupled model, we have the system of equations jGH{^t) = -2 {mwk + 2//) G u (k,t) (4.23) q q ^ t f ^ (k,t) = -2mwkHH (k,t), (4.24) where I have denoted by H^ (k,£), the Fourier transformed two-point correlation function of the total number field, Hu (k,t) = (dkd_k) = (S |dkd_k|ip). representing the Fourier transform of the total number field. For the case of a homoge-neously distributed population, one has f iVrc2,q=0 Hu (q,*) = < { 0 ,q^0 and therefore one recovers the same equation for G^p (k,t) as in the exchange-coupled model. For the rest of this thesis, we will assume a homogeneously distributed population and work only with the difference fields. This means that as far as we're concerned, the DCSS and the ECSS are equivalent. It would be interesting, however, to use the system of equations, Eqs. (4.23), to examine situations where the population is initially localized over some small region. We note that unlike many other stochastic population models, the equation of mo-tion for the two-point correlation function does not couple to higher order correlation functions. This fact allows us to find exact results for this model. Chapter 4. The Stepping-Stone Model Revisited 83 4.2.3 The steady-state. Using Eq. (4.22), we find the steady-state solution by trying the solution GW(k) = ^ — , (4.25) mu k + 2p where C is an undetermined constant. This results in an equation for C, c = l-SLy—I— (4.26) n nN ^—' mu)a + 2u q M r = ( l - C / ) / n , where / = i V - i y 1 . mwq + 2/i with the solution for C, c 1 n + I We recognize the two-point correlation function, in real space and evaluated at r = 0, as the sum Gu (*) = jj £ G** (q) and using the solution for C, N q n + I' To get the correlation function in real space, we need to perform the inverse Fourier transform. We wish to study the spatial dependence of the correlation function in the continuum limit in which details such as lattice structure do not matter. In the continuum limit a>k = k2l2, so mwq = mk2l2 = Dk2, where we have obtained the diffusion constant D = ml2. The two-point correlation function becomes Chapter 4. The Stepping-Stone Model Revisited 84 which is C times a very familiar quantity in statistical physics: the two-point correlator of the Gaussian model. Furthermore, in the continuum limit the inverse Fourier transform becomes an integral k (27T) Jo where A I"1 is a cut-off parameter. We define the correlation length and the dimensionless constant K = n" 1 / = l—: j ddk —^ and end up with the result G** (r) = rr^ exp(-^ )/iri d - 1 2 The integral for K is easily evaluated (n2D(i) -1 /2 d = 1 « = < 2TTDT), ' 2ir2Dn > d = 3 (4.28) (4.29) (4.30) (4.31) All along, we have be measuring time in generations. We can explicitly put back in the generation rg if needed, f (n2DfJLTg)-1/2 ,d=l «= < . d = 2 . (4.32) MM) H _ 9 2-KTgDn ; " — ^ A 2ir2TaDn , d = 3 I will refer to K as the "strength of local differentiation", for as K becomes large, the amplitude of correlation function becomes one and the local populations are fixed to one Chapter 4. The Stepping-Stone Model Revisited 85 allele or the other. When K is small, the local populations are polymorphic. As an alternative to K, a more appropriate measure, would be the quantity ^q -^ In population genetics, this quantity referred to as FST, where FST stands for fixation index. If we go back to the lattice description and use our results for K to calculate FST, we find 1+nWmfj. d = 1 FST = \ 1 - ^ / 1 , ( 3 ) ' d = 2 (4-33) 1 _ A — Q l+ 2 7 T 2 m n ) " — ° These results for FST should be compared to the steady-state value of FST in Wright's Island model [58], FST = —\ . (4.34) 1 + 4nm v ' The Island model of Wright also consists of colonies of size n, but in contrast to the Stepping-Stone model, each colony exchanges migrants with all the rest. We see that our results for FST f ° r the Stepping-Stone model only qualitatively agree with the Island model for d — 3 ( 4nm vs. 47r2nm, as Wright uses the genetics convention of 2N for the population size). For d < 2, FST for the Stepping-Stone model depends explicitly on the mutation rate and is therefore qualitatively different that its cousin in the Island model. Figure 4.1 shows several snapshots from computer simulations of the d = 2, Stepping-Stone model with N = 200 X 200 sites, and illustrates the effects of varying £ and n. Simulations in the same row have the same correlation length, while those in the same column have the same value of FST- The method of simulation is described in appendix B. The results up to this point, while perhaps expressed somewhat differently, are similar to those of Kimura and Weiss [1][2]. We now leave them behind. Chapter 4. The Stepping-Stone Model Revisited 86 4.3 Broken Symmetry or d < 2. For any finite size, or for spatial dimension d < 2, I diverges as p —> 0, and hence This implies fixation and the spontaneously broken symmetry likewise observed in the Wright-Fisher/Moran model. Another way to see this is to look at $ k , for k —> 0, which is found to be G u (k) = 6 (k). Perhaps the simplest way to see this, however, is to look at the results for K in the continuum limit. For d < 2, K, diverges as fi —* 0, and the fixation index FST 1-This is not the case for d > 2 as K. is independent of fj,. This result combined with the diverging correlation length £, means that no matter how large we make the system, it will eventually reach fixation. 4.4 Critical dynamics of the two-point correlation function. 4.4.1 The decay of heterozygosity. To go beyond the steady-state solution and solve for the dynamics, I return to the equation of motion 4.22, which we re-express as j G H (k, t) = -2 {Duk + 2/x) GH (k, t) + n^h (t), (4.35) with h(t) = l - $ r = 0 = l - F s r ( t ) . Chapter 4. The Stepping-Stone Model Revisited 87 H Figure 4.1: Snapshots of computer simulations of the d = 2 Stepping-Stone model on the 200 X 200 site square lattice. The mutation rate and migration rates were adjusted so that the simulations in each row have the same correlation length, while the simulations in the same column have the same value of FST-Chapter 4. The Stepping-Stone Model Revisited 88 The quantity h (t) is known in the genetics literature as the heterozygosity or more prop-erly as the local heterozygosity. From above, we note that for [i = 0, and d < 2, h (t) decays asymptotically to 0. The integral solution to 4.35 is given by rt GH (k,t) = - I dt'e-^-^h {t'). n Jo Here we have used the initial condition G^ (k,0) = 0. This means that we start out with n a = riA = n/2 at each site. Summing over both sides of Eq. (4.35) results in an equation for h (t) h(t) = 1 - ^ - ^ 2 f dt'e-^-^h(f), (4.36) k J Changing the order of summation and integration, we recognize the sum £ e-2m<t-t'>Wk = G2m (0,t - t') (4.37) k as the diffusion Green's function for a random walk in d dimensions, with initial and final position at the origin, and diffusion rate 2m (Appendix A). We have i r* h(t) = l-- dt'G2m (0,t - f) h {t1), (4.38) n Jo with G2m (0,t)for the hypercubic lattices given by G2m (0,t) = e~imdt [Jo (4mt)]d , where lo (4mt)is a modified Bessel function of the first kind. Taking the Laplace trans-form of both sides results in the equation ~h is) = s'1- ~ * w • (4.39) K 1 l + G2m (0,s) /n V ' Using the scale transformation properties of the Laplace transform, Gbm (0,s) = b^Gm (0,6s), Chapter 4. The Stepping-Stone Model Revisited 89 and taking b = 4dm, we can write G2m (0,s) = (Mm)'1 G i . (0,4dms), 2d where Gj_ is the Laplace transform of a random walk with rate ^ . As shown in appendix A, the s = 0 limit of Gj_ is fundamental and universal constant that is related to the 2d return probability Q (0) of a random walk through « ( 0 ) = 1 - < 5 7 W ( 4 ' 4 0 ) 2d v ' For d < 2, lims^oGj_ (0,s) = oo, so as noted above G$& (r = 0) —• 1 at large times. For 2d d > 2, however, (r = 0) approaches a constant. We can summarize the situation for arbitrary dimension by writing lim G u (r = 0) = \ — . (4.41) To find the dynamics of the approach of the two-point correlation function evaluated at r = 0 to its limiting value, we take the Laplace transform of the Green's function for random walks which is explicitly calculated for the linear [d = 1), square [d = 2) and cubic (d = 3) lattices in appendix A. Asymptotically, the Laplace Transform for small s has the dependence (8ms)~ 1 / 2,d= 1 G2m (0,s) ~ I Hsm/s) d 2 / 4 42) 16m7r ^ V 1 \ 12m + B V 1 2 3 m 3 ' d ~ 3 where the constants C and B for the d = 3 result are found from a simple relation to the universal return probability of a random walk. Using this result gives us the behavior of h (s) for small s, 2s ' " — L h{8)~{ 4m7r ns]n(8jn/s) ,d = 2 • (4.43) Chapter 4. The Stepping-Stone Model Revisited 90 Standard 'Tauberian' theorems [59] for Laplace transforms relate the behaviour of h(s) for small s to the behaviour of h (t) for large t. One finds r ,d = 1 n y m t h(t)~ 1 nymt' d = 2 • (4.44) l n ( m n 2 t / 4 ) 3mn j q V Zmn+C ) " — ° Furthermore, even though the situation is qualitatively different for d = 3, in that h (t) approaches some non-zero constant, it is interesting to determine how h (t) approaches this constant. From the result above, if the constant is subtracted, one finds that the small s behaviour of the Laplace transform of g (t) = h (t) — lim^oo h (t) is given by 9(s) mn2 4s ' and hence 9{t)~\[^t The integral equation for h (t), Eq. (4.36), is readily solved by numerical quadrature. This solution is presented in Figure (4.2), confirming our results for the asymptotic decay of h(t). Several remarks are in order: 1) the logarithmic dependence of this quantity for d = 2 is indicative of what is termed "upper critical dimension" in statistical physics; 2) Another interesting feature is the "crossover" region at early times; but most importantly, 3) even though strictly speaking, and as we have argued in the preceding section, fixation occurs for d = 2, the extremely slow logarithmic dependence makes one wonder whether or not, in any practical sense, the system ever gets there. In summary we have found that the asymptotic value of (0,t), known as FST M the genetics literature, is related to the return probability of a random walk through 1 lim FST (t) = : «->«> 5 J v ; l + dnm(l-Qd) Chapter 4. The Stepping-Stone Model Revisited 91 In the thermodynamic limit its approach to this value is ~ t 1 / ' 2 for d = 1, ~ 1/lnt for d = 2 and ~ r 1 / 2 for d = 3. 4.4.2 The infinite alleles model. Sawyer [3] and Nagalaki[4] both investigated a related neutral model introduced by Malecot [55] that differs from the present model by having infinitely many alleles or types. They found asymptotics similar to those of Eqn. (4.44) for d = 1 and 2, but for a related quantity. Sawyer's starting point is an equation derived by Malecot [55] that describes the time dependence of the probability of identity by descent. This quantity is denoted / ( r 1 ; r 2 ; t) and gives the probability that two randomly selected individuals located at r i and r 2 are of the same type due to one or more common ancestors. For infinitely many alleles and assuming that at t = 0 all types in the population are unique, as in done in the paper by Sawyer, the probability of identity by descent (i.b.d) equals the probability of identity by type (i.b.t.). A similar equation was studied by Nagalaki for the probability of identity by type. One can show that in a two-allele model, the probability of i.b.t. is related to the correlation function G^ (T\, r 2 ) . First we observe that when n is large, the discrete site variables become frequencies which we denote by the operators fT = n _ 1 a r . The probability of i.b.t. is given by the expectation / ( r i , r 2 ) = ( / r i / r 2 + ( l - / r i ) ( l - / r 2 ) ) , where the first (second) term is the contribution from selecting two a's (A's). The / -operators are related to the c-operators by the transformation Chapter 4. The Stepping-Stone Model Revisited 92 Time t Figure 4.2: Numerical solution for the heterozygosity h (t) = 1 — FST (t) showing the dynamics the fixation index (two-point correlation function). Note the slow logarithmic decay for d = 2. For d = 1, the asymptotic behaviour is clearly £ - 1 / 2 . Chapter 4. The Stepping-Stone Model Revisited 93 giving us / (r i ,r 2 ) = \ + \{cricr2) (4.45) = 2 + 2 ^ 0 0 ( r i ' r 2 ) • Note that J(ri,r 2 ) —» 1 iff G^(r i ,r 2 ) —»• 1. Furthermore we see that the initial conditions used to solve Eqn. (4.35) which were fr = 1/2 for all r, imply G^ (r 1 ; r 2; 0) = 0 and I (ri, r 2; 0) = 1/2, in agreement with Eqn. (4.45). Therefore, in the two-allele model studied in this thesis, / (ri, r 2; t) and G^ (r\, r 2; t) have the same asymptotics. It is quite interesting that I(r!,r 2;t) in the infinite-allele model studied by Sawyer has the same asymptotics as for the two-allele model. One is tempted to conjecture that / (r 1,r 2;t) has the same asymptotics in models with arbitrary number of alleles. An extension of the methods described in this thesis could likely answer this question. Sawyer also found that J(ri,r 2;t) for any two individuals in any bounded set was asymptotically equal to unity iff the random walk performed by the individuals was recurrent. This result is also very similar to the results of this thesis for the asymptotic value of FST a n d the existence of broken symmetry. Sawyer's result, however, is obtained by a very different method than that described in the previous subsection and pertains only to infinite alleles model. Further analysis of the infinite alleles model was carried out by Sawyer [60] who pro-vided a limit theorem for patch sizes and found that the expectation of N (t) , the average number of individuals of the same type scales as t1/2, t/\nt and t in d = 1, 2 and 3, re-spectively. 4.5 Alternative Formulations of the Stepping-Stone Model. In section 4.1 I formulated the Stepping-Stone model using the occupation number rep-resentation. Here I show that one can use either the response functional formalism or Chapter 4. The Stepping-Stone Model Revisited 94 the operator form of the Fokker-Planck equation to achieve an equivalent description. Moreover, the starting point for applying these techniques is the Langevin equation in-troduced by Kimura and Weiss [1] [2] in their original formulation of the Stepping-Stone model. I also point out that all of these reformulations of the Stepping-Stone model lead to the same path-integral and its associated action. I use this path-integral exclusively in the section 4.6 and chapter 5. To remind the reader, the Stepping-Stone model of Kimura and Weiss is spatial model where colonies (sub-populations) are arranged on a hypercubic lattice and are coupled by random migration of individuals. Apart from random migration between colonies, each sub-population evolves according to the Wright-Fisher model. We denote by fr, the frequency of the allele with a selective disadvantage, in the colony located at position r. The Langevin equation reads dfr (t) dt = m £ (/r+e + / r - e ~ 2/r) - y. (2/r - 1) - ST'1 fr ( l - / r ) (4.46) + J / r ( 1 . fr\r(t), (4.47) nr g where m is the migration rate, /i is the mutation rate, n is the total number at each site, and s is the coefficient of selection . The noise term is Gaussian and specified by the auto-correlation The amplitude of the noise in Eq. (4.46) is the result of Wright-Fisher binomial sampling at each site, and as was shown in section 3.2, is the same for Moran's dynamics. Finally, Eq. (4.46) is to be interpreted with Ito's prescription. Chapter 4. The Stepping-Stone Model Revisited 95 4.5.1 The Fokker-Planck equation and operators. The Fokker-Planck equivalent to Eq. (4.46) is found in the standard way, r where we have defined the lattice Laplacian A r , A r / r = J ] ( / r + e + / r - e - 2 / r ) . e For the operator description, we define the operator fT as multiplication by / r , and the operator gr, as the partial derivative y^-. The gives us the Liouville operator C = ygr[mAjr + 2pf\-l^+srg1Y/gJr(Kl-fr) (4.49) r r & r Equation (4.49) can also be obtained from £ECSS in Eq. (4.3) by performing a canonical transformation to eliminate the bT and 73r operators, followed by the system size expansion for Moran's model discussed in section 3.2.2, applied to each site. As discussed there and in chapter 2, this is accomplished by identifying appropriately scaled creation and annihilation operators with Fokker-Planck operators, which for the present discussion are the / r and gr operators. 4.5.2 The path-integral. One can easily proceed from Fokker-Planck equation to a path-integral description as discussed in chapter 2. The path-integral can also describe the lattice model, but here Chapter 4. The Stepping-Stone Model Revisited 96 we take the continuum limit.,Using the fields / r and fr, that correspond, respectively, to the operators fr and gr, the action is given by S[fr,fr] = ddrdt[frdtfr + fr(DV2 + 2pff-l)+sr-1frfr(l-fr)] (4.50) where for simplicity we have dropped the explicit time-dependence of the fields. The continuum limit gives us the diffusion constant D = ml2, where I is the lattice spacing. I use the path-integral description in section 4.6 below and for the renormalization group analysis in chapter 5. I should also point out that this path-integral is completely equivalent to that obtained using the coherent-state representation of the creation and annihilation operators. This procedure is also discussed in chapter 2. 4.5.3 The response functional method. We can use the response functional formalism to go directly from the Langevin Eq. (4.46) to the path integral, without making contact with operators. One follows the same procedure as the example in chapter 2.4.1. I will not show all the steps since terms in the action that come from the deterministic part of Eq. (4.46) are produced in exactly the same way as the example of chapter 2.4.1. I only highlight the difference from that example here, which is the form of the noise term. In our example of 2.4.1, the noise term was the constant, with variance T. Performing the functional integral over the noise field r/r (t) required the Gaussian integral and therefore generated the 0 2 term in the Lagrangian density. For the Stepping-Stone model defined in Eq. (4.46), (4.51) 71TV g Chapter 4. The Stepping-Stone Model Revisited 97 and the integral over r)T becomes ~ 2r [ / P / e x p ( . Vjr ~ ^fTl) dnT cx exp (±%T [fr] ) , which adds the term to the action. The full action is identical to that given in Eq. (4.50). 4.5.4 The neutral regime. For the neutral regime, s = 0. We can take either the Fokker-Planck Equation (4.48), or the Langevin Equation (4.46) and make the change of variables to the difference field, 4>r = 2ft ~ 1. This produces the Fokker-Planck equation, 8W [(f), dt ~ = ^ ^ ^ ^ ( 0 r + e + 0P-e-20 r) + 2 ^ W [ / P ] ' 2nrg ^ dcf>2r For the operator description of the Fokker-Planck equation, we define multiplication by (f)T as the operator 0r, and as the operator ipr. The Liouville operator is £ = TO A ( A r + 2^ r ) + ^ r £ # ( 1 - ^ ) , (4-52) r r and by now we should not be surprised that this operator is completely equivalent to £ECSS in Eq. (4.8), if we identify the scaled creation-annihilation operators c and with (pT and ipr. Chapter 4. The Stepping-Stone Model Revisited 98 4.6 The Broken-Symmetry State: Fate of a single mutant gene. In this section, we will assume that the system has evolved into the broken symmetry state corresponding to fixation of a particular allele, along the lines of our discussion in section 4 of chapter 2. In that section, we found that the non-spatial model, in the thermodynamic limit, was equivalent to a simple branching process. Were we to write the non-spatial model using the creation-annihilation formalism, we would see that the terms of O (-/V-1) in the master equation, Eq. (3.25), that were safely ignored as iV —> oo, correspond to interactions between the otherwise independently branching particles. In the spatial model of this chapter, these terms are 0{pTl), where n is the local density (number per site), and hence cannot be ignored in the thermodynamic limit. We are left with a type of interacting branching process, of which an exact solution appears very difficult. Fortunately, we will be able to find an exact result for the relaxation time r. This quantity is the time scale over which long wavelength perturbations relax and is an important quantity to calculate since it pertains to any avalanche or cascade, induced by a single new mutant, that lasts long enough to spread throughout the population. This is relevant since we are ultimately interested in whether or not this cascade of new alleles significantly changes the composition of the population or indeed approaches fixation. In this section, we will find that r diverges as s —* 0, and the exponent describing this divergence differs from its mean-field counterpart below d = 2. As discussed in the preceding section, we can describe the model using several different formalisms. We will need the path-integral description for this section, however. In the previous section we described how one could accomplish this starting from the Stepping-Stone Langevin equation. We also pointed out the equivalence of this procedure with the same derivation of the path-integral from the occupation number representation. We remind the reader of this latter technique. The details can be found in the example of Chapter 4. The Stepping-Stone Model Revisited 99 section 2.4.2. To get the path integral using the occupation number representation, one starts with the Liouville operator of Eq. (4.8) for the ECSS model. We then use canonical transfor-mations of the kind discussed in section 3.2.2 for Moran's model to eliminate one set of operators. For this section, we will eliminate the operators corresponding to the fixed allele. Assuming a homogeneously distributed population, as was done for the neutral theory, produces identical results for the DCSS model, in that terms that couple the total field dr = ar + br and the ar field we are interested in, do not affect the calculation of the relaxation time below. After Fourier transformation, we are left with the Liouville operator C = - E (mu;k + s) a k a k + sN 1 ^ a l c a k - q a k + q k k,q + E a k + q a k - q a k - (2nN)~1 E 4i a k 2 a kl-qak 2 +q, k,q ki,k2,q where time is in units of the generation time r g , (so r g = 1). For a selectively advan-tageous allele, s > 0, and for a disadvantageous allele, s < 0. The a and af operators here are scaled creation and annihilation operators and therefore represent the allele frequency. To facilitate the calculation, I take the continuum limit of C, and replace the wave-vector sums by integrals. This will in no way change the results. Using the standard procedure in field theory discussed in section 2.4.2, I pass from the Liouville operator, via the coherent-state representation, to the path-integral. To remind the reader, the path-integral is written as Z= j exp (-S[<M), Chapter 4. The Stepping-Stone Model Revisited 100 where the "action" S reads, S [4>, <f>] = j dt j ddr (4>T (t) dt<f> (t)r - L [4>r (t), 4>T (t)]), and where L [0 r (t) , (f)r (£)] obtained from £ [aj, ar] by replacing the a r and a r operators in £ with the fields 0 r (£) and 0 r (£), respectively. The quantity under the action integral, 0 r (t) dt(f) (t)T — L [0 r (t) , 0 r (t)], is called the Lagrangian density. As an alternative, we could stick with our Liouville operator, and use the formalism of "time-ordered" operators [44], but we would obtain the same results. Finally, I perform an additional Fourier transformation from time t to frequency to. This produces the action s [4>> <f>] = J du yA adk (iu + Dk2 + s) 4>K^ (4-53) + S I dude J ddkdd(\ 0k,u;0k-q,u,- £0q, £ + {2ny1 J dude J ddkddq fa,M-w-sKe - (2n) - 1 jdu1du2d£ddkddk2ddc$Kl,ulfa2,u2(t>k1+^^^ Here, and in the rest of this thesis, we'll use d = so ddk — jj^jd- A graphical representation of S [</>, </>] in terms of edges and vertices is given in Figure 4.3. Each vertex represents a term in the action. Lines with outgoing arrows represent 4>T (creation) fields, those with incoming lines represent </>r (annihilation) fields. The single edge (line) without a vertex represents the 0 r 0 r term. Except for this latter term, the label under each vertex gives the coupling constant that multiplies the respective term in the action. The factor {iu + Dk2 + s) 1 under the single edge is the two-point propagator of the "free theory", which comes from the solution of the path-integral for the theory without the interaction terms. In the next section, I give a short description of how these diagrams assist us in the calculations that follow. Chapter 4. The Stepping-Stone Model Revisited 101 Figure 4.3: Vertices representing the action for propagation of new mutant alleles. The propagator of the non-interacting theory represented by the line in the upper left. The vertices representing the interactions are labeled with their respective coupling con-stants. Chapter 4. The Stepping-Stone Model Revisited 102 4.6.1 Diagrammatic Methods in Field Theory. We are unable to solve the full path-integral for the action in Eq. (4.53) or alternatively Figure (4.3). At the end of section 2.4.1 of chapter 2, I mentioned that the procedure for dealing with such "interacting" systems was to break the action up into two parts: a Gaussian part So, that is bilinear in the <j6 and 0 fields, and an interacting part Si containing the higher order terms. A perturbation expansion of a given response function can then be carried out using the propagator from the action So- I will not describe the technical details of how one generates such a perturbation expansion here. The reader is encouraged to consult one of the many texts on quantum field theory or the references cited in this thesis [38] [61] [37] [33] [46] [61]. Here I give an account of the evaluation of such a perturbation expansion using diagrams generated from the vertices from Figure (4.3). A response function is defined by the expectation where for convenience, we have dropped the indices for wavevectors k i , . . . , k m + n and frequencies cui,tom+n (or alternative positions and times). The perturbation expansion for Gnim using the "bare" action So will contain a countable infinite number of terms, each of which can be represented by a Feynman diagram constructed using the vertices and propagator line of Figure (4.3). A term in the perturbation expansion for Gn,m is constructed from these vertices by connecting the outgoing lines of one vertex with the incomimg lines of another vertex in such a way that the total number of incoming lines is equal to m and the total number of outgoing lines is equal to n. The "Feynman rules" are then applied to convert the diagram into an expression. For example, the diagram in Figure (4.4) comes from combining the vertices in the upper right and lower left of Figure (4.3). This diagram has one incoming line and one outgoing line and therefore represents a correction to the two-point propagator. The Feynman rules consist of assigning a Chapter 4. The Stepping-Stone Model Revisited 103 wavevector q and frequency u to each line which then brings a factor of (iu + Dq2 + s) , the propagator of the So action, into the final expression for the perturbation term. Each vertex in the diagram contributes a factor equal to its coupling constant to this expression, i.e. the numbers under the vertices in Figure (4.3). Finally, all internal wave-vectors and frequencies are integrated and wave-vector and frequency are conserved at each vertex. Internal variables are those not associated with the m incoming or n outgoing lines; internal lines have both ends connected to a vertex. Conservation of wave-vector at each vertex means that the sum over the wave-vectors of the incoming lines at that vertex must equal the sum of the wave-vectors of the outgoing lines for that vertex. As an example, let us compute the correction to the two-point propagator represented by the diagram in Figure (4.4). The diagram gives us the following integral, 1 An2 (iu + Dk2 + s) {iuj + Dk2 + s) x li {u - e) + D (k - q) 2 + s i (u + e) + D (k + q) 2 + s The evaluation of this integral gives us the first-order correction to the full propagator, however, we will not use this result. Instead, we will use a similar perturbation expansion for E, the inverse of the relaxation time, in the next section. 4.6.2 Perturbation theory for the relaxation rate E . At the beginning of section 4.2, we pointed out that none of the interaction terms in the Liouville operator produced corrections to the mutation rate fi. Unlike the situation with the mutation rate /x, however, the interaction terms in Eq. (4.53) do produce corrections or modifications to the two-point propagator. The situation is depicted using Feynman diagrams in Figure (4.5). We see that the contribution coming from diagrams with only co-s, k-q Figure 4.4: Feynman diagram representing the 1-loop correction to the two-point prop-agator one "annihilation" vertex (the vertex in the bottom left of Figure (4.3) are O(s), while the contributions from those with m "annihilation" vertices are D(sm). For example, the bottom diagram in Figure (4.5) contributes a term O (s2). Along the same lines, a perturbation expansion for E, the inverse of the relaxation time, is given in Figure (4.6 a). As we are interested in behavior of the relaxation time for the limit of s —»• 0, we might be tempted to use the expansion in a) which uses the bare s in all of propagators of the loop integrals. This expansion leaves out diagrams of O (s2) and higher and therefore is of questionable validity in the desired limit of s —*• 0. We can remedy this however, by using the expansion in Figure (4.6 b), where we have substituted E back into the propagators to get a self-consistent equation, formally called a "Dyson" equation for E. This procedure has the effect of including all the diagrams we have left out. In the limit of s —* 0, the equation for E is exact. If the self-consistent equation can be solved, we're in luck and can extract the asymptotic divergence of r = E _ 1 with s. Chapter 4. The Stepping-Stone Model Revisited 105 Figure 4.5: Perturbation expansion for the two-point propagator. The upper diagrams produce corrections of O (s) to the relaxation rate, while the lower diagrams contributes a correction of O (s2). Chapter 4. The Stepping-Stone Model Revisited 106 Figure 4.6: Perturbation expansion and Dyson equation for E. a) perturbation expansion to O (s) using the bare propagator, b) A self-consistent Dyson equation is derived using the full propagator. Chapter 4. The Stepping-Stone Model Revisited 107 4.6.3 The effective relaxation rate and "screening" of selection in low di-mensions. Using the diagrammatic expansion in Figure (4.6b), the Dyson equation for the inverse relaxation time E = r _ 1 , results in an alternating geometric series E = s ( i - / + / 2 - . . . ) i + r with 1 f00 fA 1 I = — / du I dddq -2n J0 J M w 2 + ( D q + ( q2 + E) Here again, we'll use d = so ddq = -^SL. The integral / is given by (27T)a r-A rfdr, : + E ' i r v i which we have seen before in Eq. (4.29), and evaluates to d=l r l 4 n v ' 2 S £ ) ' I={ b m ( ^ + l ) , d = 2 d = 3 SnnD A 4 7 r 2 n D ' We now proceed to solve the Dyson equation for d = 1, 2 and 3. For d = 1, we have E _ 1 + {Xbi'DE)-" (4'55' for E —> 0 this gives us E = V32n2s2ED Hence the relaxation time diverges as 1 r 32Dn2s2 Chapter 4. The Stepping-Stone Model Revisited 108 This is to be contrasted with the mean-field result in which r ~ s _ 1. For d = 2, the Dyson equation reads E = x ^ 87r£»n m I £ J The solution is n E = SDSTV-'W - i ( - 8 s 7 r ^ e - 8 ^ ) ' where Wfe (x) is the /cth branch of Lambert's W-function, the inverse of /(W) = W e w For small x, W - i (x) ~ In a;, so E = l n /8i>»r + 1 and we find that apart from logarithmic corrections, E is linear in s. Next we consider d — 3. Here we simply have s E = 1 + 47r2nL> Note that the integral / in this case is independent of E, so this result is simply what we would have gotten using the naive perturbation approach of Figure (4.6a). The complete solution of the Dyson equation for d = 1, 2 and 3 is plotted in Figure (4.7). For d = 1, the complete solution, as opposed to the result for small s is 2 - V4 + 512sn2£> h ~ S + 128^73 Note that for larger s, E ~ s, i.e. it approaches the mean-field behaviour. The crossover occurs roughly at s* ~ (64n2D) 1 . One can think of the relaxation time r, or even the selection coefficient as being "renormalized" due to the spatial fluctuations. This is particularly acute for d = 1, as the present results demonstrate. Chapter 4. The Stepping-Stone Model Revisited 109 Divergence of Relaxation Time Selection coefficient s Figure 4.7: Divergence of the relaxation time for d = 1, 2 and 3. For d = 1, the relacxation time diverges as s~2, in constrast to s - 1 in higher dimensions. We also see for d = 2, the slight logarithmic departure from s _ 1. Chapter 4. The Stepping-Stone Model Revisited 110 4.6.4 Survival probabilities: scaling arguments. In the previous section, we found the exponent v that describes the divergence of the characteristic relaxation time as the selection coefficient s —* 0, i.e. r ~ s~v. We would like to know the power law behaviour of other quantities such as the ultimate survival probability of a mutation-induced "avalanche" as a function of s or its time-dependent decay. In fact, a substantial effort was put forth by the author to calculate these quantities directly via several methods: generalization of Kimura's approach [52] to the non-spatial model using the backward equation, closing of the moment equations for n-particle operators and direct summation of Feynman diagrams for small s. Un-fortunately, none of these approaches were successful, which perhaps is a testament to the difficulty of understanding such non-conservative, non-equilibrium models—even the "stepping-stone" discussed in this thesis, which as we have shown, admits some exact results. Fortunately, we can adapt scaling arguments developed for directed percolation (DP) and other growth models to relate the exponents we know from our exact results to the otherwise unknown ones. Here we use scaling arguments put forth by Dickman and Tretyakov [62] for the Domany-Kinzel cellular automaton (DKCA) [63]. The phase diagram of the DKCA possesses a line of critical points that separates an active phase from the absorbing vacuum state. The critical behaviour along this line is that of DP, except at one terminal point where the asymptotic behaviour is know exactly for d = 1 and is identical to that of the Voter model. Spreading of a critical process from a localized source is described by a power laws whose exponents are generally connected by two scaling relations [64]. One of these, known as "hyperscaling', at first sight fails for the DKCA at the terminal critical point corresponding to the Voter model. Dickman and Tretyakov realized that this apparent failure of hyperscaling is due the fact that for the Voter model—and as Chapter 4. The Stepping-Stone Model Revisited 111 we have discussed in the above chapters of this thesis, the "stepping-stone" model—the order parameter is discontinuous at the transition. When this is taken into account, a different hyperscaling relation is obtained. Spreading of a critical process, localized at the origin is described by a set of exponents defined by P(t)cxrs, (4.56) n(t.)ocf, (4.57) and R2{t)ocf, (4.58) where the survival probability P (t) is simply 1 — Po (0> the probability of being in the O-particle state w.r.t. new alleles and the power law behaviour of Eq. (4.56) is valid at long times. In Eq. (4.57), n (t) is the mean population size of new alleles and R 2 (t) in Eq. (4.58), the mean-squared displacement of the new particles. The ultimate survival probability is characterized by the exponent {3, P^ = lim P (t) oc sP. (4.59) t—»oo The scaling hypothesis [64] for spreading from an initial perturbation (in this single new mutant at r = 0, t = 0) postulates the existence of two scaling functions G and (f) defined by p (r, t) ~ t^dzl2G (r2/tz, st1,v) , (4.60) and Ps(t)^rs(i>{t/T) = t-B(f>{tsv). Here p(r,t) is the local order-parameter density, i.e. the density of the new alleles. The form of Eq. (4.60) follows by noting that for s = 0, integrating this equation over all space produces Eq. (4.57). Chapter 4. The Stepping-Stone Model Revisited 112 The existence of the limit in Eq. (4.59) implies (j) (x) ~ x^v = t^^s13 and hence, in order to cancel the -^dependence in Eq. (4.59), B = bv. (4.61) In the case of a selectively advantageous new mutation that reaches fixation, the local density, for fixed r, must approach the stationary density p0 which for our model is simply 1. Since p (r, t) is the average over all trials, for small but finite s we have, lim p (r, t) ~ 5 % . (4.62) t—>oo It follows that G (0, y) ~ y13 for large y. Then in order to cancel the -^dependence in Eq. (4.60), we must have (3/u = y — n, or using the first scaling relation above S + V = f- (4-63) This "hyperscaling" relation depend crucially upon the validity of our assumption that even in the limit that s —+ 0, the order parameter is discontinuous, or in other words, that the symmetry is spontaneously broken. As we have shown, this is true only for d < 2. This situation in not unlike equilibrium critical phenomena, where hyperscaling is not valid above" the upper critical dimension. In the previous section, we found v = 1 for d > 2 and v = 2 for d = 1. For the Stepping-Stone model, this is the only difficult exponent to determine. The other exponents, r\ and z are relatively easy to get. First, we note that when s = 0, the two-point propagator coming from the action of Eq. (4.53), or Figure (4.3), is simply that of a random walk. Therefore, z = 1. The next observation to make, is that, also for s = 0, the expectation of the number of particles is easily shown to be conserved. Chapter 4. The Stepping-Stone Model Revisited 113 (Formally, one can show this by going back to the Liouville operator and evaluating this expectation using the commutation relations). This tells us that rj = 0. We can us the scaling relations above to get the rest. In particular, for d = 1 but P is the same as in higher dimensions P = l. For higher dimensions, 8 takes its "mean-field" value of 1, in the sense that this is what we derived for the Wright-Fisher/Moran model in section 2.4. Figure (4.1) shows simulation results for the time-dependent survival probability in the d = 1 Stepping-Stone model. These results agree with the theoretical calculation of 6 = |. We should remark that the scaling relations above are definitely correct for d = 1, 3, but in d = 2 one expects logarithmic corrections. Chapter 4. The Stepping-Stone Model Revisited 114 10° 10 1 10 2 10 3 10 4 10 5 Time t Figure 4.1: Simulation results for the time-dependent survival probability for d = 1. Shown are simulations for the N = 128 and 256 lattices which clearly show the t'1'2 dependence of the survival propability. Chapter 5 Renormalization Group Analysis 5.1 Introduction. A very significant advance in our understanding of critical phenomena occurred with the advent of Renormalization Group (RG) . The basic ideas of the RG were introduced in a series of papers by Ken Wilson and co-workers in the early 1970's [65]. These ideas have since taken hold in practically every avenue of theoretic physics and have changed the way we understand and look at problems in many-body systems. Prior to the inception of the RG picture, phenomenological assumption about the form of thermodynamic functions near critical points lead to scaling laws, similar to those described at the end of the previous chapter. RG theory not only justifies the scaling laws, it also provides a framework to calculate the exponents in a systematic way, and for difficult problems that defy exact results. RG techniques fall roughly into two classes. The first goes by the name 'fc-space' or 'wave-vector shell' Renormalization Group, and the second 'real-space renormalization'. In this thesis, we will use the wave-vector shell Renormalization Group to analyze the Stepping-Stone model, and recover many of the results presented earlier, but from a very different point of view. The basic idea of the RG is fairly simple. One wishes to know how the model in question behaves at large scales of both space and time. One constructs a series of trans-formations on the model, First, by integrating out the small-scale, fast-time dynamics. This is followed by a re-scaling of space and time. The parameters characterizing the 115 Chapter 5. Renormalization Group Analysis 116 model are thus modified and repeating this procedure generates a flow in the parameter space of the model. The ultimate goal is to find fixed points of this transformation which can be due to 1) the fact that the fluctuations have all averaged out at large scales, or even more interesting, 2) the fluctuations have no characteristic scale and persist from the micro-scale to the macro-scale. The latter possibility is the hallmark of critical be-haviour. Fixed points corresponding to these critical points are unstable with respect to one, or perhaps multiple parameters of the model and it is the linearized eigenvalues of the RG flow that determin the critical exponents. 5.2 A Simple Example of Renormalization: The biased random walk. In this section we will provide a simple illustration of Wilson's wave-vector shell Renor-malization Group. The time-evolution operator for a random walk with diffusion constant D and bias a in the a -^direction is given by U[0,0] = e - s M (5.1) with action S [0,0] = j dt j ddr Udt ~ DV2)4>r - a {^j <j>t. In wave-vector space, one has S [4>, 4>}= j dtddk fa{dt + £>k 2 )0 k + ikia4>k(j)k. (Recall that d = so ddk = Step 1: Separate the k integral in to "fast" and "slow" modes and integrate out modes near the cutoff. This resulting action is Chapter 5. Renormalization Group Analysis 117 S[4>,<f>]=S< [0,0]+S>[^0] with r-A/6 and dt J ddk 0k(<9t + mk2)0k + ihafafa, (5.2) S> [$,</>] = f dt f ddk 4>k(dt + mk2)0k + i^afafa. J JA/b Here b > 1 is our factor describing our wave-vector shell just below the cutoff. Now we integrate out the "fast" modes. This is easy since the action, written in wave-vector space, does not mix "fast" and "slow" modes and is of gaussian form. The integration simply adds a multiplicative constant to the exponential in Eq. (5.1), and the resulting action of the "slow" fields is simply what we had before in Eq. (5.2). Step 2: Re-scale momenta and time. As the objective is to end up with an action that looks the same as before, in partic-ular, we want to restore the cutoff on wave-vector to its original value, so we let ki • b ki, and re-scale time by some as yet undetermined factor bz, i.e. t -> bzt. The resulting action for the "slow" fields becomes S< [4>A] = j' dt j d^b^Mdt + b^Dk^ + ib'-^ak^^. Step 3: Re-scale Fields to find fixed-points of the transformation. Chapter 5. Renormalization Group Analysis 118 To accomplish this we let 0k —> & V k , and now look for values of the exponents z, p and q that produce a fixed point. It is im-portant at this stage to point out the necessity of choosing independent re-scaling factors for the 0k and 0k fields. This is due to the fact that within both the creation/annihilation formalism, and the response function or Fokker-Planck based field theory, most processes, especially those that describe an interaction between diffusing particles (unlike the ones considered in the present section, but see below) contribute two terms to the action that have different numbers of creation operators, or equivalently 0 fields. First we consider the term under the action integral with the time derivative, dtddkdtc])k(j>k - » • bq+p-ddtddkdt4>kcpk. We need to keep this fixed so p + q = d. Next we consider the diffusive and drift terms, which as a result of the above choice give us new constants via the recursion relations . D' = bz~2D (5.3) and a' = ft*"1**. (5.4) From the recursion relations above, we find two different fixed point corresponding to different choices of z. The choice z = 2 produces the 'random walk' fixed point, however, the drift term grows like b at each iteration of the renormalization. In other words, this Chapter 5. Renormalization Group Analysis 119 fixed point is unstable to drift. On the other hand, the choice z = 1 gives us the 'Drift' fixed point, for which the diffusion constant D flows to zero. This is naturally what we expect for drift: at large length and time scales, the effective motion is rectilinear. Using the terminology of the renormalization group, the drift term is relevant, with respect to the random walk fixed point. Albeit simple, this example illustrates some of the basic notions of the renormaliza-tion group and provides a foundation upon which to generalize the procedure to more complicated systems. With further analysis of the above RG flow, one could extract the various exponents that characterize the random walk with and without drift, however in a sense, the results are somewhat trivial given that these are exactly known, i.e the propagators are easily derived due to the non-interacting character of the system. 5.3 Perturbations from the Random Walk Fixed Point. It will be necessary for us to understand more complicated perturbations from the random walk fixed point than considered above. In general, we cannot easily make the separation into "slow" and "fast" fields as above—the resulting action will contain terms that mix "slow" and "fast" modes. It will, however, be very instructive to consider the RG flow for small perturbations from the random walk fixed point. In this manner we can study the stability of the random walk fixed point. In the process, we will discover yet another non-interacting fixed point of central importance in population-like systems. 5.3.1 n-particle to m-particle processes. We start with a random walk action as discussed above, but omitting the drift term. Suppose we add small perturbation corresponding to the stochastic process of going from n-particles to m-particles. This contributes 2 terms to the action Chapter 5. Renormalization Group Analysis 120 //•A / m+n \ dt / c T d k m + n ^ k 1 ^ k 2 - ^ k m 0 k m + 1 0 k m + 2 - 0 k m + „ * I £ k i ) ' and //•A _ / 2n \ dt / (f dki...rf dk 2„ 0 k l 0 k 2 - 0 k „ 0 k n + 1 0 k „ + 2 - 0 k 2 „ < 5 I> • Following the three steps outlined in the example of the preceding section, and requiring that the term with the time-derivative be invariant under re-scaling, we have for the field re-scaling factors p + q = d. The random walk fixed point is accessed by choosing z = 2. For arbitrary p and q, the two terms will scale differently. Accordingly, we will denote the two new coupling constants um]n and uf)„. The resulting recursion relations describing the growth (or decay) of the perturbations are _ u2-(m+n-l)d-\rmp+nq u'm:n ~ u ""m,n ?/(2) _ h2-(2n-l)d+n(p+q) u'm,n u "m,n We can make both coupling constants behave in the same manner if we choose P = d, q = 0. The result being Vm,n = Um]n = umn, and Um,n = 0 ^ ^ um,n (5-5) So the relevance of a perturbation from the random walk fixed point depends only on the number of incoming (f> fields. Chapter 5. Renormalization Group Analysis 121 5.3.2 Relevance of independent decay and branching processes. The simplest process to consider are those which have only one annihilation operator or 4> field, for example random decay (or death), and branching (or birth). Having only one power of <fi in the action means that these processes occur independently of other particles. These processes are always relevant from random walk fixed point, as can be seen by putting n = 1 in Eq. (5.5), which tells us that these terms grow as b2. One could take z = 0 in the above RG procedure and find D' = b'2D and u m , l ~ u m , l , i.e. the diffusion constant D flows to zero. This type of fixed point is non-critical, which is easy to see since it is stable to any perturbation (except, of course other terms with n = 1, which are marginal). 5.3.3 Behavior of coagulation, or annihilation processes. Next we consider processes such as coagulation (A+A —»• A) or annihilation (A+A —> 0). For example, the former would have a term in the Lagrangian u (4>r - 4>l) tf. (5.6) and would produce a RG flow u' = b2~du, so these processes are irrelevant with respect to the random walk fixed point for d > 2. 5.3.4 A paradox. Higher order process (i.e. n > 3) would all be irrelevant at d > 2 and so one might be tempted to conclude that the upper critical dimension dc = 2 for these models. The Chapter 5. Renormalization Group Analysis 122 processes A + A —> A and A + A —> 0 have been quite well studied in the literature where both exact results, RG arguments and simulations clearly establish that dc is in fact equal to 2. Perhaps the most pervasive critical behaviour known in non-equilibrium statistical mechanics is that of directed percolation (DP), also known as Reggeon field theory in particle physics. DP is quite easily formulated using the operator/path-integral tech-niques discussed in this thesis. Its population model analog contains, in addition to births and deaths, density-dependent mortality—or effectively decay, branching and a coagulation term, all of which have been discussed above. So shall we conclude on the basis of the above argument that dc = 2 for DP also? Unfortunately this is quite wrong. DP has been quite well-studied in the literature and is found to have an upper critical dimension dc — 4. This fact is well-established by the breakdown of perturbation theory for d < 4 and non-so-transparent field-theoretic RG results. It seems unfortunate that for DP, dc = 4 does not come out of the simple wave-vector RG argument above. Is there something wrong with the theory or our argument? The resolution of this 'paradox' is that our argument about the stability of the random walk fixed point is in fact correct. What we have neglected, however, is the possibility of another critical fixed point. It is with respect to this fixed point that dc = 4 for directed percolation. 5.4 The Critical Branching Fixed Point. In the preceding section we found that perturbations described by a single annihilation operator or (f> field were classified as 'relevant' with respect to the random walk fixed point. Independent branching is described by such a term and therefore the random walk fixed point is always (in arbitrary dimension) unstable to such a branching process. Chapter 5. Renormalization Group Analysis 123 We should expect from our study of branching processes in appendix C to find another fixed point corresponding to a branching random walk where the branching rate is equal to the decay rate. The qualitative difference between a critical branching random walk and the simple random walk is that the latter conserves the number of particles. In terms of the spreading exponents introduced in the last section of chapter 3, the simple random walk has a survival probability identically equal to 1, so both the exponent (3 in Eq. (4.59) and the exponent 8 in Eq. (4.56) are identically equal to 0. For the critical branching process, however, we know from the exact results in appendix C (see also section 2.4.2) that 8 = 1 and (3 = 1. The mean number of particles in both models is conserved so both have n = 0. Of course, in both cases these exponents are consistent with the scaling laws in Eq. (4.61) and Eq. (4.63). We will now explicitly find the RG fixed point corresponding to the critical branching process and then consider perturbations as was done for the simple random walk above. In addition to the random walk term, the action for the branching random walk has both a decay term and a branching term, i.e. S — SRW + SB + SD with and and the system is critical when 7 = A. It will be convenient to shift the cpr field by l,by defining the field <Pr = 4>r- 1-Chapter 5. Renormalization Group Analysis 124 Grassberger [43] has discussed this shift in the context of the creation-annihilation op-erator formalism, showing that it corresponds to an alternative, but equally valid choice of the scalar product. The resulting action becomes S W = J dt J ddr <f>t{dt - £>V2)0r, JD = j dt J ddr 7 0 2 0 r - (7 - A) 0r0r, SB + SL At the critical point, therefore one finds only the 70 2.0R term. Perform the first two steps of the renormalization procedure, we find the recursion relations D' = b2~zD, and 7 ' = &2-2d+2p+g Once again we must have p + q — d. We can keep both the diffusion constant D and the branch rate 7 fixed by the choices q = 2, P = d-2, z = 2. Thus we have found the fixed point corresponding to the critical branching process. 5.4.1 Perturbations from the Critical Branching fixed point. We now consider perturbations from the Critical Branching fixed point. Following our standard procedure, we look at a perturbation due to a n-particle to m-particle process, and find the coupling constant behaves under our renormalization procedure as Chapter 5. Renormalization Group Analysis 125 u. 'm,n U, (5.7) Let us first consider the coagulation, or density-dependent mortality term, whose coupling constant we denote by U\^- Recall that in the last section, this term was found relevant for d < 2 with respect to the random walk fixed point, and thereby provided us with something of a paradox given the well-known upper critical dimension of DP, dc = 4. Using the above recursion relation, we see that confirming the upper critical dimension for DP. In the next section, when we renormalize the stepping-stone model in the broken symmetry state, the term corresponding to will vanish at the critical point, thus will not play a role in the critical behaviour. It will be necessary, however, to consider the process with m = 2, n = 2. One finds Hence, models that fall into this category will have an upper critical dimension, dc = 2, which agrees with the exact results presented in the last chapter. 5.5 The Gaussian Fixed Point. Before we proceed to renormalize the stepping-stone model, there is one more "easy" fixed point to consider. We shall need it when we consider renormalization of the stepping stone model for the description of the steady-state critical behaviour obtained as the mutation rate goes to zero. It is also important to consider since it represents a convenient way to incorporate dynamics into equilibrium statistical mechanics and has been thoroughly studied in the context of dynamical critical phenomena. [30] [31] [32] Ul,2 = b « 1 , 2 U, '2,2 Chapter 5. Renormalization Group Analysis 126 In equilibrium statistical mechanics one has a free-energy H, which for the present context can be considered a functional of some field 0r. The simplest way to add dynam-ics to the static picture of equilibrium statistical mechanics is to postulate a Langevin equation for (f)T t-'-W + 'M- <5'8> where n (t) is a Gaussian noise term with auto-correlation yl{t)n{t')) = 218{t-tl). Eq. (5.8) is also know as the time-dependent Landau-Ginzburg (TDLG) equation. The action S, or dynamical functional corresponding to Eq. (5.8) is found by looking at the equivalent Fokker-Planck equation, or using the response function method, and is given by 5 = J d t J d ' r ^ - ^ - ^ l (5.9) If we consider H to represent the Gaussian model, we have # = f (V0r)2 + r o ^ , (5.10) and the dynamical functional (e.g. from the response-function formalism) associated with the TDLG equation above is S = Jdtjddr 4>r{dt - DV2)cf>r + r o0 2, - 7 $ . (5.11) Here rnis proportional to the reduced temperature, and hence the critical point occurs at ro = 0. We note that the action in Eq. (5.11) can be viewed as a random walk action with a "source" term. We now have to deal the this "source" term, 4>2. Setting n = 0 in Eq. (5.5) shows that it is relevant with respect to the Random Walk fixed point in all dimensions. Chapter 5. Renormalization Group Analysis 127 Once again, one needs to find the appropriate field re-scaling exponents p and q to produce the RG fixed point. One simply repeat the above arguments for the jcp2. term. One finds the results d-2 P = 2 ' d + 2 q = 5.5.1 Perturbations from the Gaussian Fixed Point. Consider now the effect of our standard perturbation oc 0 ™ 0 r o n the RG flow near the Gaussian fixed point. Using the above values for p and q, we find U/m,n u am,n To make further contact with equilibrium statistical mechanics, we consider the RG flow arising from a </>™+1 perturbation in H. The corresponding term in the dynamical functional, would be oc 0 r 0 n . It is important to note here that within the context of equi-librium statistical mechanics, except for the (p2. term that we have kept invariant under the RG transformation, only terms with one power of 0 r will appear in the dynamical functional. This is evident from Eq. (5.9). The form of the noise term, oc 0 2 , is the result of a constant noise source (i.e. the noise is independent of </>r). The term coming from the free-energy H is given by cpT SHJ^ and only contains the explicit cpr since H [cj)r] does not depend on any <pr. The coupling constant for a 0r</>" term in S we denote by un, flows as u'n = b1+n-^n-VUn from the Gaussian fixed point. In particular, a (p>\ term in H, which corresponds to a 4>r(/>l term in the action, flows as u'3 = bA-du3. Chapter 5. Renormalization Group Analysis 128 from the Gaussian fixed point. In a similar manner, a <fiP term in H, produces the flow u5 = b 2 u5. Both of these results are identical to the standard, non-dynamical RG theory of equilib-rium statistical mechanics, (see, for example [37]). Finally, in the following section, we will need to consider a term in the Lagrangian proportional to (%<^ r- As mentioned above, this kind of term would never occur in the equilibrium theory. Using the above results, this term scales as u2,2 = b2~du2,2, further confirming dc = 2, for the stepping stone model. 5.6 S u m m a r y o f n o n - i n t e r a c t i n g fixed p o i n t s . In the preceding sections of this chapter, we went to great lengths to elaborate the various non-interacting fixed points described by the creation/annihilation formalism or alternatively the response-function/Fokker-Planck formalism. Our reasons for doing this are twofold. First, this distinction has not always been appreciated or even made in the literature. This is perhaps because all of these non-interacting models can be called "gaussian" in a sense—the propagator is of gaussian form. Yet, as we have discussed, each one describes different physics. Although they share some critical ex-ponents (specifically, those describing the two-point correlation function), they differ in other respects, for example, the survival exponents and conservation of particles. More importantly, however, the distinctions we have elucidated above demonstrate both the flexibility and breath of these operator/path-integral formalisms to describe the diverse world of stochastic processes. Chapter 5. Renormalization Group Analysis 129 Let us briefly summarize these fixed points and their stability within the RG picture. First, we have the Random Walk fixed point. A prototypical model would be a random walk (or collection of random walkers) with spontaneous decay. When the decay rate vanishes, the number of particles (in contrast to the mean number or particles) is con-served. The Random Walk fixed point is always (in arbitrary dimension) unstable to a "branching term". When the branching rate equals the decay rate, the relevant fixed point becomes what we have named the Critical Branching Process fixed point. Only the mean number of particles is conserved, and non-zero survival exponents are found. Next we found that the Gaussian fixed point of the TDLG extension of the equilibrium Gaussian model could be viewed as the inclusion of a "source" term to the Random Walk Lagrangian. This term also is relevant in arbitrary dimension, with respect to the Ran-dom Walk fixed point, giving way to the Gaussian fixed point. Indeed, we have uncovered a very interesting equivalence: fluctuations in the TDLG extension of dynamics to equi-librium statistical field theory can be described as a source of spontaneously decaying random walkers. In the Landau-Ginzburg model of equilibrium phase transitions, one has a 0^ term in H and hence the walkers would also be interacting. In general, using the random walker picture, the absolute number of walkers is not conserved. (We should note, however, that it is possible to construct a free-energy in Eq. (5.8)that does conserve the absolute number, by making the source term oc fc2. This model is called "model B" in the literature on dynamical critical phenomena [30][31]. It is quite straight-forward to see that a term like this grows like bd from the Random Walk fixed point.) Finally, we have evaluated the stability of these fixed points to perturbations due to various interactions between particles. Some of these terms, like those found in directed percolation or </>4-theory are relevant for d < 4, however, these term are not present at the critical points of the Stepping-Stone model, nor do they get generated in the renormalization procedure. We have discussed the terms that are present, finding them Chapter 5. Renormalization Group Analysis 130 relevant for d < 2. In the next section, we will go beyond the perturbation analysis in this section and explicitly find the non-trivial fixed points for d < 2 and show how they give rise to "non-classical" critical exponents. 5.7 Renormalization of the Stepping-Stone Model. Using the wave-vector shell RG introduced in the previous sections, we now proceed with full renormalization of the Stepping-Stone Model. First we must distinguish between the two different regimes of the model, already discussed in chapter 3. The first describes the steady-state of the "neutral theory", which occurs when selection is absent and becomes critical as the mutation rate vanishes. Using the results of the last few sections, mere inspection of the Lagrangian shows that the critical behavior for d > 2 is controlled by the Gaussian fixed point. On the other hand, for the description of the broken-symmetry state we have a Lagrangian describing single new alleles, with vanishingly small selective difference, propagating in a "sea" of the fixed allele. From the Langrangian of the broken-symmetry state, we see that the critical behaviour for d > 2 is controlled by the Critical Branching Process fixed point. We also note that d = 2 is marginal for both these situations, therefore one expects logarithmic corrections for d = 2 systems. 5.7.1 The Neutral Regime. Here we follow the Wilson-like procedure [65] [38] introduced in section 4.2. Because of the similar diagrammatic expansions, the results of our renormalization group analysis are somewhat similar to those of the diffusive annihilation (A + A —> 0) and diffusive coagulation (A + A —* A) models.[66] [67] First, we employ the same procedure discussed at the beginning of section 3.5 to get an action 5" from the Liouville operator for the ECSS model given in Eq. (4.8). One again, the procedure is to introduce the fields Chapter 5. Renormalization Group Analysis 131 (ico+Dk2+2|a)"1 Figure 5.2: Vertices representing the Lagrangian of the Neutral regime 0r (t) and (pr (t) and write S as S [4>, 4>]= j dt J ddV 4>r (t) 0t(P (t)r - L [0P (t) , (Pr (t)} , where L [(pT (t) , (pr (£)] obtained from C [ar, ar] by replacing the a r and a r operators of C with the fields (pr (t) and (pr (t), respectively. The resulting action reads, S[$,<t>] = J du J ddk (iu + Dk2 + 2fi) (PK„(PKu] (5.12) +v Jdu J ddk (pk^fp-k-ui -U J du1du2de J rf"dk1d'dk2(idq0kl,Wl0k2,aJ20k1+q>^1+£0fc2-q,a;2-£, which can be represent graphically as the edges and vertices in Figure (5.2). In the previous chapter, we had u = v = (2n) 1 . In the bare model, they are identical, and our choice of time measured in generations made them equal to 1. For the present purposes, we need to distinguish between them to allow for the possibility of distinct RG flows. The separation of fields into "fast" and "slow" modes results in a term that mixes the two and hence the integration over the "fast" modes produces additional contributions to the renormalized coupling constants. We must use our knowledge of field theory to evaluate these contributions—the net Chapter 5. Renormalization Group Analysis 132 < X + + XX + + ccx xxo< + + Figure 5.3: Bubble diagrams contributing to the renormalization of v (above) and u (below). result is precisely what we would arrive at from a standard perturbation expansion, with the important exception that the wave-vector q in the "loop" integrals is restricted to A/6 < |q| < A. The first important observation to make is that applying steps 1 and 2 of our RG procedure to the Lagrangian in Eq. (5.12) or Figure (5.2) does not "renormalize" the mutation rate //. (Obviously, this is also the case for the standard perturbation expan-sion.) So its behavior under the full RG transformation will be identical to its flow from the Gaussian fixed point and grow as b2. Hence, we can forget about it for the rest of the discussion. The second observation is that the RG procedure only can only generate additional terms oc 0m0 n , with n > 2 and m > 3. These terms are easily shown to be irrelevant. This leaves us with only the two original terms and these have similar renor-malizations, which we can perform to all orders by summing the "parquet" diagrams in Figure (5.3), we have already seen in the previous chapter. For the evaluation of the "bubble" integral, we can take the external wave-vector to be Chapter 5. Renormalization Group Analysis 133 - c o , - k A/b<|k|<A Figure 5.4: Diagram representing the 1-loop integral that contributes to the differential RG recursion relation. zero, see Figure (5.4). This is because any fc-dependence generated in the renormalization of u is easily show to be irrelevant. Furthermore, we only need to know the differential recursion relation for u and v to find the RG fixed point and its associated eigenvalue. We will show shortly that in the limit b = el —» 1, diagrams with n loops proportional to ln, therefore only the one-loop diagrams contribute to the differential recursion relation. The net result of our analysis, is the recursion relations u1 = ( l - , where Chapter 5. Renormalization Group Analysis We have 7T A ' -^ln&,d = 2 -^A^,d = 3 ix D b ' We now define I through el = b. In the limit I —> 0, 6 ~ 1 + Z and we have 7T A ' 7T.D ' As mentioned earlier, this tells us the diagrams with n loops are proportional to lr The differential recursion relations for d = 1 read du dl dv dl = u — D-KA' UV DTVA' These equations have the non-trivial fixed point u = v = D T T A , 0. For d = 2, the recursion relations are du u2 dl = ~^D' dv uv dl ~ ~Lhr' These solution to these equations is „ w - , * ( i + £ r ) - . Chapter 5. Renormalization Group Analysis 135 Now the only fixed point is at u* = v* = 0. For d = 3, the differential recursion relations are, once again, qualitatively different, The somewhat surprising result is a line of fixed points with u* = 0 and v dependent on the initial values. The RG flow for all of the above is plotted in Figure (5.5). At first sight the result for d = 3 might seem somewhat strange but a moment's reflection tells us that this is in complete harmony with the exact results from the previous chapter. For d = 2 and 3, we see that the w-term representing the interaction is irrelevant. This term is responsible for the anomalous decay of the two-point correlation function in d = 1. For d = 1 and 2, we see that the f-term, which loosely speaking is the 'source term' for the fluctuations, flows to zero at large length and time scales. This reflects an absence of fluctuations in the steady-state and is precisely what happens in the absorbing state of fixation that occurs in d = 1 and 2, where as we have argued the symmetry is spontaneously broken. For d = 3, v remains finite and the effective description a large scales is that of the Gaussian model. 5.7.2 The Broken-Symmetry State. We now follow a similar procedure and apply the RG to the action that describes prop-agation of single, new mutant alleles. We re-write the action from Eq. (4.53) here, making the parameters of the model explicit by introducing the constant u as was done in the previous sub-section, du H dv uv dl Dir' Chapter 5. Renormalization Group Analysis 136 v Figure 5.5: Renormalization group flow for the neutral theory. Top: The flow for d = 1 shows the non-trivial fixed point at u* = DirA. Top and Centre: v flows to zero for d = 1 and 2 and indicates that the absorbing state is eventually reached. Centre and Bottom: u flows to zero but v reaches a line of fixed points that correspond to the Gaussian model. Chapter 5. Renormalization Group Analysis 137 ?k,uj(Pk,uj S [cj), 0] = j duj j ddk (iu) + Dk2 + su) fa + SU j dude j ddkdd<\ 0k,w0k-q,u)-e0q,£ + 2~nj d u j d e J ddkddq (pk^k-^-e^e J du1duJ2deddkddk2ddq4>kliW14>k2:OJ2(fan + q , u , ! +£0fc2-q,c u 2n Once again, we can graphically represent this action using the vertices and edges of Figure (4.3). Using the same arguments as above, we find the recursion relations s' = b 2 s { l - ^ I d ) , «' = b 2 - d u ( l - ^ I d ) , with Ij, the same as above. For d = 1, we have the differential recursion relations ds su = 2s dl TV DA' du u 2 = u dl TTDA These have the non-trivial fixed point s* = 0, u* = DTTA. For d = 2, the RG flow is ds' us ~dl = S_7r73' du' u2 III = _vr/3' and the only fixed point is at s* = 0, u* — 0. For d = 3, we find ds' 1 suA. ~dl = S~2lrD' Chapter 5. Renormalization Group Analysis 138 dv^ _ _ l u 2 A ~dl ~ ~U~2:KD' and like the previous result, the only fixed point is at s* = 0, u* = 0 . The RG flow for this regime of the Stepping-Stone model, and for d = 1 and 2 is shown in Figure (??). The flow for d = 3 is qualitatively similar to that for d = 2. These results also correspond nicely to what we found in the previous chapter, namely, that the behavior of the model is qualitatively different for d = 1. 5.8 Determining critical exponents from R G eigenvalues. In the previous section, we deduced the RG flow and found non-trivial fixed point for d = 1. Here we show how this result can be used to find v, the relaxation time exponent calculated in the previous chapter. Recall that the relaxation time r is a divergent function of s, T = T(S)^S-U. (5.13) Under a renormalization transformation parametrized by b, time and hence r gets re-scaled by bz. On the other hand, the RG analysis above also tells us how s gets re-scaled. This allows us to write T(s) =bzT(sby*), (5.14) where ys is the eigenvalue coming from the linearized RG flow near the non-trivial fixed point u*. (At the Critical Branching fixed point, ys = 2). We are free to choose b, and making the choice b = s _ 1^ s results in the relation r(s) = b-zlysr{l), ' (5.15) and so we see that v = - . (5.16) Figure 5.6: RG flow for the broken symmetry state, flow around u* determines the critical exponent v = Top: For d = 2. 1, the linearized RG Chapter 5. Renormalization Group Analysis 140 We have related an eigenvalue of the RG to a critical exponent. The RG we performed in the previous section had z = 2. For the non-trivial fixed point, ys is determined from the eigenvalues of the matrix M ds'(u,s) ds'(u,s) ds du du'(u,s) du'(u,s) ds du where s' (u, s) = ds/dl, u' (u, s) = du/dl. One finds (5.17) M = 1 0 0 -1 therefore ys = 1, and hence v = 2, (5.18) (5.19) in complete agreement with the calculation of chapter 4. Chapter 6 Summary of Results and Conclusions 6.1 The Wright-Fisher and Moran Models. In chapter 3 we introduced the Wright-Fisher Model, which for large population size is described by a Fokker-Planck equation with a heterogeneous diffusion term. We also introduced Moran's model, formulated using a master equation, and used a system size expansion to show that for large population size it obeys the same Fokker-Planck equation. For all intensive purposes, therefore, the Moran model is identical to the Wright-Fisher model. By examing the Fokker-Planck equation from the point of view of the effective potential, we showed that both the steady-state and dynamical behaviour were qualitatively different, depending on the value of 2N/i, where N is the population size and [X is the mutation rate. We also showed that for 2Np <C 1, the mean-field picture that results from ignoring the fluctuations becomes irrelevant. This observation is the basis of Kimura's Neutral Theory[8]. When there is no selective advantage of one allele relative to the other, and when the mutation rate vanishes, the eventual fixation of the population to one of the two equivalent absorbing states was shown to be an example of what is termed in statistical physics, a "spontaneously broken symmetry". Broken symmetry in the Wright-Fisher model was also shown to be associated with a kind of critical behaviour, that of the branching process. In contrast to the familiar examples in statistical physics, however, the critical behaviour in the Wright-Fisher model is accompanied by a discontinuous jump of the order parameter. Although we presented 141 Chapter 6. Summary of Results and Conclusions 142 no new results for the Wright-Fisher-Moran model, a review of the main features of this model from perhaps a slightly different point of view than is discussed in the population genetics literature was necessary in preparation for our analysis of the Stepping-Stone model. 6.2 The Stepping-Stone Model. The main results of this thesis were exact calculations for steady-state and dynamical properties of the Stepping-Stone model. We examined two regimes of the model: the neutral regime and the broken symmetry state. 6.2.1 The neutral regime. In the neutral regime there is no selective advantage of one allele over the other, and for non-zero mutatation rate, a steady state is achieved that is characterized by the correla-tion length and the fixation index FST- The fixation index is the two-point correlation function evaluated at vanishing separation of the field that relates the difference between frequencies of the two alleles. In other words, FST is the variance of the frequency differ-ence averaged over the sub-populations. When the mutation rate becomes vanishingly small, the correlation length diverges for all spatial dimensions d, but FST ~*• 1 o n l y for d < 2. This indicates a broken symmetry of the type found in the Wright-Fisher model. The asymptotics of FST, including its time-dependence, were explicitly shown to be related to the return probability of a random walk and are also indicative of critical behaviour. However, for d > 2, the critical behaviour is not associated with symmetry breaking. The above results for Stepping-Stone model are reminiscent of the Voter model, where "consensus" is reached for infinite lattices with d < 2 and is also a consequence of the Chapter 6. Summary of Results and Conclusions 143 eventual return of the random walk[68]. In contrast to the Stepping-Stone model, FST is identically 1 in the Voter model since the site variable consists of a "spin", cr = ±1 . The fixation index FST is actually measured in natural populations (see ch.5 of [69] and references therein). Therefore the Voter model is of limited use for population genetics. Exact results have been given for the dynamics of reactive interfaces in the Voter model that are in rough agreement with our results, in that logarithmic behaviour is observed in d = 2[70][71]. My results are also similar to those found by Sawyer[3] and Nagalaki[4] for the infinite alleles verions of the Stepping-Stone model and were discussed in section 4.4.2. Most measurements of FST f ° r natural populations are interpreted using Wright's formula from his "Island model", a model with infinite range migration[58]. The results I have found for the steady-state value of FST for the Stepping-Stone model suggest that this may not be appropriate for large population size and d < 2. I should point out that natural populations with d > 2 are quite unrealistic. The results of Eq. (4.33) are presented here using the genetics convention for a population size of 2N, i d = 1 FST — \ I l + 8 7 r m n / l n ( ^ ) d = 3 d = 2 (6.1) l+ 4 7 r 2nm' and are to be compared with the result from the Island model, f - = I T ^ ' ( 6 ' 2 ) in Wright's Island model. Only for d = 3 do the results qualitatively agree. In the limit of vanishing mutation rate p, the results for d < 2 are quite different. One should also be careful in applying these results to d = 2 populations for very small p, due to the extremely slow logarithmic relaxation of FST near the critical point. A final Chapter 6. Summary of Results and Conclusions 144 point concerning FST and long-range migration is also in order. Long-range migration, for example Levy flights, will push the critical dimension for symmetry breaking lower than d = 2. (The infinite range migration of the Island model is an extreme limit where symmetry is not broken even for d = 1). Therefore one should expect long-range migration to have a very strong effect on FST for two-dimensional populations. In other words, species with long-range migration should exhibit greater polymorphism than those with finite-range migration. 6.2.2 The broken symmetry state. I have calculated the relaxation time r for long-wavelength perturbations due to rare, new mutant alleles with selective disadvantage s. The relaxation time r diverges with s, r ~ s-". (6.3) The exponent v characterizing the divergence of r with s is identical to the non-spatial result v = 1 for d > 2. For d — 1, however v = 2. For small selective advantage, r in the Stepping-Stone model is higher than that of the non-spatial Wright-Fisher model, the result of a "screening" of selective advantage due to the spatial structure of the population. The shift of r is finite for d = 3, logarithmically dependent on s for d = 2, and very different for d = 1 as indicated by exponent v = 2 in Eq. (6.3). Other critical exponents describing the spread of new alleles at the critical point s = 0 were easily obtained and found to be equal to their random walk counterparts. Scaling arguments were used to find the asymptotic behaviour of the survival probability Ps (t) of the cascade of new neutral alleles. For d = 1 the non-classical value of v leads to the result Chapter 6. Summary of Results and Conclusions 145 while for d > 2, the asymptotic behaviour of P s (t) is that of a critical branching process, Ps(*)~t - 1 . Although the relaxation time and survival probability are characterized by non-classical exponents for d = 1, the ultimate survival probability = lim^oo Ps (£), or fixation probability for alleles with a selective advantage s' grows as Poo ~ (s'f , (6.4) with j3 = 1 for all dimensions. This "classical" value of ft describes the fixation proba-bility in the Wright-Fisher model as well as the ultimate survival probability for super-critical branching processes. One can summarize these results: Selectively disadvanta-geous new alleles tend to last longer and selectively advantageous alleles tend to take longer to reach fixation, and much more for d = 1 populations, but the fixation probabil-ity scales the same in all dimensions as well as in the non-spatial model. Unfortunately, I was unable to calculate the much desired prefactor in Eq. (6.4). To remind the reader, in the Wright-Fisher model P ^ ~ 2s for small s and large population size [54] [52]. This would be a nice factor to know as the fixation probability is a key ingredient to Kimura's neutral theory. 6.3 Renormalization Group For Reaction-Diffusion and Related Stochastic Processes. In preparation for the renormalization group analysis of the Stepping-Stone model, a good deal of effort was spent examining the general structure of the functional integral description of stochastic processes. The complexity of this elegant description, from the point of view of the renormalization group, manisfested itself in the existence of several non-interacting fixed points. This fact has not been appreciated or even noticed in Chapter 6. Summary of Results and Conclusions 146 the literature concerning the application of renormalization group methods to spatially-distributed stochastic models. The need for this distinction between non-interacting fixed points of the RG may be due to the specific wave-vector shell renormalization used by the author and therefore may not be necessary for others type of renormalization group approaches, e.g. real-space RG[72], field-theoretic RG[33], operator-product expansion RG[32] or "dynamic RG"[73]. On the other hand, understanding the structure of the general theory aids the application of RG methods to not only the Stepping-Stone model, but other more complicated population and reaction-diffusion models that defy exact results. 6.3.1 Renormalization to all orders. Here I reiterate my argument of chapter 5 that the renormalization group analysis for the Stepping-Stone model is valid to all orders of perturbation theory and therefore represents a novel, exact application of the renormalization group . This is contrasted with the usual situation in which the RG is an expansion in £ = dc — d. The reason for this is the existence of a small parameter, s that cuts off a large number of diagrams at the critical point s = 0. These diagrams are generated from the s<f>4>2 term represented by the vertex in the bottom left of Figure (4.3). For example, in the model of directed percolation this vertex is non-zero at the critical point and hence the differential RG flow obtained by including only one loop diagrams is only valid for small e = dc — d. 6.3.2 Modern Theory of Stochastic Processes. In addition to finding exact results for the Stepping-Stone model, I have decribed the synthesis of several theoretical formalisms of stochastic processes. These ranged from the classical Fokker-Planck and Langevin equations to the more modern description using operators and path integrals. Using Moran's model and the Stepping-Stone model, I Chapter 6. Summary of Results and Conclusions 147 demonstrated the equivalence of the several of these approaches. Hopefully this exercise will further clarify the extension of these formalisms to more complicated models. Bibliography [1] M. Kimura and G. H. Weiss, Genetics 49, 561 (1964). [2] G. H. Weiss and M. Kimura, J. Appl. Prob. 2, 129 (1965). [3] S. Sawyer, Ann. Prob. 4, 699 (1976). [4] T. Nagylaki, Theoretical Population Biology 10, 70 (1976). [5] M. Kimura, Cold Spring Harbor Symposium 20, 33 (1955). [6] S. Wright, Genetics 16, 97 (1931). [7] R. A. Fisher, Proc. Roy. Soc. Edin. 42, 321 (1922). [8] M. Kimura, The Neutral Theory of Molecular Evolution, Cambridge University Press, 1983. [9] M. Kimura, Scientific American 241, 94 (1979). [10] J. F. Crow, Science 253, 973 (1991). [11] J. A. Coyne, N. H. Barton, and M. Turelli, Evolution 51, 634 (1997). [12] S. Rouhani and N. Barton, Theor. Pop. Biol. 31, 465 (1987). [13] M. Kimura, Ann. Math. Stat. 28, 882 (1957). [14] A. D. Fokker, Ann. Physik 43, 810 (1914). [15] M. Planck, Sitzber. Preufi Acad. Wiss. , 324 (1917). [16] H. Risken, The Fokker-Planck Equation, 2nd ed., Springer, 1989. [17] R. L. Stratonovich, Topics in the Theory of Random Noise, Gordon and Breach, New York, 1963. [18] N. S. Goel and N. Richter-Dyn, Stochastic Models in Biology, Academic, New York, 1974. [19] N. G. van Kampen, Adv. Chem. Phys. 34, 245 (1976). 148 Bibliography 149 [20] N. G. van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, 1981. [21] L. Reichl, A Modern Course in Statistical Physics, Edward Arnold Ltd., 1981. [22] W. Horsthemke and L. Brenig, Zeitschrift fur Physik B 27, 341 (1977). [23] W. Horsthemke, M. Malke-Mansour, and L. Brenig, Zeitschrift fur Physik B 28, 135 (1977). [24] T. G. Kurtz, J. Appl. Prob. 8, 344 (1971). [25] T. G. Kurtz, J. Chem. Phys. 57, 2976 (1972). [26] P. R. Langevin, C. R. Acad. Sci. (Paris) 146, 530 (1908). [27] N. G. van Kampen, J. Stat. Phys. 24, 175 (1981). [28] K. Ito, Proc. Imp. Acad. 20, 519 (1944). [29] R. L. Stratonovich, Conditional Markov Process and Their Application to the Theory of Optimal Control, Elsevier, New York, 1968. [30] S. K. Ma, Modern Theory of Critical Phenomena, Addison-Wesley, 1976. [31] P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49 , 435 (1976). [32] J. Cardy, Scaling and Renormalization in Statistical Physics, Cambridge University Press, 1996. [33] J. Zinn-Justin, Quantum Field Theory and Critical Phenomena,2nd ed., Oxford University Press, 1993. [34] P. C. Martin, E. D. Siggia, and H. H. Rose, Phys. Rev. A 8, 423 (1973). [35] H. K. Janssen, Z. Physik B 23, 377 (1976). [36] R. Bausch, H. K. Janssen, and H. Wagner, Z. Physik B 24, 113 (1976). [37] J. J. Binney, N. J. Dowrick, A. J. Fisher, and M. E. J. Newman, The Theory of Critical Phenomena, Oxford University Press, 1992. [38] N. Goldenfeld, Lectures on Phase Transitions and the Renormalization Group, Addison-Wesley, 1992. [39] R. D. Mattuck, A guide to Feynman diagrams in the many-body problem, 2nd ed., McGraw-Hill, 1976. Bibliography- ISO [40] M. Doi, J. Phys. A 9, 1465 (1976). [41] M. Doi, J. Phys. A 9, 1479 (1976). [42] Y. B. Zel'dovich and A. Ovchinnikov, Sov. Phys. JETP 47, 829 (1978). [43] P. Grassberger and A. de la Torre, Fortschr. Phys. 28, 547 (1980). [44] D. C. Mattis and M. L. Glasser, RMP 70, 979 (1998). [45] A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics, Prentice-Hall Inc., 1963. [46] J. W. Negele and H. Orland, Quantum Many-Particle Systems, Addison-Wesley, 1978. [47] D. L. Hartl and A. G. Clark, Principles of Population Genetics, Sinauer Associates Inc., 1997. [48] R. A. Fisher, The Genetical Theory of Natural Selection, Oxford Claredon Press, 1930. [49] S. Wright, Ecology 26, 415 (1945). [50] P. A. P. Moran, Proc. Camb. Phil. Soc. 54, 60 (1958). [51] M. Kimura, Proc. Natl. Acad. Sci. USA 41, 144 (1955). [52] M. Kimura, Genetics 28, 713 (1962). [53] J. D. Murray, Mathematical Biology, Springer-Verlag, N.Y., 1993. [54] J. B. S. Haldane, Proc. Cambridge Phil. Soc. 23, 838 (1927). [55] G. Malecot, Proc. Fifth Berkeley Symp. Math. Statist. Prob. 4, 317 (1967). [56] S. Itatsu, Nagoy Math. J. 114, 143 (1989). [57] T. Shiga and K. Uchiyama, Probab. Th. Rel. Fields 73, 87 (1986). [58] S. Wright, Genetics 28, 114 (1943). [59] W. Feller, An Introduction to Probability Theory and Its Applications, Wiley, New York, 1966. [60] S. Sawyer, J. Appl. Prob. 16, 482 (1979). [61] M. L. Bellac, Quantum and Statistical Field Theory, Oxford University Press, 1991. Bibliography 151 [62] R. Dickman and A. Y. Tretyakov, Phys. Rev. E 52, 3218 (1995). [63] E. Domany and W. Kinzel, Phys. Rev. Lett. 53, 447 (1984). [64] P. Grassberger and M. Scheunert, Ann. Phys. (NY) 122, 373 (1979). [65] K. G. Wilson and J. Kogut, Phys. Rep. 12C, 75 (1974). [66] L. Peliti, J. Phys. A 19, L365 (1986). [67] M. Droz and L. Sasvari, Phys. Rev. E. 48, 2343 (1993). [68] T. M. Liggett, Interacting Particle Systems, Springer-Verlag, NY, 1985. [69] J. B. Mitton, Selection in Natural Populations, Oxford University Press, 1997. [70] L. Frachebourg and P. L. Krapivsky, Phys. Rev. E. 53, R3009 (1996). [71] E. Ben-Nairn, L. Frachebourg, and P. L. Krapivsky, Phys. Rev. E. 53, 3078 (1996). [72] T. W. Burkhardt and J. M. J. van Leeuwen, Real-Space Renormalization, Springer-Verlag, 1982. [73] A.-L. Barabaasi and H. E. Stanley, Fractal Concepts in Surface Growth, Cambridge University Press, 1995. [74] G. S. Joyce, J. Phys. A. 5, L65 (1972). [75] G. N. Milshtein, Theory Prob. Appl. 23, 396 (1978). [76] J. Bentley, Comm. ACM 28, 245 (1985). Appendix A Random Walks on the Linear, Square and Cubic Lattices. A . l Random Walk Green's Functions. In this appendix, we calculate the Green's function and return probability for continuous time random walks. Most treatments of random walks in the mathematical physics literature deal with discrete time random walks, and even though these have similar asymptotics, there are some subtle differences to continuous time walks. We need the properties of continuous random walks for the analysis of the Stepping-Stone model in chapter 4. Written in terms of creation-annihilation operators, a random walk on a lattice is described by the Liouville operator For hypercubic lattices, r = Yliniei> a n ( ^ the lattice vector e^ is simply I times the unit vector in the ith direction. Here a is the migration rate or jump rate. We identify the a in Eq. (A.l) with D/l2, so that in the continuum limit D becomes the diffusion constant. We will consider a lattice with a total number of sites N = Yii Ni and assume periodic boundary conditions. We will also take Ni = N2 = ... = L. Fourier transforming to wave-vector space using (A.1) r e al = N- —ik-r (A.2) k k in which , k = ^ with l b ; 2TT one finds 152 Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 153 £ = - a ^ w k a k k a k , where ^ = Z (2 - e*ke -e~ e On hypercubic lattices in d dimensions, d k I uk =4^s in 2 -^ - ; 3=1 For small kl, uik ~ k2l2, so in the continuum limit ike ( 2 7 T ) The Heisenberg equations don't couple higher order operators, dak dt dt = [a k ,£ ] = -aukak a k , £ = au)kak So (A-3) (A.4) (A.5) (A.6) (A.7) (A-8) ak{t) = a k(0)e-^ k t, 4(0 = 4(0)e"*'. (A.9) Thus the propagator, in wave-vector space is Ga{k,tf-ti) = ^|a k(t/)4(*i) = (s 0 , ak (0) 4 (0) e -«"x(*/- t 0 - a ^ k ( * / - * i ) /g a k (0) 4 (0) 0 oi, _ g - a u J k ( * / - * i ) (A.10) Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 154 Now we Fourier transform back to get k Written in terms of indices, Ga({n},t) = ^ £ exp (-aukjt) exp (2^) , (A.12) j=l mj=0 ^ d L = N n £ e x p (~ 4 c r f s i n 2 ^jf) e x p ^V^YI) , j=l m,j=0 Y[ £ exp ^ 2atcos ^E^j e x p ^2n—^j . N ]=1 m,=0 The random walk on the infinite lattice corresponds to the limit as L —>• 00. The multiple sum becomes a multiple integral by introduction of the angles dj = 2n™3, with ddj = 2j^. Then I d r-2n Ga ({n} ,t) = 2 TT / ddj exp (-4at sin2 0/2) exp(irij9) d0j, (A.13) (2TT) fJiJo e-2at d r-2-K = [ / dQj exp (2at cos 9) exp (irijO) ddj. (2TT) f^Jo Recalling the integral representation of the modified Bessel function of the first kind 1 r2n In (x) = — / ein9excos9de, (A.14) 2?r JO one has d Ga ({n} ,t) = e~2dat HIn. (2at), (A.15) i=i c2ix 1 fZn In(x) = —\ eineexcos9d6 (A.16) 2TT JO First we note that at any time Ga (r,t) is normalized to 1. This follows from the identities ez = I0 (z) + 2h (z) + 2I2 {z) + (A.17) Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 155 and In(z) = I.n{z). (A. 18) We also see for large r and t, the anticipated result for the diffusion propagator GD (r,t) ~ {4ntD)~d/2 exp (-T2/UD) , (A.19) which follows from the asymptotic properties of the modified Bessel function. A.2 Relation of the first-passage time to the Green's function. Using the Green's function of a continuous-time random walk, Ga (r,i) gives the prob-ability of finding a random walker at the lattice site r at time t, given that the walker started at the origin at t = 0. Since this is not necessarily the first visit to r, we define the first-passage time distribution Fa (r,t) dt, as the probability that the walker reaches r for the first time in the interval [t, t + dt]. The probability of eventually reaching the origin, or the return probability^ obtained by integrating the first-passage time, /oo Fa{0,t)dt. (A.20) _ We need to relate Fa (r,t) to the Green's function Ga (r,t). To do this we note that if the walker is a point r at time t, it must have been there for the first time at some t' < t. Therefore, one has the convolution Ga(r,t)= ! Ga(0,t-t')Fa(r,t')dt', (A.21) Jo which is valid for r ^ 0. For r = 0, we must add a term for the probability P (t) that the walker has never left r = 0, i.e. Ga (r,t) = P(t)6r+ f Ga (0,t - t') Fa {r,t') dt', (A.22) Jo Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 156 Taking the Laplace Transform of Eq. (A.22) for r = 0 results in Fa {0,8) = 1 - (A.23) Ga (0,s) And since Q (0) = Fa (0,s), w-'-im- (A'24) For the Markovian random walks considered in the previous section, the time interval between successive jumps is exponentially distributed with the characteristic time az, where z is the number of nearest neighbors and equal to 2d for hypercubic lattices con-sidered here. The probability of never leaving the origin is therefore P (t) = e~2adt, (A.25) and its Laplace transform The result for Q (0) is then Q ^ 1 2adGa(0,0) ^A'27^ Under a change of time scale, or equivalently, a —> ba, the Laplace transform behaves as Gba(0,s)=b-1Ga(0,s/b), (A.28) Taking b'1 = 2ad, we can write for the return probability Q(0) = l-c^m (A'29) We see that the propagator for a random walk with rate (2d) -1 is quite special, so in what follows we simply drop the subscript and refer to it and its Laplace transform as G and G. The fundamental quantity G (0,0), which is of course independent of the jump rate, gives us the return probability through Eq. (A.29). Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 157 A.3 The Critical Dimension of the Random Walk. The return probability is unity only if the integral of the Green's function G (0,t) diverges. The result (first proved by Polya for a random walk in discreet time) is that for d < 2, Q (0) = 1, but Q (0) = const for d > 2. We will recover this result for our continuous-time random walks by explicitly calculating G (0,0) below. So we see that d = 2 is special for random walks. A.4 Asymptotic Behaviour of the Random Walk Green's Function Ga (0,t). It is not difficult to explicitly calculate Ga (0,s) in d = 1 and d = 2, for the linear and square lattice, respectively. The case of d = 3 is a much harder integral but has been done by Joyce. Nevertheless, we first show how one can extract the asymptotic dependence of Ga (0,s), for s —> 0, in arbitrary dimension, without explicitly calculating the return probability (for d > 2). This is accomplished by looking at the dependence of Ga (0,t) for large t. It is independent of the particular choice of lattice, e.g. triangular, hexagonal, fee, and depends only on dimension d through Ga(0,t)~(at)-d/2. (A.30) We also observes that for small t, Ga (0,t) is bounded, so any divergence from the integral over t, comes from the large t part of the integral. So Ga (0,0) diverges in d = 1 as limr-joo T 1 / 2 , and in d = 2 as l imr^oolnT, but approaches a constant C for d > 2. The asymptotic behaviour of the Laplace transform of Ga (0,t) is obtained from so-called 'Abelian' theorems and only depends on behaviour for large t given above. One Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 158 ln(4a/s) 8o7T c i B r~s «r .<* = 2 • (A.31) finds (4as)"1 / 2,d = 1 Ga (0,s) ~ < , 6^ + HaVlia'd = 3 The constant C in A.31 is of course related to the return probability by c = T^Wy (A'32> A.5 The Laplace Transform of Ga (0,t). Above we defined G (r,t) as the propagator of a random walk with rate (2d)_1and found that its time integral G (0,0) was intimately related to the return probability of the random walk. Here we explicitly calculate the Laplace transform G (0,s) for hypercubic lattices in d = 1,2 and 3. A change of scale will then give us e - ^ = 2 S e ( 0 - 2 s ) ' From above, G (0,t) = e"' [Jo (t/d)}d First one uses the integral representation of the modified Bessel function 1 f2n Io(z) = 7T / d9 e z c o s d , (A.33) 2vr Jo to obtain poo G(0,s) = / e-ste-tIo{t/d)dt (A.34) l /^°° /*27r ( t ^ ^ —^ / dt / dd# exp I - t + - cos 0j - si TT) Jo Jo \ ^ (27T) Integrating over £ results in the multiple integral 1 f 1 G ( 0 ' s ) = 7 T W ^ , i , _ 1 V , r (A.35) (27r) JO S + 1 — d 1 2^ i=i c o s ^ Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 159 A.5 .1 The Linear Lattice. For d = 1, we have 1 f2n 1 G(0,s) = — / -d6 (A.36) V ' ; 2vr./0 s + l-cosfl V ; 1 vM* + 2) The corresponding quantity for the random walk with rate a is Ga {0,a) = 1 (A.37) Vs (« + 4a) We see that for small s, G Q (0,s) - (4asy1/2 . (A.38) For large s, Ga{0,s)~s-\ (A.39) reflecting the fact that Ga (0,s) —» 1 with slope 0 as t —> 0. A.5.2 The Square Lattice. For d = 2, f27T pi-K 1 1 G(0,s) = —J / , , l f , 7 T T ^ 2 (A.40) 4 7 T 2 Jo Jo « + 1 - i (cos 6>i - cos 6>2) 1 /' 2 7 r 1 = — / rfg2 2 7 r ' / o y/(2s + 2 - c o s 0 2 ) 2 - l 1 -K( 1 7T(S + 1) V^ + f, where K (z) is the Complete Ellliptic Integral of the first kind, so Ga (0,s) = . 1 A K (^—) , (A.41) Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 160 For small s G (0,s) ~ In (2V2) /TT + In (s"1) /2vr (A.42) and hence Ga (0,3) ~ (4«TT)- 1 ^ln (2^) + I In (4a/s)^ (A.43) A.5.3 The Cubic Lattice. For d = 3, we have one of the three so-called 'Watson integrals' 1 r-1-K p2n p2-K 1 G ^ = ^ / / / • 1 w A a r:d61dB3de3 (AAA) ovrd Jo Jo Jo s + 1 _ 3 ( c o s Vi - c o s #2 - cos 03) The other two 'Watson integrals' give the corresponding random walk propagator for the body centered and face centered cubic lattices. These integrals originally were derived to analyze the discreet time random walk, but as we have found, also pertain to the continuous time random walk discussed in this thesis. A closed form expressions for this integral has been given by Joyce [74], with and 4 1 V 1 ~ I2* G ( ° ' s ) = 1 K ) *(*-) , (A.45) kl = \ ± ^ 2 V 4 - x2 - - (2 - x2) VT^x~2 , (A.46) x (s + l)2 + l-Ws(3s + 2)(3s + A)(s + 2) 1 2(s + l) 2 ' [ ' ' Xi X2 Xi~ 1 For small s, G(0,s) ~C + By/s, Appendix A. Random Walks on the Linear, Square and Cubic Lattices. 161 where C = *gK ^ 1 -l-V6-V3) K (jmV6-Vz) (A.48) = 1.516386059 and B = |c = 2.2745790885 (A.49) Then for d = 3 and small s , we have Appendix B Computer Simulations This appendix provides additional details on the computer simulations presented in chap-ter 3 and 4. These simulations were written using C++ and compiled for the Windows, Unix and Linux operating systems. B . l Stochastic Differential Equations. Stochastic differential equations (s.d.e.'s) are easily simulated on the computer using techniques similar to those developed for integration of ordinary differential equations (o.d.e.'s). The only complications involve questions of accuracy and efficiency. The simulations in Figure (3.5) were produced using a method of second-order accuracy de-scribed by Milshtein.[75]. B.2 Continuous-time Reaction Diffusion Processes. For simulations of the Stepping-Stone model, a C++ library for the simulation of continu-ous time reaction-diffusion processes was developed by the author. This library consists of classes representing various fundamental reaction processes (for example A + B —»• A + A, A + B —> B + B, A —> A + A and B 0), in addition to random walk processes (Ar —> AT+e), and exchange processes (Ar + Br+e —» Ar+e + Br). Not all of these pro-cesses were necessary for the Stepping-Stone model, however, the library was designed in such a way that programs for the simulation of other lattice reaction-diffusion processes 162 Appendix B. Computer Simulations 163 could easily be implemented without rewriting the core code. We will not present all of the details of the class library in this appendix. Instead, we will point out the important differences of our simulations for continuous time reaction-diffusion models with those commonly used in statistical physics. The main differences are two-fold. First of all, most simulations in the literature use a spin-like variable at each site which has at most several possible values. The theoretical formulation of birth-death type processes we use is the occupation number formalism, which describes classical "Bosons". By this we mean that there can be an arbitrary number at each site. Therefore we need simulations where the site variable is an integer, and more importantly, not fixed to a particular value. Two-species models have two integers for site variables. Furthermore, the fundamental processes divide naturally into two classes: those that are purely on-site—the reaction processes, and those that couple sites—diffusion and exchange. The second major difference involves time. In the standard Monte-Carlo method, time is discrete and measured in steps. At each step, one tries to change the configu-ration according to some allowed "rules". The new configuration is either accepted or discarded according to probabilities that depend on the difference of the new configu-ration with the old configuration. For example, when simulating equilibrium systems like the Ising model, at each step a site is selected and the spin is flipped. If the new configuration has a lower energy (AE < 0), it is accepted. If it has a higher energy (AE < 0), it is accepted with a probability e ~ A E l k b T . Once again, we strive to make the simulation fit the theoretical description, which is for continuous time Markov processes where the "event" times are exponentially distributed with a characteristic time for each event that depends on the local configuration. Using the Monte-Carlo method, events cannot happen simultaneously. Furthermore, the Monte-Carlo method is inefficient for simulations near a critical point where 'critical slowing down' makes the rejection rate extremely high. In reaction-diffusion models, the critical point is often associated with Appendix B. Computer Simulations 164 an absorbing state, so using the Monte-Carlo method is wasteful because "dead" sites are still picked at random and rejected. For these reasons, a technique was developed for the simulation of continuous time processes. The main idea for the simulation engine, was to keep track of the next event for each site. The simulation is first initialized by iterating through the lattice, site by site, and picking an allowed event from set of the fundamental processes that specify the model. The occurence time of the event is also recorded and the times of the events from all the sites are sorted. After this is completed, the simulation begins by executing the earliest timed event. This event changes the configuration of its associated site, in addition to connected sites, depending on the nature of the event. A new event is generated for this site as well as those coupled to it. The times are again sorted and the procedure is iterated. The reordering of times, especially for the whole lattice seems somewhat inefficient: the time cost in sorting an unsorted collection of N items grows as O (N2) . The solution is to use a "heap", or "priority queue" to maintain the items such that the earliest time is on top. A heap is a special type of binary tree that maintains a set of elements under the operations of inserting new elements and extracting the smallest element [76]. The worst case initial sort using a heap is O (N log N) , but more importantly, the operation of removing the smallest, inserting a new element and reordering the tree requires a time cost of O (log N). This provides an efficient way to implement the continuous time simulations. For this purpose, a class known as a "dynamic priority que" was developed as the core component of the simulation engine. Appendix C The Branching Process The branching process is perhaps the simplest stochastic birth-death process, however, it is necessary to understand, its properties, especially its critical behaviour. This reason is that many critical points in stochastic population models behave as a critical branching processes in the thermodynamic or large system size limit. A even more compelling motivation is that the critical points of spatial versions of these models above their upper critical dimension exhibit the same critical exponents as the simple branching process. Thus, the RG fixed point corresponding to the critical branching process plays the same role in these transitions as does the Gaussian fixed point in equilibrium statistical mechanics. C l The model. Consider a system of 'A'-particles or individuals subject to the Poisson-distributed ran-dom events Death : A 0 rate A Birth : A -> 2A rate 7 The state of zero individuals is an absorbing state, so the system hits this state it stays there. 165 Appendix C. The Branching Process 166 The corresponding master equation is dPn (t) dt = 7 (ra - 1) P n _ i (t) + A (ra + 1) P n + i (£) - (A + 7) nPn (t) (C.l) Suppose we start with one individual at t = 0. We'd like to know something about -the sizes and durations of the resulting cascades. The interesting quantities therefore are 1. The survival probability. Ps (t) = l- Po (t) , or the number of cascades (or "avalanches") lasting longer than time t, and in particular, its asymptotic value Poo = lim P s (t) t—>oo 2. Cumulative size distribution n(s), the number of cascades larger than size s, or the probability that a cascade involves at least s individuals. One finds critical behaviour when A = 7 — A = 0, i.e. all of these quantities behave as power laws Ps{t) - t~s Poo - A* n (s) ~ s~a C.2 Survival probability. The survival probability can be found by using the generating function 00 F{z)t) = YJPn{t)zn (C.2) n=0 Appendix C. The Branching Process 167 From the generating function, we can find the moments by taking derivative w.r.t z and then setting z = 1 . For example dF{z) (n) = dz z = l and often the moments < „ ) - < B ) _ _ A 2 z = l = (n(n - l)(n - 2) • • • (ra - m + 1)) = dmF(z) dzT' z = l are referred to as 'factorial' moments. From the factorial moments, one can always get back the probabilities m=n (C.3) (C.4) We can use Eq. (Cl) to obtain an equation for the generating function for this birth-death process. By multiplying by zn and summing over ra, one finds This can be solved using the method of characteristics. (C.5) C.2.1 The method of characteristics. Let [z (r), t (r)] be some curve in the (z, t) plane. The change in F along this curve is (C.6) dF__df_dt_ dF dz dr dt dr dz dr We want to choose a curve such that S ) / ( £ ) - - e - ^ - * > Appendix C. The Branching Process 168 Eq. (C.5) implies F (z,t) is constant along this characteristic. Integrating dz - = - ( ; - l ) ( 7 2 - A ) we find C = z - 1 3 ( 7 - A ) t 7 Z — A where C is some constant of the curve z (t). Then F (z,t) must be some function of C.If we consider initial conditions such that at t = 0, the number of individuals was exactly m, or Pn (0) = 5n>m, we have F(z,0) = zm One finds the solution - e ( A - # t ( 7 2 _ A) + A ( l - z ) ' F(z,t) = (C.7) This complicated looking expression tells us everything we need to know about this birth death process. Expanding in a power series, we can find Pn (t). The survival probability is (t) = A — 7 e(A-7)*A 7 (C.8) For 7 < A, Poo =0 the population becomes extinct. For A < 7 ^ = 1-A 7 For the critical case, A = 7, one finds A 7 (C.9) (C.10) (t) = 1 + 7* ( C H ) Appendix C. The Branching Process 169 The result for the critical point can also be found by solving the (factorial) moment equations. Using Eq. (C.3), one can show that d(j>m (t) _ *yk\ dt ~(k-2)\nm-1 and that for 1-initial particle boundary conditions , these equations can be recursively integrated to find 4>m (t) = ro! (/ft)""1 , which is then substituted into Eq. (C.4) to obtain Eq. (C.ll). C.2.2 The backward equat ion. There's actually an easier way to get the duration results. For this method, we need conditional probabilities. Pn/m (^2, ti) = prob. of finding n at time t<i given ro at time t\. The forward equation reads, — = > Wnyrni/m \ t2,ti) — Wnitnrn/m [t2,ti) , 0^2 , TV or with t = t2 — ti, dPn/m {t) \ ^ r . , s , > ^ = 2-^Wn<n n'/m ^ ' ' Wn'^nlm \t) • ri The backward equation is, dPn/m{t2,ti) D /, , >. . D /. , N — = — / Wm1 ,m"n/m' \^2,tl) + Wmi^m^n/m {l2,^l) , Appendix C. The Branching Process 170 or with t = t2 — t\, dPn/m jt) _ V - ^ p / ,x _ p / , v ^ — / J l^m',m* n/m1 \L) wm',m-*n/m \<') • m' Note that we can obtain the master equation by multiplying the both sides of the forward equation by Pm and summing over TO, so that one only need worry about state probabil-ities at time t2 and not conditional probabilities. For this simple branching process, the backward equation for 1-initial individual boundary conditions is dPn,i (t) dt Using the generating functions = - (A + 7 ) Pn/1 + 7 P n / 2 + XPn/0 n=0 00 F2(z) = yPn/2(t)zn n=0 the backward equation becomes dFi {z,t) = - ( A + 7 ) F ! + 7 F 2 (z,t) + A, dt but since Fm {z,t) = [F1(z,t)}7 d F l { Z ' t } = - ( A + 7) ^ + 7 ^ ( 3 , 0 +A, at subject to the initial condition Fi{z,0) = z hence e ( A - 7 ) t ( 7 Z _ A) + A ( l - z ) e(A-7)t (72; — A) + 7 (1 — z) Appendix C. The Branching Process Figure C.l: A typical tree from a branching process C.3 Size distribution. In order to find the size distribution, we look at a generic tree The probability of a tree with n leaves is # of trees Qr with n leaves X [Prob(death)]# l e a v e s x [Prob(branch)]# = (Cayley's formula) x 2n- 1 x A + 7 y v A + 7 7 n - l 2n - 1 n A7 (A + 7 ) ' A + 7 7 Using Stirling's formula Qr - 3 / 2 —=n n 2v^ L(A + 7) . 4A7 A + 7 2 7 or Appendix C. The Branching Process 172 with the characteristic size C = r l - A (C.13) So the cumulative size distribution ra(s) ~ s " 1 / 2 We briefly mention another way to calculate the size distribution by using a slightly modified model. Suppose we write a master equation for the process corresponding to Death : A—> R rate A , Birth : A -»• 2A rate 7 , where R is a 'removed' individual. Using a two species generating function F(y,z) = y/Pn,lylzn, n,l one can use either the backward or forward equation and the relevant method above to solve for the t —* 00, R probabilities, defined by Qt = J2n Pn,i-
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Broken symmetry and critical phenomena in population...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Broken symmetry and critical phenomena in population genetics : the stepping-stone model Duty, Timothy Lee 2000
pdf
Page Metadata
Item Metadata
Title | Broken symmetry and critical phenomena in population genetics : the stepping-stone model |
Creator |
Duty, Timothy Lee |
Date Issued | 2000 |
Description | In this thesis, I study the behaviour of the Stepping-Stone model: a stochastic model from theoretical population genetics. It was introduced by Kimura and Weiss [1][2] as a simple model to investigate the interplay of the evolutionary processes of random genetic drift, mutation, migration and selection. In particular, they were interested in the behaviour of spatially-structured populations when these processes were mediated by local interactions. From the point of view of statistical physics, the Stepping-Stone model can be viewed as a 2-species, non-equilibrium reaction-diffusion model with spatial degrees of freedom and unique fluctuations that arise from irreversible processes at the microscale. My thesis begins with a brief overview of the theory of stochastic processes that includes both the classical treatment using master equations, the Fokker-Planck equation and the Langevin equation, and modern formalisms that map these equations to operators and functional integrals. This is followed by a discussion of the steady-state and critical phenomena in the Wright-Fisher and Moran models, the non-spatial predecessors of the Stepping-Stone model. Critical phenomena in these models is shown to be associated with the breaking of a discrete symmetry that results in a discontinuous order parameter. Next, the Stepping-Stone model is introduced and reformulated using operators and pathintegrals. Kimura and Weiss were able to solve for the steady-state of the Stepping-Stone model, under certain conditions, but the dynamics of the model remained elusive. This model is related however, to the "Voter" model, which has been well-studied and is known to have a critical spatial dimension of 2 in that only for d ≤ 2 does the system asymptotically reach one of the 2 degenerate absorbing states. I extend the results of Kimura and Weiss and obtain exact results for the dynamics of the Stepping-Stone model. Two regimes of the model are analyzed: 1) the neutral regime, where selection has vanished and mutation is viewed as the control parameter. Here a steady state exists, characterized by a correlation length that diverges as the mutation rate, μ, becomes zero; and 2) the broken symmetry or "fixed" state, where selection is small but finite, and mutations so rare that the relevant description concerns the dynamics of "avalanches" or "cascades" of new alleles induced by the initial mutation and perturbing the system from the absorbing state. A unique kind of dynamical critical phenomenon occurs when selective advantage and mutation become negligible. Like the Voter model, it is qualitatively different both above and below a critical dimension d[sub c] = 2. For spatial dimension d ≤ 2 , the critical behaviour is associated with the breaking of a discrete symmetry and corresponds to fixation of one of the two alleles (genotypes). Symmetry breaking in the Stepping-Stone model is shown to be a consequence of the asymptotic return probability of a random walk—identically one for d ≤ 2 and strictly less than one for d > 2. In addition to the correlation length, the steady-state of the neutral regime is further characterized by a measure of variance defined as the amplitude of the two-point correlation function evaluated at vanishing separation. In genetics this measure of variance is known as F[sub ST], the fixation index. Exact results are derived for both the steady-state value of FST, and its asymptotic time-dependence at the critical point, which approaches 1 in both d = 1 and d = 2. At the critical point, 1 — F[sub ST]~ t[sub -1/2] for d = 1. For d = 2, a much slower decay is found, 1 — F[sub ST] ~ 1/ln(t). The d = 3 critical behaviour of F[sub ST] is that it approaches a constant C < 1 as (C — F[sub ST]) ~ t[sub -1/2]. The constant C is non-universal and related to the return probability of a random walk. It approaches 1 for very large values of the dimensionless constant K = [Λ/(2 π² Dn[sub Tg)] where Λ⁻¹ is the spatial scale of the interaction, rg is the generation time, D is the diffusion constant and n is the population density. A related infinite-alleles model was studied by both Sawyer [3] and Nagalaki[4]. These authors found similar asymptotics for the probability that two randomly choses individuals are genetically identical under certain assumptions. Sawyer[3] also proved that fixation in the infinite alleles model occurs iff the random walk followed by the individuals is recurrent. The broken symmetry or "fixed" regime of the Stepping-Stone model has been explored from the point of view of survival of rare mutant alleles, here parameterized by a coefficient of selection s. The exponent u, governing the divergence of the relaxation time as s —> 0 is calculated and found to be v = 2 for d — 1, while for d ≥ 2 it is given by the mean field value of v = 1. Two other exponents for critical spreading processes are determined and scaling arguments are presented and used to find the decay exponents that characterize the time-dependent survival probability and its asymptotic value for very long times. These results also establish the upper critical dimension of the model, d[sub c] = 2. Finally, these results are rederived and supplemented by a dynamical renormalization group (RG) analysis. The critical behaviour in d = 1 is found in both regimes to be controlled by non-trivial fixed points. The critical behaviour of the broken symmetry regime for d ≥ 2 and its RG fixed point is that of a critical branching process. For the neutral regime, however, the renormalization group flow is qualitatively different in rf = 2 and 3, reflecting the existence of broken symmetry for rf = 2. The RG flow for the rf = 3 neutral regime contains a line of fixed points with the effective description at large scales given by a Gaussian version of the time-dependent Landau-Ginsburg model. Finally, the renormalization results are found to be valid to all orders of perturbation theory. |
Extent | 7917538 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-27 |
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. |
IsShownAt | 10.14288/1.0085696 |
URI | http://hdl.handle.net/2429/11348 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2000-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2000-56536X.pdf [ 7.55MB ]
- Metadata
- JSON: 831-1.0085696.json
- JSON-LD: 831-1.0085696-ld.json
- RDF/XML (Pretty): 831-1.0085696-rdf.xml
- RDF/JSON: 831-1.0085696-rdf.json
- Turtle: 831-1.0085696-turtle.txt
- N-Triples: 831-1.0085696-rdf-ntriples.txt
- Original Record: 831-1.0085696-source.json
- Full Text
- 831-1.0085696-fulltext.txt
- Citation
- 831-1.0085696.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-0085696/manifest