Models of genetic drift as limiting forms of the Lotka-Volterra competition model

The relationship between the Moran model and stochastic Lotka-Volterra competition (SLVC) model is explored via timescale separation arguments. For neutral systems the two are found to be equivalent at long times. For systems with selective pressure, their behavior differs. It is argued that the SLVC is preferable to the Moran model since in the SLVC population size is regulated by competition, rather than arbitrarily fixed as in the Moran model. As a consequence, ambiguities found in the Moran model associated with the introduction of more complex processes, such as selection, are avoided.

The modeling of genetic drift -the mechanism by which the genetic makeup of a population can change due to random fluctuations -is frequently viewed in a different way to the other key genetic processes, such as mutation, migration and selection. It requires an inherently stochastic approach, and many texts explicitly or implicitly take the view that it is best explained in a way that illustrates the key idea, avoiding too much realism that might result in complex stochastic models. The standard way to do this is through the use of the Wright-Fisher model in which a population of N individuals carrying one of two alleles is sampled at a given time t in order to create a new generation, also with N individuals, at time t + 1 [1]. It is recognised that this does little more than illustrate how genetic drift can arise in populations when N is finite, and some authors do preface their presentation of the model with a discussion of its possible relation to more realistic situations [2]. Of course, the basic model may be extended in many ways: by using the diffusion approximation [3]; by introducing an effective population size to replace variations in population size, for instance, by an average [4,5]; by replacing non-overlapping generations by overlapping generations (the Moran model [6,7]), amongst others. Yet in all these examples, the effect is to retain the fixed size of the population.
In the case of the Moran model, which has been extensively used in many areas, especially in the last decade or so [8][9][10], the assumption that N is fixed is hardly ever questioned. However this results in a single rate parameter encompassing birth and death, making it impossible to tease apart effects resulting from the various processes. The ambiguities inherent in this approach are particularly apparent when trying to include the effects of selection in the model.
In this Letter we will advocate a different starting point which allows us to address some of these questions. We adopt a more ecologically-orientated approach and begin from a population of n 1 haploid individuals which carry allele A 1 and n 2 haploid individuals which carry allele A 2 . They will reproduce at rates b 1 and b 2 respectively and die at rates d 1 and d 2 . We will also allow for competition between individuals of type A i and A j , at a rate c ij . This will tend to regulate the population size, without imposing the condition n 1 + n 2 = N . The model will be formulated as an individual based model (IBM), since the stochastic aspects are central to the discussion. Beginning from an IBM used to be less common in ecology than in population genetics, to some extent because the effects of stochasticity were felt to be less important, however it is now becoming more widespread [11,12].
Although it is quite easy to sketch out this idea, showing how precisely this model relates to conventional models in population genetics, such as the Moran model, is not so straightforward, and we are aware of only a few studies which touch on this issue. In [13], an exact mapping between the two models was sought via the introduction of unconventional fitness weightings, while [14][15][16] were concerned with significant deviations from Moran phenomenology. To provide a systematic understanding of the relationship between the two approaches, we will apply an approximation procedure which we recently developed based on the elimination of fast modes [17,18].
The system we will investigate will be well-mixed, and so its state will be completely specified by n = (n 1 , n 2 ). This state will be able to change because of births, deaths or competition between individuals, defined by the rules given in Fig. 1. To define a dynamics we need to specify the rates at which the allowed changes in Fig. 1 take place. We assume mass action for the competitive interactions leading to transition rates given by The parameter V is not the total number of individuals in the system, which is free to vary. Rather it is a measure of the size of the system. For instance, in a terrestrial ecology it would be the area of land which the individuals inhabit. So typically it would be an area or a volume, but its precise value or even its dimensions can be left unspecified, as they can be absorbed into the rates b i , d i and c ij . The situation is analogous to classical thermodynamics, where if boundary effects are neglected, extensive quantities such as n i scale with the volume of the system, so that if the system doubles in size, the expected values of these quantities would also double [19]. The probability of finding the system in state n at time t, P n (t), may be found from the master equation [20] dP where ν µ describes how many individuals of one type are transformed during the reaction µ = 1, . . . , 4. So, Equations (1) and (2), together with an initial condition for P n , allow us in principle to find P n (t) for all t.
In practice, the master equation is intractable. To make progress the diffusion approximation is made, that is V is assumed sufficiently large that x i ≡ n i /V is approximately continuous [21]. The variables x i then specify the state of the system. We can rewrite the master equation (2) in terms of the variables x i as (3) Here P (x, t) is P n with n replaced by V x and f µ (x) is similarly equal to T µ (V x + ν µ |V x). For large V we may expand the functions P and f as a Taylor series about x to obtain the Fokker-Planck equation (FPE) [22] ∂P where τ = t/V is a rescaled time and where we have neglected higher order terms in V −1 . The functions A i and B ij can be expressed in terms of the ν i,µ and the f µ as [23] A where i, j = 1, 2. The diffusion approximation, made popular by Kimura and others in the context of population genetics [21], is usually expressed in the form of an FPE such as Eq. (4), however for our purposes it is preferable to work with the entirely equivalent Itō stochastic differential equations (SDEs) [22] dx where η i (τ ) is a Gaussian noise with The precise form of the functions A i (x) and B ij (x) for the system of interest to us can be read off from Eq. (5) using Eq. (1) and the ν i,µ given earlier. One finds where i = 1, 2 and j = i. The off-diagonal entries of the matrix B are equal to zero. In the limit V → ∞, Eq. (6) reduces to the two deterministic differential equations dx i /dτ = A i (x), with A i (x) given by Eq. (8), which are the familiar Lotka-Volterra equations for two competing species [24,25]. We begin the analysis by assuming that individuals of type A 1 and A 2 have equal fitness. Thus the theory is neutral, and A 1 and A 2 have equal birth, death and competition rates: Simulations of the original IBM defined by Fig. 1 and Eq.(1) are shown in Fig. 2, where it is seen that the trajectories quickly collapse onto a line in the x 1 -x 2 plane. This can be understood by first considering the deterministic trajectories (shown in gray in Fig. 2). We begin by looking for fixed points of the dynamics by setting A i (x) = 0, i = 1, 2. Taking the combinations A 1 ± A 2 we find that the fixed points are solutions of the two equations We see that, apart from the trivial fixed point x 1 = x 2 = 0, there is a line of fixed points given by This is the equation of the blue line shown in Fig. 2. Further insight can be gained by calculating the Jacobian at points on Eq. (10). One finds that it has eigenvalues λ (1) = 0 and λ (2) = −(b 0 − d 0 ), with corresponding eigenvectors and where u (i) and v (i) are respectively the left-and righteigenvectors corresponding to the eigenvalue λ (i) (normalised such that , and x 2 is given by Eq. (10). This shows that Eq. (10) defines a center manifold to which the deterministic system quickly collapses; there is then no further motion along this line. The timescale for the collapse of this fast mode is given by |λ (2) The stochastic dynamics (shown in red in Fig. 2) are dominated by the deterministic dynamics far from the center manifold, and there is a rapid collapse to its vicinity. Fluctuations taking the system too far away from the center manifold are similarly countered by the deterministic dynamics dragging the system back. The net result is a drift along the center manifold until either of the axes are reached and fixation of one of the types is achieved. This effect has also been noted and exploited in [26,27] under the investigation of the evolution of dispersion. To encapsulate this behavior in a mathematical form we use a methodology which we have recently used to reduce stochastic metapopulation models to effective models involving only one island [17,18]. The essential idea is to restrict the system to the center manifold, and to obtain the effective stochastic dynamics along the manifold by applying the projection operator to the SDEs (6) to eliminate the fast mode, keeping the slow mode intact. To carry out this program, we first rescale the x i and time in order to eliminate various constants from the calculation. Writing the equation of the center manifold becomes y 2 = 1 − y 1 and the non-zero eigenvalue of the Jacobian is now equal to −1. Applying the projection operator to the first term on the right-hand side of Eq. (6) gives P A = 0, confirming that there is no deterministic dynamics along the center manifold. We denote the coordinate along the center manifold as z, and choose this to be equal to y 1 , although of course many other choices are possible. Sinceẏ 2 = −ẏ 1 on y 2 = 1 − y 1 , application of the projection operator to the left-hand side of Eq. (6), and to the noise term on the right-hand side givesż = ζ(τ )/ √ V , where ζ(τ ) = P 11 η 1 (τ ) + P 12 η 2 (τ ). It should be noted that since the projection operator depends on y 1 , the direction of the dominant noise component changes, as can be seen from Fig. 2. From the properties of η i , we see that the effective noise ζ is Gaussian with zero mean and with correlator with the B ij being evaluated on y 2 = 1−y 1 . A calculation of the term in square brackets in Eq. (14), allows us to arrive at the following form for the SDE describing the dynamics after the fast-mode elimination: zero mean and correlator (16) If we define N = (b 0 − d 0 )V /c 0 (the size of the population on the center manifold, Eq. (10)) then Eqs. (15) and (16) are exactly the Moran model in rescaled timē where z is the fraction of type A 1 alleles and N the total population size. We therefore conclude that the neutral form of the SLVC model reduces to precisely the Moran model, under the fast-mode elimination procedure described in [17,18].
It is now natural to ask what model is obtained by the elimination of the fast modes of the non-neutral SLVC model. As usual in population genetics, we work to linear order in the selection strength, s, and so begin by writing where is a small parameter which will later be related to s. The constants β i , δ i and γ ij are assumed to be of order one. Although for = 0, there will not be a center manifold, we still expect there to be separation of timescales which will allow us to identify fast and slow variables. We pick out the slow subspace, and so eliminate the fast deterministic dynamics, by setting the product u (2) · A(x) equal to zero [17]. This leads to an equation of the form y 2 = 1 − y 1 + f (y 1 ) + O( 2 ), where f (y 1 ) is quadratic in y 1 . In order to make a comparison to the Moran model, we ask that this line passes through the points y = (1, 0) and y = (0, 1), which implies that f (0) = 0 and f (1) = 0, which leads to the two conditions With this choice of the birth rates, the slow subspace takes the form It is interesting to note that the conditions in Eq. (18) have also eliminated any reference to the death rates δ i , and that as long as the competition rates do not satisfy Γ ≡ γ 11 + γ 22 − γ 12 − γ 21 = 0, the slow subspace will be curved. Simulations show that to an excellent approximation, the deterministic system collapses down to a line given by Eq. (19), as shown in Fig. 3. The agreement persists even if we do not impose the conditions (18), so that the line does not pass directly through the points y = (1, 0) and y = (0, 1).
The effective stochastic dynamics on the slow subspace is found by applying the same arguments as in the neutral case. A key aspect of the approximation is that the same form of the projection operator will be used when = 0 as was used when = 0. In previous applications [17,18] this was found to be a very good approximation for small values of s, and we will see that a similar conclusion applies in the current case. Therefore applying P ij given by Eq. (13) to Eq. (6) gives Eq. (15), but now with To test the validity of the fast-mode elimination procedure we compare the results for the probability of fixation and the time to fixation found from the reduced model to Gillespie simulations of the original IBM [28]. Both of these quantities can be found, either analytically or numerically, from the backward FPE corresponding to the reduced SDE (15) [29]. Fig. 4 shows the reduced model captures the properties of the full model extremely well.
From Eq. (20) we find that, in terms of the time variableτ , the reduced model is equivalent to the Moran model with selection strength s = (b 0 −d 0 )(γ 11 −γ 12 ) /b 0 , only if Γ = 0. Thus the condition of fixed population size, used for historical reasons and for mathematical convenience, does not give the same results as the more physically motivated SLVC model. This difference allows for fundamentally distinct dynamics. The discrepancy is expected to become more marked in more complex versions of the model. We hope to report on this and give further details elsewhere.
G.W.A.C. thanks the Faculty of Engineering and Physical Sciences, University of Manchester for funding through a Dean's Scholarship.