J Stat Phys (2018) 172:3–43 https://doi.org/10.1007/s10955-017-1900-1 Exploiting Fast-Variables to Understand Population Dynamics and Evolution 1 2 George W. A. Constable · Alan J. McKane Received: 14 July 2017 / Accepted: 11 October 2017 / Published online: 1 November 2017 © The Author(s) 2017. This article is an open access publication Abstract We describe a continuous-time modelling framework for biological population dynamics that accounts for demographic noise. In the spirit of the methodology used by statistical physicists, transitions between the states of the system are caused by individual events while the dynamics are described in terms of the time-evolution of a probability density function. In general, the application of the diffusion approximation still leaves a description that is quite complex. However, in many biological applications one or more of the processes happen slowly relative to the system’s other processes, and the dynamics can be approximated as occurring within a slow low-dimensional subspace. We review these time-scale separation arguments and analyse the more simple stochastic dynamics that result in a number of cases. We stress that it is important to retain the demographic noise derived in this way, and emphasise this point by showing that it can alter the direction of selection compared to the prediction made from an analysis of the corresponding deterministic model. Keywords Stochastic models · Population genetics · Time-scale separation · Effective models · Noise-induced selection · Population dynamics 1 Introduction In a series of articles on biological evolution published in the Journal of Statistical Physics, it is natural to ask what expertise and insights statistical physicists can bring to the study of evolution, and in what way might their approach to the subject differ from biologists. Alan J. McKane alan.mckane@manchester.ac.uk George W. A. Constable g.w.a.constable@bath.ac.uk Department of Mathematical Sciences, University of Bath, Bath BA2 7AY, UK Theoretical Physics Division, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK 123 4 G.W.A.Constable, A. J. McKane If the subject is the one that will largely interest us in this paper—the study of evolution within the framework of population genetics—these questions are more easily answered. This is because in a system containing a large, but ﬁnite, number of individuals with given genetic characteristics, genetic drift leads to a stochastic dynamics which has many of the features which allow the application of the ideas and techniques of non-equilibrium statistical physics. We will use this formalism, but in addition our approach will have parallels with the traditional methodology of theoretical physicists. Firstly, we will stress the fundamental nature of the microscopic description. That is, we will start with genes as basic constituents, which can be in different states according to their type (allele), location (if the individual carrying the gene is located on an island), the sex of the individual carrying the gene, etc. Secondly, since the microscopic description contains too much detail which is irrelevant at the macroscale (or in our case, the mesoscale, where some stochastic element is retained) we will derive a reduced or effective model which contains parameters which depend on the parameters of the microscopic description and so encapsulate the relevant aspects of the microscopic description. Thirdly, we will be interested in generic behaviour. By this we mean that we will attempt to formulate a microscopic description that does not have inbuilt assumptions that make the model more easily solvable. Instead we will try to formulate the model in such a way that it could be generalised to include many more effects, without changing its structure. In this philosophy, the simplifying assumptions are brought in during the process of obtaining the effective model, and should be clearly stated. Finally, although the whole basis of our work is mathematical, we will use intuition to explore the admissibility of the techniques we use outside of their regime of strict applicability, and check their correctness through the use of computer simulations. Although these ideas are familiar to the theoretical physics community, they tend to be utilised less in the biological sciences. For example, many biologists may use quite complex verbal arguments to gain insights. Conversely, our methodology may differ from that of many mathematical biologists, since rigorous justiﬁcation will not be a central feature of our approach. In addition, many mathematical biologists are focussed on the deterministic dynamics found at the macroscale. Nevertheless, we view the approach we will discuss here as able to form a bridge between the intuition gleaned by biologists and the more analytic investigations of mathematicians. In this way we hope that our methods prove of interest to a wider audience outside of the theoretical physics community. In a previous paper [43], we have reviewed the process of setting-up a description of this class of biological systems in terms of its basic constituents, and from this deriving the mesoscopic equations governing the dynamics which generalise what might be the more familiar macroscopic equations. In particular, in Ref. [43], we give formulae for writing down the form of the mesoscopic dynamics in terms of quantities which appear in the microscopic formulation. This essentially is the ﬁrst point of our methodology described above, and so while we will discuss it here, we will refer the reader to this earlier paper for more details. Instead we will focus on the second point above, namely obtaining an effective theory that is more amenable to analysis than the original. There may be several ways of reducing the complexity of the model, but here we will concentrate on one which is based on time-scale separation arguments. That is, we will seek to identify fast modes which die away relatively quickly, and slow modes which endure at long times. The dynamics of systems featuring such timescale separation are illustrated in Fig. 1. This is, of course, a well-known procedure, perhaps the most famous example being in hydrodynamics, where the microscopic molecular dynamics can be replaced by a macro- 123 Exploiting Fast-Variables to Understand Population Dynamics… 5 scopic dynamics with a few long-lasting variables. Although this dynamics is macroscopic, a mesoscopic extension can also be derived along the same lines [23]. In the theory of dynam- ical systems, the concept of a centre manifold (CM) is another manifestation of these ideas. During our discussion of this methodology, there will be several illustrations of the third and fourth points discussed above, namely the wish to use generic structures and the use of numerical simulations to check the precision of the approximations we utilise. As we have already mentioned, one of our aims in writing this article is to make the ideas and techniques available to a larger audience. To help to achieve this we will present an application of the method in a pedagogical manner in Sect. 2. We have chosen one of the simplest possible systems: haploid individuals on one of two islands of equal size which can migrate from one island to the other. It will be assumed that are only two possible alleles which are modelled by a Moran process. After this informal, and hopefully easily accessible, introduction to the method, we will describe its application to a number of models in Sect. 3. These include: a haploid Moran model on an arbitrary number of islands with selection and mutation; a stochastic Lotka–Volterra competition model with an arbitrary number of islands and a stochastic Lotka–Volterra competition model with an arbitrary number of species; a derivation of the Hardy–Weinberg approximation from ﬁrst principles; a model of epidemic spread on a network. In Sect. 4, we will illustrate the method in the slightly more technical case where noise-induced dynamics are present. We will see that noise-induced selection can cause selection for genotypes that are neutral in a deterministic setting, and that further, this noise induced selection can, under certain conditions, be strong enough to reverse the direction of deterministic selection. We will illustrate this behaviour with reference to Lotka–Volterra competition models, where we will see that this effect can help alleviate the dilemma of cooperation, and a model of transitions between sex-chromosome systems. Finally, in Sect. 5 we conclude with a discussion. 2 Pedagogical Example In this section we will explain as simply as possible how to apply the ideas discussed in the Introduction to a concrete example. The example we choose is a Moran model with migration. We ask how this can be reduced to an effective one-island model. 2.1 Setting-Up the Model We set up the model at the microscale, that is, at the level of haploid individuals who each carry an allele of type 1 or of type 2. The individuals can reside on one of two islands, both of which can only carry a ﬁxed number of individuals, which we denote by N. We therefore denote by n the number of individuals carrying allele 1 on island 1, by N − n the number 1 1 of individuals carrying allele 2 on island 1, by n the number of individuals carrying allele 1onisland2,and by N − n the number of individuals carrying allele 2 on island 2. So the state of the whole system is given by only two variables, which we can form into the two- dimensional vector n = (n , n ). We would like to reduce this description to one involving 1 2 only one variable, which gives the fraction of individuals in the system carrying allele 1. This would allow us to calculate, for example, the probability that allele 1 or allele 2 ﬁx, and the mean time to ﬁxation. There is not enough information about the birth, death and migration of individuals to model them in any other way than as random processes, so the model is speciﬁed by giving the probability per unit time that the system transitions from its current state, given by the 123 6 G.W.A.Constable, A. J. McKane 1.0 0.8 0.6 0.4 0.2 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.8 0.6 0 10 0.4 0.2 0.0 0 500 1000 1500 2000 Fig. 1 Illustration of four systems featuring timescale separation that can be analysed with the methods reviewed in this paper. Top left panel: Phase plot for a haploid Moran model with two alleles on two islands with strong migration (described in Sect. 2). The deterministic dynamics rapidly collapse to a slow subspace, indicated here by the blue dashed line. Top right panel: Deterministic trajectories (grey arrows) for a system similar to that in the top left panel but with three islands, addressed in Sect. 3.1. Again, the deterministic dynamics rapidly collapse to a one-dimensional subspace indicated by the blue dashed line. Bottom left panel: Genotype frequencies as a function of time for a population genetic model, described in Sect. 4.2. Stochastic trajectories (solid lines) initially rapidly relax on quasi-deterministic trajectories (inset) before reaching a one-dimensional slow subspace along which the system moves on a slower timescale. Bottom Right panel: A neutral three-species Lotka–Volterra model, addressed in Sect. 3.2.2. Stochastic trajectories (orange) rapidly collapse along quasi-deterministic trajectories onto a two-dimensional slow subspace (blue surface), about which they are conﬁned (Color ﬁgure online) vector n, to a new state n . We write these transition rates as T (n |n), with the initial state on the right and the ﬁnal state on the left (some authors use the reverse convention). The probability distribution function (pdf) of the system, p(n, t ), is then given by the master equation d p(n, t ) = T (n|n ) p(n , t ) − T (n |n) p(n, t ) . (1) dt n =n 2 Exploiting Fast-Variables to Understand Population Dynamics… 7 This is relatively easy to understand. The ﬁrst term on the right-hand side is made up of the probability of being in state n multiplied by the probability of being in that state and making a transition to state n. It therefore represents the probability of starting in state n and making a transition to state n. In the same way the second term on the right-hand side represents the probability of starting in state n and making a transition to state n . Their difference, summed over all states n , different to n, gives the rate of increase of p(n, t ) with time. The form we take for the T (n |n) depends on the model choice. Here we choose the Moran model, because it is simple: it amalgamates births and deaths and asks that birth, death and migration events happen in such a way that the population size of each island, N, is kept ﬁxed. These are not the most realistic assumptions, and we discuss ways to relax them later in the paper, but they have the merit that the number of model parameters is kept to a minimum. The method of constructing T (n |n) also seems a little more complicated, due to the requirement of keeping a ﬁxed number of individuals on each island. This is done as follows: (i) Pick an island (with probability 1/2, since the islands are identical) and then pick an individual randomly from that island. Allow the individual to reproduce by duplicating the individual. (ii) With probability m, the progeny migrates to the other island. In this case choose an individual on the other island at random to die. (iii) With probability (1 − m), the progeny remains on the same island. In this case choose an individual on the same island at random to die. It should be noticed that the processes of birth and death are inextricably linked and that they are assumed to happen at rate 1, this choice being possible through a choice of time units. On top of this process, migration is imposed with a probability of occurrence equal to m (0 ≤ m ≤ 1). Following these rules, if the model is neutral the transition rate from the state (n , n ) to 1 2 the state (n + 1, n ) is 1 2 1 n (N − n ) 1 n (N − n ) 1 1 2 1 T (n + 1, n |n , n ) = (1 − m) + m . (2) 1 2 1 2 2 N N 2 N N Similar expressions can be be found for T (n −1, n |n , n ) and for T (n , n ±1|n , n ). 1 2 1 2 1 2 1 2 However, we would like to include selection in the model. In this case we have to weight the choice of picking an allele by the relative ﬁtness of that allele on a particular island. Since we are aiming to be as simple as possible to illustrate the basic ideas, we will assume that this ﬁtness weighting is the same on both islands, though it is simple enough to relax this (1) condition. Therefore we will denote the ﬁtness weighting of allele 1 to be W (n) and of (2) allele 2 to be W (n). Then the four transition rates from state (n , n ) to the new state are 1 2 (1) 1 W (n)n (N − n ) 1 1 T (n + 1, n |n , n ) = (1 − m) 1 2 1 2 2 W (n) N (1) 1 W (n)n (N − n ) 2 1 + m , 2 W (n) N (2) 1 W (n)(N − n ) n 1 1 T (n − 1, n |n , n ) = (1 − m) 1 2 1 2 2 W (n) N (2) 1 W (n)(N − n ) n 2 1 + m , 2 W (n) N (1) 1 W (n)n (N − n ) 2 2 T (n , n + 1|n , n ) = (1 − m) 1 2 1 2 2 W (n) N 123 8 G.W.A.Constable, A. J. McKane (1) 1 W (n)n (N − n ) 1 2 + m , 2 W (n) N (2) 1 W (n)(N − n ) n 2 2 T (n , n − 1|n , n ) = (1 − m) 1 2 1 2 2 W (n) N (2) 1 W (n)(N − n ) n 1 2 + m , (3) 2 W (n) N (1) (2) where W (n) = W (n)n + W (n)(N − n ), i = 1, 2, is the ﬁtness of the individuals on i i i island i. Here superscripts are a label for the two different alleles, whereas subscripts are a label for the two different islands. For further background on how to arrive precisely at the forms given by Eq. (3) the reader is referred to previous discussions in the literature [1,6,43]. (1) (2) If W and W are independent of n, then the selection is known as frequency independent selection. This will be assumed in this pedagogical treatment, and we write (1) (1) 2 (2) (2) 2 W = 1 + α s + O(s ) and W = 1 + α s + O(s ),where s is a selection coefﬁ- (1) (2) cient and α ,α are constants. Since s is typically very small, we do not expect it will be 2 (1) (2) necessary to include the order s terms in the expressions for W and W . Equations (1)and (3) deﬁne the microscopic model, and once an initial condition p(n, 0) for the pdf has been given, specify the dynamics for all t. All other systems discussed in this paper will have a similar structure; the form of the transition rates will differ depending on the system, but in all cases the substitution of these rates into the master equation (1) will give the dynamics. As we discussed in the Introduction, the validity of our methods and approximations are checked via computer simulations, and these use the microscopic model deﬁned by Eqs. (1)and (3). The simulations use the Gillepsie algorithm [25,26]which is developed within the same formalism discussed in this section. However, the master equation is difﬁcult to study analytically. It is for this reason that we make the diffusion approximation, replacing the microscopic model with a mesoscopic version. The diffusion approximation was applied very early on in the development of population genetics [21] and is widely used [15]. We do insist, however, that it should be derived from an underlying microscopic model, since there are potentially many microscopic models that give the same mesoscopic model, and so simply deﬁning the model at the mesoscale can lead to ambiguities. The idea itself is very simple: if N is large enough, the ratios n /N, which are formally fractions, are assumed to be continuous, and denoted by x .Atthe same −1 −3 time the master equation is expanded in powers of N , and terms in N and higher are neglected. This process can be carried out directly, but formulae exist for the analogues of the transi- tion rates which appear in the mesoscopic equations [43]. To use these we need to introduce what are in effect stoichiometric vectors corresponding to the four “reactions” in Eq. (3). In other words we write the ﬁnal state, n , in terms of the initial state, n,as n = n + ν ,where μ = 1,..., 4 specify each of the four reactions. So for example, in the ﬁrst reaction of Eq. (3), n increases by 1 and n does not change, so ν = (1, 0). Similarly, ν = (−1, 0), ν = (0, 1) 1 2 1 2 3 and ν = (0, −1). These identiﬁcations allow us to use Eqs. (18) and (21) of Ref. [43]to show that the A and B functions which appear in the mesoscopic equations are A (x) =− m (x − x ) 1 1 2 (1) (2) 2 + α − α s [(1 − m)x (1 − x ) + mx (1 − x )] + O s , 1 1 2 2 123 Exploiting Fast-Variables to Understand Population Dynamics… 9 A (x) =− m (x − x ) 2 2 1 (1) (2) 2 + α − α s [(1 − m)x (1 − x ) + mx (1 − x )] + O s , 2 2 1 1 (4) and B (x) = x (1 − x ) + m (x − x )(2x − 1) + O(s), 11 1 1 1 2 1 B (x) = x (1 − x ) + m (x − x )(2x − 1) + O(s), (5) 22 2 2 2 1 2 with B = B = 0. These then specify the model, and the general dynamical equations 12 21 which allow us to ﬁnd the dynamics are either the Fokker-Planck equation (FPE) 2 2 ∂ P(x, t ) 1 ∂ 1 ∂ =− [ A (x) P(x, t )] + B (x) P(x, t ) , (6) i ij ∂t N ∂ x 2N ∂ x ∂ x i i j i =1 i, j =1 or the equivalent Ito¯ stochastic differential equation (SDE) dx 1 = A (x) + η (τ ), i = 1, 2, (7) i i dτ where τ = t /N is a rescaled time and η (τ ) is a Gaussian white noise with zero mean and with a correlator η (τ )η (τ ) = B (x)δ τ − τ , i, j = 1, 2. (8) i j ij As for the microscopic model, substitution of the speciﬁc forms given by Eqs. (4)and (5) into the generic forms of the FPE or SDE gives the behaviour of the mesoscopic model for all time, provided an initial condition P(x, 0) is given. For more details on derivation and meaning of these equations standard texts on the theory of stochastic processes [24,53]maybe consulted, or previous articles on the application of these ideas in a biological context [1,6,43]. We end this section with two general points. First, if we take the limit N →∞ of Eq. (7) we obtain the macroscopic model, which is deterministic, since the noise is eliminated by taking the limit. The dynamics of the deterministic model is then given by dx /dτ = A (x), i i i = 1, 2. Second, we keep terms of order s in A, but neglect them in B, since we are envisaging keeping terms of order s/N or 1/N in the FPE, but discarding terms of order 2 2 3 −1 s /N , s/N and 1/N . This essentially assumes that s ∼ N , although we will keep s and N to be independent variables throughout the paper. 2.2 First Stage of the Reduction Process: Identifying the Fast and Slow Variables Although the mesoscopic equations (6)and (7) are potentially more manageable than the differential-difference equation (1), they are still formidably difﬁcult to analyse—equation (6) is a partial differential equation in three variables, and in most of the other systems discussed later in this paper, the corresponding equation may have tens or even hundreds of variables. To ﬁnd a simpler, or reduced, form we want to eliminate the variables which decay away quickly, since they are not relevant to making predictions about the medium- to long-term behaviour. This subset of slow variables will form a slow-subspace (SS), so that instead of allowing the system to explore the whole space of variables, we only allow it to move within this subspace. 123 10 G. W. A. Constable, A. J. McKane In practice, instead of searching for a SS directly, we frequently search for a CM which is composed not of slow-variables, which hardly change with time, but of conserved variables that do not change with time at all. In population genetics, for example, neutral theories may contain conserved quantities, due to symmetries in the system (the different alleles behave in the same way), and the effects of selection can be added as perturbative corrections, given the extremely small size of selection coefﬁcients. The CM is found by looking for ﬁxed points of the macroscopic equation (the macroscopic limit of the SDE with no noise term present). This ﬁrst stage of the reduction therefore consists of the following steps: 1. Identify a CM, perhaps by setting some parameters to zero in order to increase the sym- metry of the deterministic equations (this could be the neutral limit of the deterministic dynamics, for example). 2. Find the Jacobian at the ﬁxed points that constitute the CM, and so ﬁnd the eigenvalues and eigenvectors of the Jacobian evaluated on the CM. 3. Form a projection operator from the eigenvectors found in 2, which is used to operate on quantities in the full mesoscopic system in order to eliminate the fast variables. 4. Use this projection operator, or use conservation laws, to ﬁnd the point where the system ﬁrst reaches the CM. This will be the new initial condition for the reduced system. We will now illustrate these four steps on the pedagogical example. (i) Setting s = 0in Eq. (4), we see that there is a line of ﬁxed points x = x .Thisisthe 1 2 CM. (ii) The Jacobian of the deterministic system dx /dτ = A (x), i = 1, 2, with s = 0, is i i −m/2 m/2 J = . (9) m/2 −m/2 This matrix has zero determinant and a trace equal to −m, which immediately gives {1} {2} its two eigenvalues as λ = 0and λ =−m. Two typical features we would expect are illustrated here: the number of zero eigenvalues equals the number of dimensions of the CM (since there is no dynamics at all on the CM—it is comprised only of ﬁxed points) and the non-zero eigenvalue has a real part which is negative (so that it can be identiﬁed as the fast mode which dies away quickly). In this case the non-zero eigenvalue is real, which is a reﬂection of the fact that the Jacobian, J , is symmetric. Another consequence of J being symmetric is that we would normally be required to ﬁnd right- and left-eigenvectors of J , but these coincide for symmetric matrices and are given by 1 1 1 1 {1} {2} v = √ , v = √ . (10) 1 −1 2 2 We would expect that the eigenvector corresponding to the zero eigenvalue would lie in {1} the CM, and indeed v lies on the line x = x . The normalisation of the eigenvectors 1 2 2 {μ} {ν} has been chosen so that they are orthonormal: v v = δ ,where μ, ν = 1, 2 μν i =1 i i and δ is the Kronecker delta. μν {1} {1} (iii) The projection operator is deﬁned by P = v v , constructed only from the eigen- ij i j vectors of the zero mode(s). To illustrate its use, we operate with it on a general vector {1} {2} of the full system given by φ = C v + C v ,where C and C are constants. Then i 1 2 1 2 i i {1} {2} P φ = C v , that is, the term involving the fast mode(s) in φ , C v ,has ij j 1 i 2 i =1 i i been wiped out using the orthogonality of the eigenvectors. 123 Exploiting Fast-Variables to Understand Population Dynamics… 11 IC (iv) If the point at which the system begins is x , then we would expect it to reach the CMIC 2 IC CM at x = P x , since only the slow (zero) modes will have survived ij i j =1 j by this time. Here i = 1, 2 and ‘IC’ and ‘CMIC’ stand for ‘initial condition’ and ‘centre-manifold initial condition’ respectively. Applying the projection operator one CMIC CMIC IC IC ﬁnds that x = x = (x + x )/2. Another way to obtain this result is to use 1 2 1 2 the conserved quantity which exists in this degenerate system. From Eq. (4), one sees that d(x + x )/dτ = 0when s = 0. Therefore, x + x is unchanged in time, and so 1 2 1 2 IC IC CMIC CMIC CMIC CMIC x + x = x + x = 2x or 2x . 1 2 1 2 1 2 These results will be used to construct the reduced model. As an illustration of the fourth general point made in the Introduction relating to intuition and the checking of approxima- IC CMIC tions through simulations, we note that the trajectory from x to x is stochastic, not CMIC deterministic. Nevertheless, we will use x as the initial condition for the reduced system on the CM, even though it has been deduced through a deterministic argument. We expect the IC deterministic dynamics to dominate the collapse from x to the CM under the condition that the rate of migration (which controls the collapse of the system to the CM) is much stronger than the rate of genetic drift (which causes deviations from the deterministic collapse to the CM). Since the rate of genetic drift grows linearly with the population size [see Eq. (6)], this −1 condition can be expressed m N . In addition, the projection of the IC onto the CM as described is only strictly true when trajectories to the CM are linear (as in this case) or when the initial condition is close to the CM (in which case the linear approximation is applicable). If the deterministic trajectories to the CM are non-linear, more sophisticated mathematical techniques may be necessary CMIC IC to calculate x when x liesfar fromthe CM [54]. We will also continue to use the eigenvectors of the neutral model when constructing the s = 0 reduced model. It would be possible to ﬁnd perturbative corrections to these in s, but we expect the effects to be sufﬁciently small as to be completely negligible. These are judgements made on the basis of intuition. Their validity will be examined through numerical simulations, and a comparison made of analytic results found on the basis of these assumptions with results obtained by simulations of the original model. Finally it is worth noting that, as addressed in point (iv) above, an alternative line of attack is possible in which one transforms into the fast-slow basis of the problem, removes the fast-variables and then transforms back into the original, biologically relevant, variables of the problem (for an illustration of such an approach, see [11]). For the system at hand, the fast-slow basis is straightforward to obtain and is given by (x − x ) and (x + x ).However 1 2 1 2 in more general problems such a basis may not be straightforward to obtain analytically (we will explore such a scenario in Sect. 4.2) while the projection method that we develop here will continue to yield insight. We therefore continue to explore the current pedagological example using the projection formalism. 2.3 Second Stage of the Reduction Process: Construction of the Reduced Model In the ﬁrst stage we worked with a version of the model which had a CM, but our real interest is in the actual, realistic, model which will typically only have a SS. This will not change materially the process of collapse described in Sect. 2.2: we still expect what is in effect a deterministic collapse onto the SS. However, we would like to be able to assume that the system would then stay in the subspace which is deﬁned by the slow variables. This is true (at least deterministically) when there is a CM, since the system ceases to move once it reaches the CM (because A = 0). But this is not true when there is a SS. We therefore demand that A 123 12 G. W. A. Constable, A. J. McKane has no component in the fast directions. These conditions give the equation for the SS. There will still be a (weak) dynamics in the direction of the slow variables. We will also ask that there is no noise in the fast directions, only in the slow directions. In this way, the system is effectively constrained to evolve in the SS. This second stage of the reduction therefore consists of the following steps: 1. Ask that A has no components in the fast directions. This gives the equation that the SS must have for this to be true. 2. Apply the projection operator to the SDE of the full system, to obtain the SDE of the reduced system. We can now illustrate these steps on the pedagogical example. {2} (i) The condition that A has no component in the fast direction is v · A = 0, where A is the full (s = 0) form given by Eq. (4). This condition is simply A (x) − A (x) = 0. 1 2 (0) (1) Substituting x = X + sX + O(s ) into this equation we can determine the SS. i i (0) (0) We know that when s = 0 the SS is the CM and that X = X ,however using 1 2 (1) (1) A (x) − A (x) we ﬁndthatitisalsotruethat X = X , so the SS is deﬁned by 1 2 1 2 x = x to order s. It therefore has exactly the same form as the CM, and is linear. This 1 2 is not true in general, and is only a feature of this simple pedagogical example. The variable x ,or x , will be denoted by z on the SS. 1 2 (ii) Applying the projection operator to the SDE (7)gives (2/ 2)(dz/dτ) for the left-hand side, since x = x ≡ z on the SS. The ﬁrst term on the right-hand side gives (2/ 2) A(z) 1 2 for a similar reason: A (x) = A (x) on the SS x = x , and we write A or A on the 1 2 1 2 1 2 SS in terms of z as A(z). Therefore the reduced SDE is given by dz 1 = A(z) + √ ζ(τ ), (11) dτ where (1) (2) 2 A(z) = α − α sz (1 − z) + O s , (12) and where 1 1 1 {1} {1} ζ(τ ) = √ √ v η (τ ) + v η (τ ) = √ [η (τ ) + η (τ )] . (13) 1 2 1 2 1 2 2 2 2 2 From the properties of the noise η(τ ), including the correlation function in Eq. (8), it follows that ζ(τ ) is a Gaussian white noise with zero mean and with a correlator ζ(τ )ζ(τ ) = [B + B ] δ τ − τ ≡ B(z)δ τ − τ . (14) 11 22 where B(z) = z (1 − z) + O (s) . (15) The reduction we have described here has produced an effective mesoscopic model deﬁned by Eqs. (11), (12), (14)and (15), which has only one variable, and which is sufﬁciently simple that it can be analysed mathematically. The pedagogical example was chosen on the grounds that its simplicity allowed the method to be clearly explained, and unfortunately this means that the reduction does not give that much useful information. For instance, the ﬁnal results do not depend on m, and all we have done is to verify a known result [42] that with this type 123 Exploiting Fast-Variables to Understand Population Dynamics… 13 of selection and with the same selection pressure on all islands, the population behaves in a similar way to a well-mixed population of size equal to that of the two islands added together. Actually, if one takes the calculation to higher order in s, then one can ﬁnd an exception to this: migration-selection balance can occur if selection acts in opposing directions on the different islands. In later sections we will describe applications of the method which will yield informative reduced models. Although the essentials of the reduction method will be the same as those that we have described in this section, a few aspects will make it seem slightly more complicated. Apart from the number of variables being greater, requiring additional indices, we mention (a) The generalisation of the Jacobian (9) will typically not be a symmetric matrix. This means that the non-zero eigenvalues will be complex in general, and the left- and right- {μ} {μ} eigenvectors will not coincide. We will use the notation u and v respectively for the {μ} left- and right-eigenvectors corresponding to the eigenvalue λ ,where μ = 1, 2,.... {μ} {ν} They will be chosen to be orthonormal, that is, u · v = δ . μν (b) As a consequence of this a typical term in the projection operator will involve left- and right-eigenvectors and the condition for there to be no deterministic dynamics in the {μ} {μ} {μ} μ-direction will be u · A(x) = 0, that is, will involve u rather than v . These are simply small technical details, and the method as discussed here is in essence that used in the other applications which we will now go on to discuss. However, having ¯ ¯ accounted for these points, we can deﬁne more general forms of A(z) and B(z) that will be relevant for later problems. In particular, for a system with initially M variables we ﬁnd an effective one-dimensional approximation of form Eq. (11) with; {1} A(z) = u A (x)| (16) i CM i =1 and {1} {1} B(z) = u u B (x)| , (17) ij CM i j i, j =1 {1} and where we recall that u is the left-eigenvector corresponding to the zero eigenvalue. {1} Since u is perpendicular to all of the fast directions [6], these two terms can be viewed as the deterministic and noisy components of the full problem respectively projected onto the {1} {1} SS; that is using P = v u . ij i j Finally, since one of the themes of this paper relates to the utilisation of techniques from theoretical physics to population genetics, and other areas of theoretical biology, we could import the use of the bra-ket notation from quantum mechanics [17]. In this notation the right- {ν} {μ} eigenvector v is written as the ket |ν and the left-eigenvector u as the bra μ|. Then the {1} {1} orthogonality relation becomes μ|ν= δ and the projection operator P = v u may μν ij i j be written as P =|1 1|. Though undoubtedly more elegant, for consistency with earlier work we shall not use the bra-ket notation here. 3 Applications of the Fast-Mode Elimination Procedure In this section we will discuss some applications of the formalism which we have described in Sect. 2, making reference to previous work, notable features of the various models and possible future work. 123 14 G. W. A. Constable, A. J. McKane 3.1 The D-Island Moran Model The modelling and analysis of migration effects in population genetics has always been chal- lenging since it involves a spatial aspect in an essential way. Historically, it was Wright [65] who ﬁrst studied migration models in population genetics, however he in fact did not assume a spatial structure, since migratory individuals were chosen from a global, well-mixed, popu- lation. The stepping stone model [36] was one of the ﬁrst which did have real spatial structure. It consisted of a line of islands, but where migration could only take place between an island and its neighbours on either side. If we view migration as an interaction event between two islands, this is a one-dimensional model with nearest-neighbour interactions. Expressed in this way, an obvious generalisation is to a network of islands with interaction strengths between the islands proportional to the probability of migration between the islands. Such a model was investigated by Nagylaki [45], although with discrete generations and strong migration, where the probability of a migration event is of the same order as a birth or a death. Assumptions either made the analysis difﬁcult to follow or were not thought to be widely applicable, and many further studies of this kind followed which made a different set of assumptions. These, and further discussions and analysis, can be found in the book by Rousset [57], while more recent results that utilise probability generating functions [52]are developed in [34]. In Sect. 2 a simple two island model was introduced. The D-island model is a generalisation of this but with added features. Details are given in earlier papers of ours [6,7,9], but examples of more general features are the migration probability to island i from island j, m ,which ij will in general not be symmetric (m = m ), and the fact that islands will be allowed to ij ji contain different number of individuals (β N on island i). The migration matrix is not so far deﬁned for i = j, however if one sets the probability of the chosen individual not migrating equal to m = 1 − m , then this deﬁnes m for all i, j. ii ji ij j =i The factor N which appears in the popoulation size is the only large parameter in the model, and the approximations made in Sect. 2 are reliant on this. So, for example, β should not be so small or large that we still cannot treat the factor β N as being of a similar magnitude to N. Similarly, the number of islands D should not be so large that it can be thought of as being of order N. It is likely that in some cases the approximations will continue to be good outside of their strict range of validity, but at present this can only be tested by comparing the analytic results with simulations. We will discuss the model with and without mutation separately, since typical questions which we are interested in answering differ. We begin with the model with no mutation. 3.1.1 The D-Island Model Without Mutation This model has been analysed in detail in Refs. [6,7], where further details may be found. Here we only note that, in addition to the points already made, that the probability of choosing an island on which an individual is then chosen to die or migrate, f , has to be done with a probability proportional to β , if we are to get sensible results. When the approximation described in Sect. 2 are made, one again ﬁnds that the reduced model has only one degree ¯ ¯ of freedom, and is given by Eqs. (11)and (14). Even the functional forms of A(z) and B(z) are unchanged, although now they contain parameters which are functions of most of the parameters of the starting model: ¯ ¯ A(z) = a sz (1 − z) + O s , B(z) = b z (1 − z) + O (s) , (18) 1 1 123 Exploiting Fast-Variables to Understand Population Dynamics… 15 Fig. 2 Plots for the probability of ﬁxation (left panel) and mean time to ﬁxation (right panel) as a function of the projected initial conditions for the D = 3 island Moran model. The solid lines are calculated from the reduced model and the various symbols indicate the results obtained from simulations. For each different colour/symbol different α vectors are used; green squares, α = (1, 1, −1), red triangles, α = (−1, −1, 1), and blue circles, α = (1, −2, −1). All other parameters are kept constant; s = 0.005, N = 200, β = (3, 2, 1) and the migration matrix m is ﬁxed, though not given here. Simulation results are the average of 5000 runs where D D m f m f {1} ij j {1} ij j a = u α , b = u . (19) 1 j 1 i i i β i, j =1 i, j =1 {1} Here α is the relative ﬁtness of the ﬁrst allele over the second on island i and u is the left-eigenvector corresponding to the zero eigenvalue. So even with the added complexity of each island differing in size, arbitrary migration probabilities, selective advantage of allele 1 over allele 2 varying from island to island, and the number of islands itself arbitrary, the reduction gives a standard (i.e. non-spatial model) Moran model with selection with effective parameters which can be seen from Eq. (19) to contain information from virtually all the parameters of the original model. It should be noted that in order to prove results pertaining to the spectrum of the eigenvalues of the Jacobian, we assumed that the migration matrix had a structure which in effect meant that no subgroup of islands were isolated from any other [6], and so the results displayed in Eq. (19) are subject to this restriction. This rules out, for example, the case where the islands may be divided into two subgroups, with no migration between the subgroups. As mentioned the reduced model is the mesoscopic version of the Moran (or, in fact, the Wright-Fisher model) with selection. The probability of either allele 1 or allele 2 ﬁxing for a given initial state and also the mean time to ﬁxation may be straightforwardly found [15], given as they are by the solution to second-order ordinary differential equations [7]. These are denoted by Q(z ) and T (z ) respectively. Here z is the initial condition for the reduced 0 0 0 CMIC system, which was referred to as x in the ﬁnal paragraph of Sect. 2.2 on the ﬁrst stage of the reduction process. We have changed notation to the new coordinate z and also indicated that it is an initial condition through use of the subscript 0, rather than the superscript CMIC. It is clear that both Q and T have to depend on precisely where the reduced system is initialised. In Fig. 2, ﬁxation times calculated from the reduced model and simulations of the under- lying microscopic model, from which it was derived, are shown. The very good agreement between the two gives us conﬁdence in the method, with more comparisons [6,7] also giving further support. The calculation can also be taken to next order in s, and an effective term of order s found in the expression for A(z) giventoﬁrstorder in s in Eq. (18). In this −1/2 case s is tentaively assumed to be of order N (although a direct indentiﬁcation is not −2 made), so that terms of order higher than N have been ignored in the expansion of the 123 16 G. W. A. Constable, A. J. McKane master equation. Novel effects can be investigated through use of the s term, and a greater range of parameters explored. We refer the reader to the original literature for a discussion of these [6,7]. 3.1.2 The D-Island Model with Mutation So far in this article we have examined the effects of migration, selection and genetic drift, but not that other important process of population genetics: mutation. The inclusion of mutation has a drastic effect on the long term behaviour of the system, since it is now possible in principle for one allele to mutate to another at any time, and so concepts such as ﬁxation probabilities and mean time to ﬁxation do not apply. On the other hand, the pdf of the allele frequencies is now non-trivial and becomes stationary at long times. This is an interesting quantity which characterises the model and we use it to test the accuracy of the effective model found through the reduction procedure. The microscopic model is constructed in a similar way to that described in Sect. 2.There is more than one way that mutation can be included; here we follow the procedure described in Ref. [9]. In this case we allow birth/death/migration events to happen a fraction ξ of the time, and mutation events a fraction (1 −ξ) of the time. The mutation rate from the ﬁrst allele (1) to the second on island i is denoted by κ and from the second to the ﬁrst on island i by (2) κ . It may be a little unusual to allow rates to vary from one island to the other, but allowing for this does not markedly increase the complexity of the calculation and so we include it. If we denote the rates without mutation (such as are shown in Eq. (3)for thecaseoftwo islands) as T (n |n), where the subscript S indicates that selection has been included, but not mutation, then the corresponding rates with mutation and selection included are (β N − n ) (1) i i T (n + 1|n ) = ξ T (n + 1|n ) + (1 − ξ ) κ , MS i i S i i β N (2) T (n − 1|n ) = ξ T (n − 1|n ) + (1 − ξ ) κ . (20) MS i i S i i β N Here the dependence of the transition rates, T (n |n),onelementsof n that do not change in the transition has been suppressed. We may rescale time in the master equation by a factor of ξ and absorb a factor of (1 − ξ)/ξ in the mutation rates. This effectively means that we can drop the factors ξ and (1 − ξ) from Eq. (20). Making the diffusion approximation, the mesoscopic model is given by Eq. (7) with the noise correlator given by Eq. (8). In this case (1) (1) (2) A (x)| = A (x)| + κ − κ + κ x , i i i MS S i i i (1) (2) B (x) = B (x) + O κ , κ , (21) ij ij MS S (1) (1) (2) (2) (1) (2) where κ = (κ ,...,κ ) and κ = (κ ,...,κ ). Since mutation has been modelled 1 D 1 D as a linear process, the κ dependence in A (x) in Eq. (21) is exact. We will neglect the κ dependence in B (x) for precisely the same reasons that we neglected the dependence on ij (1) (2) the selection coefﬁcients: we are assuming that the elements of κ and κ are so small that −1 they can be thought to be of the same order as N . We therefore only keep terms of order (1) (2) 2 κ /N, κ /N and 1/N in the FPE. The reduction process itself is similar to that discussed previously since, as just mentioned, mutation rates are generally very small, and so they can be treated as perturbations of the 123 Exploiting Fast-Variables to Understand Population Dynamics… 17 (a) (b) (c) (d) Fig. 3 The stationary pdf for the D-island Moran model on the slow subspace for a range of systems with various parameters which are omitted here for brevity but which can be found in Ref. [9]. The solid black line is obtained from an analysis of the reduced model, the orange histogram from simulations of the original microscopic model, and the dashed line from a well mixed model with the same total system size and average mutation rates (weighted by island size) neutral model in exactly the same way as was done for selection strengths. Therefore the reduced model is given by Eqs. (7)and (8) with B(z) unchanged from the form given in Eq. (18). Perhaps not surprisingly A(z) is modiﬁed by the addition of an extra term depending on the mutation rates [9]: (1) (1) (2) (1) (1) (2) ¯ ¯ A (z) = A (z) +ˆ κ − κ ˆ +ˆ κ z = a sz (1 − z) +ˆ κ − κ ˆ +ˆ κ z, MS S 1 (22) where D {1} (1) D {1} (2) u κ u κ (1) i i (2) i i κ ˆ = , κ ˆ = . (23) β β i i i =1 i =1 and where a is deﬁned in Eq. (19). Here, as in Sect. 3.1.1, we retain terms only to order s. The stationary pdf of the effective theory can be straightforwardly found from the FPE ¯ ¯ corresponding to the SDE (7). Using the explicit forms for A(z) and B(z) one ﬁnds that [9] c c 1 2 p (z) = N z (1 − z) exp (c z), (24) st 3 where N is a normalisation constant and N N N (1) (2) c = κ ˆ − 1, c = κ ˆ − 1, c = a s . (25) 1 2 3 1 b b b 1 1 1 A comparison between simulations of the original model and calculations from the reduced model in Fig. 3 shows that the reduced model captures well the features of the full model. There are several ways in which the work discussed here and in the previous section could be taken forward. There is scope for biologists to tailor the technique to their own interests, perhaps including additional processes with their own set of parameters and dropping others. Mathematicians may be able to provide conditions under which the approximations made would be expected to be valid, perhaps giving upper bounds for the negative real parts of the non-zero eigenvalues. Another extension is to perform a similar analysis, but starting not from the Moran model, but from one which is closer to those used in ecological modelling. We now go on to discuss this. 123 18 G. W. A. Constable, A. J. McKane 3.2 The Stochastic Lotka–Volterra Competition Model The theory of evolution has had a very convoluted history, and a reﬂection of this was the signiﬁcant contribution made to the subject by theoretical studies, at least compared to other areas of the biological sciences. This had repercussions on the nature of the mathematical models used in these studies: they tended to be unrelated to broader questions relating to the organism, and more focussed on the combinatorics of allele selection. This was a good strategy when trying to test the ideas of Darwinian evolution, but it tended to isolate the theoretical development of the subject from developments elsewhere. An example is the Wright-Fisher model [22,65], one of the ﬁrst models of genetic drift, as well as the precursor of the Moran model [44]. In this model there is no competition between individuals—a trait which is obviously a feature of Darwinian evolution. The population therefore grows very quickly, but is kept under control by sampling from the very large pool of individuals that come into existence, in order to form the next generation. Only a ﬁxed number, N, of individuals are retained to form each new generation (Wright-Fisher) or births and deaths are coupled so the population at any given time is always equal to N (Moran). This leads to an artiﬁciality in the way that the models are set-up. In this section we will use as a starting point models in which the population is regulated by competition, rather than by a ﬁctitious constraint which ﬁxes the population size. This is a closer reﬂection of reality, and indeed the formulation of aspects of the models seem less contrived. It is a constant theme in the biological modelling literature that models of evolution should have a more ecological ﬂavour, and this approach conforms to these views. An apparent disadvantage is that the number of variables is increased. To see that, recall that the Moran model of Sect. 3.1 had only D variables, since the number of individuals carrying the second allele on island i could be expressed in terms of the number of individuals carrying the ﬁrst allele on the same island i.e. N − n . If the population of an island is not ﬁxed, the number of individuals carrying the ﬁrst and second allele can independently vary, and so the number of variables will double. Below we will describe how application of the reduction methods to models with competition show that they reduce to Moran-like models in the medium-to-long term [8,10]. The competition will be chosen to be of the simple Lotka–Volterra type [56], but in principle more complex competitive processes could be utilised. Since these are stochastic models, we will refer to them as stochastic Lotka–Volterra competition (SLVC) models. 3.2.1 The D-Island SLVC Model (α) As indicated above this system has 2D variables which in the microscopic model are n , where i = 1,..., D labels the islands and α = 1, 2 labels the allele. The comment made above concerning SLVC models being less contrived can be illustrated here in the way that migration is modelled. The procedure for doing this in the Moran model involves consider- able care in making sure that there are no biases built into the way the transition rates are constructed [see Eq. (3) in the simple two-island case] while keeping the population of both islands ﬁxed. For example, one has to ensure that a death on island j occurs before a migra- tion from island i to this island is allowed. In the SLVC model, one simply speciﬁes birth, (α) (α) (αβ) death and competition rates, respectively b , d and c , which are all independent of i i i each other. For the neutral version of the model, the birth, death and competition rates are the same for all alleles (and are denoted by a superscript 0); selection is introduced through a small perturbation in ,where is the selection strength: (α) (0) (α) (α) (0) (α) (αβ) (0) (αβ) ˆ ˆ b = b 1 + b , d = d 1 + d , c = c 1 + cˆ . (26) i i i i i i i i i 123 Exploiting Fast-Variables to Understand Population Dynamics… 19 Fig. 4 Fixation probability of allele 1 (left panel) and mean unconditional time to ﬁxation (right panel) as a function of the projected initial condition z for an SLVC model with D = 4 islands in the neutral case (blue) and in the case with selection (red, = 0.05). Symbols: mean obtained from 3000 stochastic simulations of the microscopic system; lines: theoretical predictions for the ﬁxation probability and mean time to ﬁxation obtained from the reduced model. Here the parameter V is equal to 150 (Color ﬁgure online) The diffusion approximation is applied in the same way as in the Moran model, although now the large parameter is not N, which is no longer present in the deﬁnition of the model, but V , which is some measure of the size of the system, such as the volume. Although there are 2D variables initially, 2D − 1 of these are fast, and so the reduced model has again only one variable. It may be possible in some parameter regimes to see a clear cut decay ﬁrst to the D variables of a Moran-type model with ﬁxed populations on each island, and then a slower decay to an effective one-island model, which parallels the discussion in Sect. 3.1,but in many cases these time-scales will be similar or will overlap. The time-scales are related to the inverse of the eigenvalues of the Jacobian and, in general, these are complicated functions of all the model parameters. While it is true that the reduced SLVC D-island model does reduce to a system which has a mesoscopic description given by Eqs. (11)and (14)—although with N replaced by V —and ¯ ¯ with B(z) = b z(1 − z) to leading order in , the form of A(z) is a little different. It is found to be given by [49] A(z) = z (1 − z)(a − a z) + O . (27) 1 2 In this case a , a and b are functions of the rates given which appear in Eq. (26), β ,and 1 2 1 i the left-eigenvector of the Jacobian corresponding to the zero eigenvalue. This change in the form of A(z) may be slight, but it could give signiﬁcantly different ﬁxation probabilities and mean time to ﬁxation. One reason for this is that it is now possible to have a ﬁxed point in the deterministic dynamics. This dynamics will be given by Eq. (11) without the noise term, and ¯ ¯ so ﬁxed points are solutions of A(z) = 0. When A(z) has the structure shown in Eq. (27), an internal ﬁxed point (i.e. one not at the boundaries z = 0or z = 1) is possible if a = 0: z = a /a , where the asterisk denotes a ﬁxed point. It will only exist if 0 < a /a < 1, but 1 2 1 2 if it is stable, it may prolong the time taken for the system to ﬁx (reach the points z = 0or z = 1). Similarly if it is unstable, it may lead to a shorter mean time to ﬁxation. To test the approximation we again compare the ﬁxation probabilities and mean ﬁxation time derived from the reduced model and those found from simulations of the original model. Although the form of A(z) is slightly more complicated than before, it is nevertheless still straightforward to work with the ordinary differential equations for Q(z ) and T (z ).The 0 0 results shown in Fig. 4 indicate that the reduction method is again working well in this case. 123 20 G. W. A. Constable, A. J. McKane 3.2.2 The M-Allele SLVC Model So far we have only discussed individuals which are haploid and carry one of two possible alleles. Here we discuss the generalisation to M alleles. This is interesting for a number of reasons, not least because we are able to recognise patterns that are not apparent in the two (α) allele case. As in Sect. 3.2.1, the variables in the microscopic model are n , where now α = 1,..., M and there is no island label, because we will assume that the population is well- mixed. In the corresponding haploid multiallelic Moran model there are only M −1 variables (a) (α) n , a = 1,..., M − 1, since the ﬁxed population constraint n = N, means that α=1 M −1 (M ) (M ) (a) n can be expressed in terms of the other M − 1 variables: n = N − n .Here a=1 the Greek indices α, β, ... will always run from 1 to M and the Roman indices a, b,... will always run from 1 to M − 1. Previously we used the reduction method to obtain an effective model which was amenable to analysis. Here will have a different perspective: we will ask if we can perform a reduction on the multiallelic SLVC model with M variables to the multiallelic Moran model with M − 1 variables. If this is so, then the more natural SLVC model will give the same results as the Moran model at medium and long times. From a mathematical point of view, the difference between the reduction described in Sect. 3.1 and in Sect. 3.2.1 is that previously there were D −1or2D −1 fast modes, and a single slow mode, whereas here there is one fast mode and M − 1 slow modes, thus giving an effective model which is (M − 1) dimensional [10]. The reduction procedure has been carried out in Ref. [10]. An equation analogous to Eq. (11) is found, but now in (M − 1) variables: (a) dz 1 (a) (a) = A (z) + √ ζ (τ ), a = 1,..., M − 1, (28) dτ (a) where ζ (τ ) is a Gaussian noise with zero mean and with a correlator (0) (0) 2b c (a) (b) (a) (a) (b) ζ (τ )ζ (τ ) = z δ − z z δ τ − τ . (29) ab (0) (0) b − d Here, the notation of underlining a vector is used for (M − 1)-dimensional vectors and bold for M-dimensional vectors. The correlation function is given to zeroth order in ,which is the neutral model result. It is exactly the form found in the M-allele Moran model, up to some rescalings which are absorbed into the new time τ [10]. In the neutral case A = 0,and the SLVC model reduces exactly to the Moran model at medium to long times, after rescaling the time. If selection is included, A is no longer equal to zero. In order to aid the comparison with the Moran model, it is useful to introduce two quantities: (ab) (ab) (aM ) (Mb) (MM ) C ≡ˆ c −ˆ c −ˆ c +ˆ c , (30) and (0) (α) (0) (α) ˆ ˆ b b − d d (α) ≡ . (31) (0) (0) b − d 123 Exploiting Fast-Variables to Understand Population Dynamics… 21 Then A (z) takes the form (a) (a) (a) (aM ) (M ) (MM ) A (z) = z −ˆ c − −ˆ c M −1 (b) (bM ) (M ) (MM ) (b) − −ˆ c − −ˆ c z b=1 M −1 M −1 (ab) (b) (bc) (b) (c) 2 ˆ ˆ − C z + C z z + O( ). (32) b=1 b,c=1 For now, we note that, just as in Eq. (27), A is cubic in the components of z.Thisisan important point in mappings between reduced SLVC and Moran models, as we will now see. We begin our discussion with the M-allele Moran model in the case where the selection (α) is frequency independent, that is, when the weight functions W , analogous to those intro- (α) (α) duced in Eq. (3), are independent of n. Speciﬁcally we assume that W = 1 + ρ s,where (α) the ρ are constants. This is the case that we have examined so far in this paper. Then we ﬁnd, after making the diffusion approximation, a model given by Eqs. (28), (29)and (32), (ab) but only if C = 0for all a and b [10]. If this condition holds, the reduced SLVC model and the Moran model with frequrncy independent selection match, provided that we make (α) (α) (α M ) the identiﬁcation ρ = −ˆ c . We also need to match up the selection strength used in the SLVC model () to the one used in the Moran model (s). The relation between them is: (0) (0) (0) s = (b − d )/b . Although some care has to be taken with making the identiﬁcation between the two models [10], one can note that the function A in the Moran model with frequency dependent selection is quadratic, and Eq. (32) is cubic in general, so the condition (ab) C = 0 gives the possibility of a direct mapping between the two models. We have assumed that selection is frequency independent so far, since this is the usual supposition made by many population geneticists and historically was the standard assump- tion used. However this may simply be a theoretical prejudice, since if one wishes to allow (α) the ﬁtness weightings W to depend on the composition of the population, one has to devise a model for this dependence, and so frequency independence is the simplest and most convenient choice. In addition, there are hints from experimental investigations that even if there are attempts to suppress factors that might lead to frequency dependent selection, it still seems to emerge [41]. Therefore it seems important to devise a natural way of including frequency dependence in modelling selection. Fortunately, there does exist a methodology to do this. It is based on ideas from game theory, where each allele “plays” a game with every other allele in the population [47]. In the way we choose to implement this [10], the ﬁtness weightings are taken to have the form M −1 M −1 (b) (b) n n (α) (αb) (α M ) W (n) = 1 + s g + g 1 − , (33) N N b=1 b=1 (αβ) where g is the payoff to allele α from interacting with type β. We can now make the diffusion approximation, just as in the frequency independent (α) (α) case, but now using W (n) given by Eq. (33), rather than the n-independent form W = (α) (α) 1+ρ s.Clearly,thestructureof W (n) in Eq. (33) can potentially lead to more complicated (a) z dependence in A, and indeed A (z) is now found to be cubic and given by [10] 123 22 G. W. A. Constable, A. J. McKane ⎡ ⎤ M −1 M −1 M −1 (a) (aM ) (ab) (b) (bM ) (b) (bc) (b) (c) ⎣ ⎦ sz G + G z − G z − G z z , (34) b=1 b=1 b,c=1 which is of exactly the same form as that given by Eq. (32). To get the precise correspondence (ab) (ab) (ab) between the two models one must take G =−C for all a and b,where G has the (ab) same structure as is displayed for C in Eq. (30), namely (aβ) (aβ) (Mβ) (ab) (ab) (aM ) G ≡ g − g ; G ≡ G − G . (35) (α M ) (α) (α M ) In addition the identiﬁcation g = −ˆ c hastobemade[10]. The fact that it is (aM ) (ab) (αβ) only G and G , and not g alone, that appear in the expression for A is interesting, (aβ) since the quantity G can be interpreted as a relative ﬁtness, namely the payoff to allele a against an opponent β relative to the payoff to allele M against the same opponent. Similarly, (ab) G is a relative relative ﬁtness. Therefore, as one would expect, it is not the actual payoffs which are important, but their values relative to the other payoffs. In Sect. 3.2.1 we discussed how existence of an interior ﬁxed point, that is, one not on the boundaries, could lead to different ﬁxation probabilities and mean times to ﬁxation. To investigate the possible existence of such ﬁxed points in the frequency dependent M-allele (a) case, we set A (z),given by Eq.(34), to zero. Now summing this expression over a gives ⎧ ⎫ M −1 M −1 M −1 ⎨ ⎬ (a) (bM ) (b) (bc) (b) (c) 0 = 1 − z G z + G z z . (36) ⎩ ⎭ a=1 b=1 b,c=1 M −1 (a) If the ﬁxed point is not to be on the boundary, then z = 1 and so the second a=1 bracket in Eq. (36) must vanish. Substituting this condition into the expression (34), which is M −1 (aM ) (ab) (b) itself taken to be zero, gives the ﬁxed point equation to be G + G z = 0, since b=1 (a) z = 0 for internal ﬁxed points. Since this non-boundary ﬁxed point equation is linear, there can generically be at most one ﬁxed point. The position of this ﬁxed point can therefore easily be found, and a determination made as to whether it lies in the SS and is therefore admissible. (a) (M ) A similar analysis for the frequency independent case yields the condition ρ = ρ for (α) all a. However, if all the ρ are equal there is no selection, so in the case of frequency independent selection there are no interior ﬁxed points. The ﬁnding that the more realistic SLVC model reduces to the Moran model with frequency dependent selection is another reason to use frequency dependence in the modelling of selection in the Moran model. Although, as we have already remarked, the resulting Moran model is still (M − 1)-dimensional, and so difﬁcult to analyse, some progress can still be made in some cases [10]. In this way the SLVC model may, in effect, be analysed. An example of such a situation is shown in Fig. 5. 3.3 Diploid Moran Model with Sexual Reproduction: The Hardy–Weinberg Assumption from First Principles In the Moran model discussed in Sects. 2 and 3.1, individuals are haploid (they carry only a single allele) and reproduction occurs asexually (individuals simply duplicate themselves). While this is a relevant case for certain simple organisms, it is less so for many complex organisms which are diploid and reproduce sexually, such as animals. Suppose we want to model such a system, with diploid individuals and two possible alleles at a single locus, reproducing sexually with any other individual in the population (here, there are no sexes). A mechanistic approach might be to attempt to model the three 123 Exploiting Fast-Variables to Understand Population Dynamics… 23 Fig. 5 Plots of the unconditional mean time until the ﬁxation of a single allele/species and the probability of the ﬁxation of an allele/species for the Moran and SLVC models with frequency-independent selection in the case M = 6 alleles/species. In these plots all alleles in the Moran model are under one of two selection pressures, while in the SLVC model all species have differing parameters that combine to give two selection pressures, making the system mappable to the Moran model presented. Analytical results are only available for the probability of ﬁxation. Simulations results are the mean of 10 stochastic simulations of the Moran (0) and SLVC models. Parameters used are given Ref. [10], where the parameterization of x in terms of κ is also described (1) (1) possible genotypes in the population; here we will denote the homozygotes by A A and (2) (2) (1) (2) A A , while the heterozygotes will be denoted by A A . If we ﬁx the population size to be N, this leaves two free variables. As we have previously discussed, this makes obtaining analytic quantities in the model far more difﬁcult than in the asexual haploid case, where the system was described by a single variable. Classic Approach Classic studies in population genetics circumvented this complexity by building single variable models that implicitly exploited a separation of timescales. It was noticed early on in theoretical population genetics that if there were no ﬁtness differences between the genotypes (i.e. the system was neutral) the frequency of genotypes in such a diploid system would quickly relax to Hardy–Weinberg frequencies [31,63], where the number of each genotype could be described in terms of a single variable, the frequency of one of the alleles. In the terminology of our present paper, the system would quickly relax (1) (2) (1) (A ) (2) (A ) to a CM. Denoting the allele frequencies by x = n /(2N ) and x = n /2N = (1) (1) (1) (2) (1) (1) (A A ) (2) (A A ) (1 − x ), and the genotype frequencies by y = n /N, y = n /N, (2) (2) (3) (A A ) (1) (2) y = n /N = (1 − y − y ), this is given by [20] 2 2 (1) (1) (2) (1) (1) (3) (1) y = x , y = 2x 1 − x , y = 1 − x . (37) Rather than model the dynamics of the diploid population, the dynamics of the alleles were modelled with the assumption that they existed at Hardy–Weinberg frequencies. This was assumed to also hold when selection was sufﬁciently weak that the deviations from these “equilibrium” frequencies were not too great [44] (note the conceptual similarities with our (1) (1) approach). Genotypes AA are assumed to be under selective pressure A A , genotypes (1) (2) (2) (2) A A under (1 + sh) and genotypes A A under 1 [20]. Note then that choosing h > 1 corresponds to overdominance, while h ≤ 1 corresponds to underdominance [20]. The details are given in Appendix A, however upon applying the diffusion approximation one obtains an FPE (6)orSDE (7), with [20] (1) (1) (1) (1) (1) A x = sx 1 − x x + h 1 − 2x , (1) (1) (1) B x = 2x 1 − x . (38) 123 24 G. W. A. Constable, A. J. McKane Mechanistic Approach Whereas Eq. (38) was developed using an apriori assumption that the system lay at Hardy–Weinberg equilibrium, we may now use the methods detailed in Sect. 2 to formally obtain an approximation for the dynamics on the CM. We note that a similar approach was taken recently in [32] and that this separation of timescales has been long noted and exploited [18,19,62]. We begin by modelling the genotypes themselves; genotypes encounter each other at a rate proportional to their frequency in the population, (αβ) weighted by a joint probability W that genotype α successfully mates with genotype β. In this way we account for selection. In a similar fashion to Sect. 2.1, we assume selection is small and formalize this by setting (αβ) (αβ) W = 1 + sα . (39) Expanding the master equation for small s as in Sect. 2.1, we obtain a two dimensional description of the dynamics. We now seek to understand how this two-dimensional description is related to the classic description. Having set up the model, we can proceed to the next stage of the analysis; identifying the fast and slow variables, as described in Sect. 2.2. We can obtain a CM by setting selection equal to zero (s = 0). In this case the CM is given by Eq. (37). We can now linearise the {1} {1} system about this CM to obtain the Jacobian, and use this matrix to obtain u and v ,the left and right eigenvectors of the Jacobian corresponding to the zero eigenvalue. Finally, using Eqs. (16)and (17), we can obtain an effective description for the system dynamics in terms of (1) z = y on the CM. However, in order to make a comparison between this effective theory (1) and the classic theory, we must express the effective theory in terms of x , the frequency (1) (1) (1) 2 of A alleles [see Eq. (38)]. Since y = (x ) on the CM, we must therefore make (1) the transformation x = z. The full calculation is given in Appendix A. We note that while there are some subtle mathematical points that need to be attended to, these should not distract from the key points of the method. Our ﬁnal result is that we can approximate the dynamics of the mechanistic model by an SDE of type Eq. (7) with terms (1) (1) (1) (23) (33) (13) (22) (23) A(x ) = sx (1 − x ) α − α + α + 2α − 6α (33) (1) (12) (13) (22) (23) (33) (1) 2 + 3α x + 3 α − α − 2α + 3α − α (x ) (11) (12) (13) (22) (23) (33) (1) 3 + α − 4α + 2α + 4α − 4α + α (x ) , (40) B(z) = z(1 − z). (41) While the form of Eq. (40) appears a little complicated, we can in fact show that this (αβ) becomes identical to Eq. (38) under the assumption that each term α can be decomposed into the sum of contributions from each allele; (αβ) (α) (β) α = α + α , (42) (α) and that α take the precise values (1) (2) (3) α = 1,α = h,α = 0 . (43) Since these values give a precise mapping from Eqs. (40)and (41)toEq. (38), we can (αβ) now ask what consequences this has for W . This matrix now takes the form ⎛ ⎞ 1 + 2s 1 + s + sh 1 + s ⎝ ⎠ W = 1 + s + hs 1 + 2sh 1 + sh . (44) 1 + ss 1 + hs 1 123 Exploiting Fast-Variables to Understand Population Dynamics… 25 Fig. 6 Left panel: phase diagram of the dynamics for the mechanistic diploid Moran model described in Sect. 3.3. Non-neutral dynamics are depicted by grey arrows. The blue dashed line shows the form of the neutral CM described in Eq. (37). The orange arrow shows a single deterministic trajectory in the neutral system. Since the population size is ﬁxed such that y + y < N, the area above the dashed black line is 1 2 outside the dynamical region. Right panel: ﬁxation probabilities as a function of z , the initial condition on the CM, for the diploid Moran model. In solid colours the ﬁtness of genotypes is parameterised using Eq. (39). (11) (13) (23) (12) (22) The dashed line gives the dynamics for the parameterisation α = α = α = 1, α = α =−1 (11) and α = 2. In all cases, N = 100 This is in fact exactly the form of W that we would expect at leading order in s if (αβ) (α) (β) (1) (2) (3) W = W W , W = 1 + s , W = 1 + sh , W = 1 , (45) i.e. if the ﬁtness of each genotype is the same as that postulated in Eq. (38) and the success of genotype pairings is a multiplicative function of the individual genotype ﬁtness. A similar mapping exists if we assume an additive interaction of the genotypes ﬁtness on their pair sexual success or an averaging effect (see Appendix A). It is testament to the great physical intuition of the founders of population genetics that this captures their assumptions. While we have essentially recovered here a known result, it is worth noting that of course the approach we have taken is not without merits. In particular, it allows us to explore a far richer ﬁtness landscape than the original model (see Fig. 6). 3.4 Stochastic Epidemics on Networks The networks which we have discussed so far in this paper have described the way in which islands interact with each other through migration. The islands have had populations of individuals which are already themselves large (in the sense that they are of order N,the large parameter in the system). These are metapopulation models, since the whole system is a population of populations [30]. However, in many network models the nodes are composed of a single individual, rather than a population. It is not immediately obvious how to apply the diffusion approximation in this case, and hence the reduction method of Sect. 2.However there are situations when the programme we have outlined so far can be pursued, and we now describe such a case. The example is taken from stochastic models of infectious diseases, rather than population genetics, but it illustrates the points we wish to make clearly. We will follow Ref. [48], where more details can be found. The model is an SIR model, which means that the individuals are either susceptible to the disease (in which case they are in class S), infected with the disease (in which case they are in class I ), or recovered from having the disease (in which 123 26 G. W. A. Constable, A. J. McKane case they are in class R). We assume that the disease is such that individuals recover, and do not die as a consequence of having the disease. Our interest will be in the properties of an epidemic which might occur, and we assume that this takes place over a much shorter timeframe than the demographic processes of birth and death. So the only processes which occur in the model are (i) infection, and (ii) recovery. The network structure enters because we assume that each individual has a ﬁxed number of contacts from which the infection is acquired, even though the contacts themselves will change. Therefore we may imagine all individuals located at the nodes of a network, where the degree of that node is equal to the number of contacts characteristic of that individual. The network is considered in the so-called dynamic limit, in which the network structure is assumed to evolve much more quickly than the epidemic, so that the only role of the network is to encode the number of connections a given individual has to other individuals. The variables at the microscale are therefore S , I and R where k labels the degree of the node k k k on which individuals of the type S, I and R are located. As is often done in SIR models, we assume that the population is closed, so that at time tS (t ) + I (t ) + R (t ) = N ,where N is k k k k k independent of time. This implies that one set of individuals—for example those which have recovered—can be removed: R (t ) = N − S (t ) − I (t ). The large parameter in the model k k k k is the total number of individuals: N = N ,where K is the maximum degree of the k=1 network. The speciﬁc network of interest in Ref. [48] was a truncated Zipf distribution where the probability of an individual having degree k is given by d ∝ k , with −3 <α < −2 and k = 1,..., K , but the method we describe is also applicable to other distributions. We have stressed throughout this article that one needs to start with the microscopic description, at the level of single individuals, to avoid ambiguities. However here, to avoid too much formalism, we will move directly to the mesoscopic model which can be derived from the microscopic description [48]: K (k) ds η (τ ) =−βks li + √ , (46) k l dτ l=1 K (k) di η (τ ) = βks li − γ i + , (47) k l k dτ l=1 where (k) (l) (k) η (τ )η (τ ) = B (x)δ δ(τ − τ ), (48) kl μ ν μν and where k = 1,..., K and x is a vector of all the 2K variables. Here β and γ are the infection and recovery rates respectively, and N is assumed to be sufﬁciently large that s = S /N and i = I /N can be assumed to be continuous. We will not give the precise k k k k (k) form of B (x), μ, ν = 1, 2, here, nor the form of the rescaled time τ . μν The problem we face is clear from Eqs. (46)and (47): this is a stochastic system with 2K variables, where in many cases of interest K will be not be small. This is then another case where we wish to reduce the number of variables, and where the reduction method described above may be useful. In fact, one can effect a partial reduction before applying the technique of Sect. 2 by making the ansatz s (τ ) = d θ(τ ) , which replaces the K equations (46)by k k the single equation dθ 1 =−βθ li + √ ξ(τ ), (49) dτ l=1 123 Exploiting Fast-Variables to Understand Population Dynamics… 27 Fig. 7 The distribution of epidemic sizes for N = 20000, K = 1000 and α =−2.5. The histogram gives the result for the full model, the solid blue line is the result for the reduced mesoscopic model and the red dashed line the result for the semi-deterministic mesoscopic model given by Eq. (50) (Color ﬁgure online) (k) where ξ(τ ) is a Gaussian noise with zero mean which is related to η (τ ). After this initial manoeuvre we attempt to further reduce this (K +1)-dimensional system by searching for fast and slow modes. Examination of the Jacobian [48] reveals that there is one K -fold degenerate eigenvalue which may be signiﬁcantly greater in magnitude than the remaining eigenvalue. If we denote the ratio of the magnitude of the former to the magnitude of the latter by , then as long as is small there will be effectively K fast modes and 1 slow mode, so we would expect to be able to reduce the motion to a two-dimensional SS where one of the variables is θ and the other is the slow one just mentioned. It is found that this additional variable is λ(τ ) = li (τ ) [48], so that Eq. (49) simply becomes l=1 dθ/dτ =−βθλ + ξ(τ )/ N. The equation for λ(τ ) is more complicated, involving as it does three different noise terms, and we do not give it here. Although the system has been reduced to a manageable two dimensional model, it is always worthwhile to perform numerical investigations to see if it is possible to make further reductions, since it is still the case that two dimensional stochastic systems are difﬁcult to study analytically. In this case, it is found that typically the noise term in the theta equation is much smaller in magnitude to the noise terms in the lambda equation. Therefore we can try to omit the noise term in Eq. (49) and combine the three noises in the λ equation together to obtain the further simpliﬁcation: dθ =−βθλ, dτ dλ 1 = λ (βφ(θ) − γ ) + √ σ ¯ ζ (τ ), (50) dτ 2 k where φ(θ) = k d θ , ζ(τ ) is a Gaussian noise with zero mean and ζ(τ )ζ(τ )= k=1 δ(τ − τ ),and σ ¯ is a function of θ and λ which is not given here [48]. The model deﬁned by Eq. (50) is a semi-deterministic mesoscopic model, in the sense that one of the dynamical equations is deterministic, having no noise term. In fact it can be identiﬁed as the Cox- Ingersoll-Ross (CIR) model [14], for which some analytic results are known. As a test of the various reductions that have been made on this model we will investigate how well they capture the distribution of epidemic sizes as given by r , the number of recovered individuals at the end of the epidemic [48]. The results are shown for a particular case in Fig. 7, and indicate that even the very much simpliﬁed model given by Eq. (50)gives a good representation of the result. 123 28 G. W. A. Constable, A. J. McKane 4 Revealing Noise-Induced Selection via Fast Variable Elimination In Sect. 2, we laid the groundwork for conducting fast-variable elimination in stochastic systems, while in Sect. 3 the utility of this approach was demonstrated. We have essentially shown that, given a separation of timescales exists, the dynamics of a high dimensional system can be approximated by the projection of these dynamics onto an SS. Thus far, this is true for both the deterministic and stochastic component of the dynamics [see Eqs. (16)and (17)]. However, there is a further consideration that we have until now omitted; noise-induced dynamics in the slow subspace. While accounting for such dynamics is more technically involved than the method outlined in Sect. 2, the resultant dynamical behaviour can be very interesting, and so it is worth describing the key aspects here. Noise-induced dynamics are the phenomena whereby the projection of the deterministic dynamics onto the slow-subspace [see Eq. (16)] does not entirely describe the time-evolution of the average dynamics when demographic noise is accounted for. In the context of a system that features a CM, this means that rather than there being no average dynamics along the CM, a bias in fact emerges in some direction. Although this is a second order effect (it scales inversely with the population size and disappears in the inﬁnite population size limit), it can in fact completely govern the dynamics along the CM, where the ﬁrst order terms that control the collapse of the system to the CM disappear. If the system under consideration is one featuring competing organisms, this bias can be interpreted as noise-induced selection (or drift-induced selection in a population-genetics context). The origin of these bias terms can be interpreted in various ways. First, and perhaps most intuitively, noise-induced dynamics can be graphically understood as resulting from a bias in how ﬂuctuations taking the system off the CM (or SS) return to the CM. Under certain scenarios, such as when the CM is strongly curved or the trajectories to the CM are divergent, ﬂuctuations off the CM do not return, on average, to the point on the CM from which they originated (see Fig. 8). This introduces a bias that stochastically ‘ratchets’ the dynamics in a preferred direction. Second, these noise-induced dynamics can be understood as a mathematical consequence of making a non-linear change of variables into the slow-subspace of the stochastic system. As we have mentioned, the SDEs that we work with should be interpreted in the Itos ¯ ense. In this context, different rules of stochastic calculus apply, and accounting for the additional terms that arise in the analysis gives rise to additional terms in the reduced description on the CM. For readers not familiar with the analysis of SDEs, Ito¯ calculus can in some sense be loosely understood as arising in the same way as Jensen’s inequality; the average of a function of a random variable is not necessarily the same as the function evaluated at the average of that random variable. Finally, in a more biological context, noise-induced dynamics in systems of compet- ing organisms can be understood resulting from a selective pressure to reduce variance in reproductive output (see Gillespie’s Criterion [27,29]). However, in contrast to Gillespie’s original study, variance between the reproductive output of organisms can arise as a result of the dynamics of the organisms, rather than being assumed apriori. Of course, the existence of these noise induced dynamics is not in itself dependent on a system exhibiting a separation of timescales. We have already mentioned Gillespie’s Crite- rion, which demonstrates that selection for a genotype with a lower variance in reproductive output takes place in a model with only a single variable (measuring the frequency of one of two genotypes). Further, in a similar single-variable population genetics context that assumes a well-mixed population of ﬁnite size, it has been shown that noise-induced selection favours 123 Exploiting Fast-Variables to Understand Population Dynamics… 29 genotypes that increase the population size [33,35]. However, while fast-variable elimination does not generate noise-induced dynamics itself, it does often give us a means of quantifying and analytically understanding the effect. As in the earlier sections, rather than transforming into the fast-slow basis of the problem, removing the fast variables and then transforming back into the original variables, we take a shortcut by using a non-linear projection of the dynamics onto the CM. The reason for this is two-fold. Firstly, it is more straightforward practically; it is in the original, biologically relevant variables that we can better understand the behaviour of the system. Secondly, as has long been recognised [13,55], the non-linear transform that yields the slow-fast basis is not always obvious. In order to obtain the non-linear projection, second order perturbation techniques can be used to deduce the effective dynamics. The full calculation is too long to reproduce here, however a clear and coherent explanation is given in Ref. [51]. There it is shown that if a system has M variables and a one-dimensional slow-subspace, then the equation for the reduced/effective dynamics can be expressed by dz 1 ¯ ¯ = A(z) + A (z) + ζ(t). (51) dτ N The term A(z) and the white noise term ζ(t ) retain the forms given in Eqs. (16)–(17) (that is, they are respectively components of the deterministic dynamics and noise projected onto the SS). Most interesting is the term A (z) as it is this that controls the noise induced dynamics. It is of order 1/N (induced as it is by demographic noise) and takes the explicit form {1} M {1} M {1} du du du i {1} {1} {1} {1} {1} k A (z) = u + u − 2u u v j i i j k N dz dz dz i, j =1 k=1 M {1} dv {1} {1} k {1} −u u u B (x) , (52) ij i j k CM dz k=1 {1} {1} where we recall that v and u are the right and left eigenvectors of the neutral Jacobian on the CM corresponding to the zero eigenvalue (see Sect. 2). The terms which feature {1} du /dz capture how quickly the fast directions change as a function of z, and can thus be understood as the component of the noise-induced dynamics that arises from the non- {1} linearity of trajectories to the CM. Meanwhile the term dv /dz is the component of the noise-induced dynamics that arises from curvature of the CM itself [51]. Note that if both {1} {1} v and u are independent of z (that is, the CM is linear and the direction of fast trajectories k k to the SS do not vary along the SS to ﬁrst order in the selection strength), then there can be no noise-induced dynamics at this order. This is the case in the models discussed in Sects. 2 and 3.1. Finally, there are also cases where, despite the CM being curved or the trajectories to the CM being divergent, there is still no noise induced selection as the terms in Eq. (52) all cancel. This is the case in both Sects. 3.2 and 3.3 for the parameterisations given in those sections. In what follows we will illustrate how the effective description provided by Eq. (51) can be used to analytically tackle some particular problems of interest. In Sect. 4.1 we will address a two-species Lotka–Volterra model. Unlike in Sect. 3.2 however, we will allow the species to have distinct birth, death and competition rates at leading order. Calculation of the effective system will reveal both that slow-living species and species that increase the global carrying capacity are stochastically selected for. In Sect. 4.2 we will describe a population genetic 123 30 G. W. A. Constable, A. J. McKane model of transitions between modes of sex determination. Although this population genetic model is much more complicated than the haploid Moran model, the presence of a CM in the neutral limit allows us to analytically characterise a noise-induced bias favouring the substitution of dominant neutral mutations. 4.1 Two-Species Lotka–Volterra Type Models We begin with a two-species Lotka–Volterra model of a similar type to that described in Sect. 3.2.2 [see Eq. (26)], except we now make the restriction (1) (1) (11) (21) (1) b − d = b(1 + ) , c = c ≡ c , (2) (2) (21) (22) (2) b − d = b , c = c ≡ c . (53) By taking the limit of large system size, we can apply the diffusion approximation, as described in Sect. 2.1. The system is now approximated by an SDE of the same form as Eq. (7). However, since the system is in two variables, it is difﬁcult to analyse. We next follow the approach taken in Sect. 2.2 and identify a CM. A CM exists under the above parameterisation if is equal to zero. It now takes the form (2) (1) (1) x = b − c x . (54) (2) (1) (1) (2) (2) ˜ ˜ Notice that in isolation, species 1 exists at x = b/c and species 2 at x = b/c . (2) (1) Further, if we assume that c > c , then increasing the frequency of species 1 in the population increases the joint carrying capacity of the species, as can be seen in Fig. 8.If additionally > 0, such that species 1 reproduces at a lower rate that species 2, then this can be interpreted as a model of public good production. Species 1 pays a cost to its reproductive rate to increase the carrying capacity of both species. Species 2 pays no cost, but still enjoys a reduced death rate in the presence of species 1 as a result of the lower competition parameter of species 1 (interpreted here as resulting from the production of a public good). However, in isolation, species 2 exists at lower numbers than species 1, interpreted here as resulting from the non-production of the public good. We can now brieﬂy describe the deterministic dynamics. In the neutral limit, where = 0, the population grows to a point on the CM described by Eq. (54). We now move away from the neutral limit, such that 1 > 0; the system grows to a point on the SS (which is equal to the CM at leading order in ) after which the system moves along the SS until species 1 is driven to extinction. We now use fast-variable elimination to characterise the dynamics when demographic noise is accounted for. As discussed in Sect. 2.2, our ﬁrst task is to identify the left and right eigenvectors of the system evaluated on the CM, Eq. (54). Note that while the right-eigenvectors are not so important in Sect. 2.3 (where there was no noise-induced selection), they become very important here, featuring as they do in Eq. (52). We ﬁnd (1) (1) ˜ ˜ 1 1 b − c z b − c z {1} {2} u = , u =− . (55) (2) (2) (1) (2) ˜ ˜ −c z bc /c + c z b b (1) 1 c z {1} {2} v = , v =− . (56) (1) (2) (1) (2) (1) −c /c ˜ (b − c z)/c b − c z 123 Exploiting Fast-Variables to Understand Population Dynamics… 31 Fig. 8 Phase plots for the deterministic dynamics of the neutral Lotka–Volterra model described in Sect. 4.1 under three different parameter scenarios. In each the CM is plotted as a blue dashed line. The surface of the orange ellipses is indicative of the distribution of ﬂuctuations arising from the centre of the ellipse. Arrows from the centre of the ellipse (going up, down, left and right) show possible ﬂuctuations away from the CM that occur with equal probability. Arrows travelling from the end of these ﬂuctuations back to the CM show the trajectories along which these ﬂuctuations are quenched. Bias in the direction to which ﬂuctuations are (1) (2) mapped are illustrated by green and purple lines. Top left: System with two species identical (b = b , (1) (2) (1) (2) d = d , c = c ). Fluctuations down and right are projected back to the CM in the direction of (1) x , but this is perfectly countered by ﬂuctuations up and left which are projected back to the CM with in (2) the direction of x . Top right: System with two species, the second reproducing and dying at a faster rate (1) (1) (2) (2) (2) (1) (2) (1) (1) (2) (b − d = b − d b > b , d > d , c = c ). As a result of this asymmetry, species 2 (2) experiences greater demographic ﬂuctuations (the ellipse is larger in direction x ). Consequently, ﬂuctuations (1) off the CM are projected back to the CM with a bias that favours x . Bottom: System with two species, the (1) (2) (1) (2) (2) (1) ﬁrst increasing the carrying capacity of the system (b = b , d = d , c > c ). The asymmetry in (1) the angle of the CM now induces a bias in the projection of ﬂuctuations back to the CM that favours x We can use these expressions for the left and right eigenvectors to obtain an approximate ¯ ¯ description of the dynamics in the SS using Eq. (51) with expressions for A(z), A (z) and (1) B(z) directly from Eqs. (16), (52)and (17), where z is the value of x on the SS; 123 32 G. W. A. Constable, A. J. McKane (1) ¯ ˜ A(z) =−z b − c z , (57) 1 2 S (1) (2) (2) (1) (1) ¯ ˜ ˜ ˜ A (z) = z b − c z c (b + d ) − c (b + d ) , (58) N ˜ (1) (2) (2) (1) (1) (1) ¯ ˜ ˜ ˜ ˜ B(z) = z b − c z z c (b + d ) − c (b + d ) + b(d + β) . (59) An analysis of these terms reveals a number of points. Firstly, and as we might expect, in (1) (2) (0) (1) (2) (0) (1) (2) (0) the limit c = c ≡ c , b = b ≡ b and d = d ≡ d the noise induced term A (z) vanishes and we are left with a system that takes the same form as that in Sect. 3.2. Secondly, we examine the limit in which = 0. Now A(z) = 0 as there are no dynamics along the CM in the inﬁnite population size limit. However, if the population has a ﬁnite size the noise induced term A (z) is not zero. Thirdly, continuing with the system when (1) (2) (0) = 0, and now additionally asking that c = c ≡ c (i.e. that the carrying capacity is unchanged by the composition of the population), we ﬁnd that A (z) is positive for all z along theCMif b < b . The species with the lower birth rate (and death rate, since b is 1 2 ﬁxed) is therefore selected for. This insight, made in Refs. [5,38,50], is a result of the fact that phenotypes which are reproducing and dying more quickly are subject to greater population ﬂuctuations (they have a larger rate of population turnover). Consequently, it is easier for the longer lived phenotype (lower birth/death rates), to invade and ﬁxate. This phenomena can be viewed as analogous to Gillespie’s Criterion. (1) (2) (0) (1) (2) (0) Turning instead to the limit = 0, b = b ≡ b and d = d ≡ d reveals a similar, but biologically distinct insight. In this case though the rate of population turnover of both species is identical, one species exists at greater numbers in isolation than the other (α) S ¯ ¯ (this species has a lower value of c ). Again while A(z) = 0, the noise-induced term A (z) is non-zero and drives the dynamics in the ﬁnite system. We now ﬁnd that A(z)> 0for all z (2) (1) on the CM if c > c . That is, the species that increases the joint carrying capacity of the species is selected for. One interpretation of this noise-induced term is that it is easier for a novel mutant species to invade a small population than a large one; thus species that increase the carrying capacity of the population receive a beneﬁt by being more stochastically robust to invasion attempts [12]. Note that in the above context, if > 0but N is ﬁnite, it is possible that A (z)> A(z) along the length of the SS. Biologically, this provides a mechanism that allows for the evolution of public good producing behaviour despite the evolution of such behaviours being forbidden in the deterministic limit. This has been noted in more typical population genetic models [33,35] as well as models of the form described here. If the population is not well mixed but exists in space, it has been shown that for weak migration the noise induced selection for public good production is ampliﬁed both in metapopulation [5,12]and continuous space models [28]. In [12] it was also shown that this behaviour (whereby a species that increases the joint population carrying capacity is stochastically selected for) is generic and robust to the inclusion of a suite of environmental variables in the model that can modify population size. 4.2 Models of Transitions Between Male and Female Heterogamety In this model, a diploid population is considered in which sex is genetically determined by a dominant mutation at a single locus. In mammals, sex is determined by the dominant Y chro- mosome, so that XY individuals are male while XX individuals are female. Such a system is termed male heterogamety. In birds the situation is somewhat reversed. Their sex determi- 123 Exploiting Fast-Variables to Understand Population Dynamics… 33 nation system features a dominant feminising W chromosome, such that ZW individuals are female, while ZZ individuals are male. This scenario is termed female heterogamety. Intriguingly, while these systems are relatively static in both mammals and birds, transi- tions between male and female heterogamety can occur in reptiles and amphibians. In this section, we will discuss how fast-variable elimination can be exploited to understand the impact of genetic drift on these transitions. We consider a system of male heterogamety comprised of XX females and XY males. A mutation arises on the X chromosome, changing it to an X and rendering it dominant to the Y, such that X Y individuals are female. Genotypes X X can also be produced, which are female, as can YY genotypes, which are male. This renders a system of ﬁve genotypes. Along with the absorbing state in which the system is entirely XX and XY (male heterogamety), there is also and absorbing state when the system is entirely X Y and YY (female heterogamety). Note that this state is analogous to that found in birds up to some relabelling of the chromosomes. We wish to understand the transitions between these states. We can construct a population genetic model of this process in a very similar way to the diploid Moran model (see Appendix A), except with matings restricted to being between males and females. We assume a ﬁxed population size N. Males and females of each geno- type encounter each other proportional to their frequency in the population. They produce a progeny, that inherits its chromosomes (alternatively, alleles at a single locus) in a Mendel- lian fashion from its parents. Simultaneously, in order that the population size is ﬁxed to N, another individual is picked to die. In order to account for selection against certain geno- types, the probability that each genotype is selected to die can be a weighted probability. The detailed model set-up is given in [61], however here we will discuss only the main results. Let the frequency of XX, X Xand X Y (female) and XY and YY genotypes (male) (1) (2) (3) (4) (5) be respectively given by x = (x , x , x , x , x ). Due to the condition of ﬁxed (5) (5) (1) (2) (3) (4) population size, we can then express x as x = 1 − x − x − x − x .Ifall genotypes are equally ﬁt, then this four-dimensional system exhibits a one-dimensional CM that connects the absorbing states, characterised in [3]. It is given by (1 − z) z(1 − z) z (1) (2) (3) x = , x = , x = , 2 2 2(1 + z) (1 + z) 1 + z 1 z (4) (5) x = (1 − z), x = , (60) 2 2 and illustrated in Fig. 9. By deﬁnition there are no deterministic dynamics along this line and so one might expect that transitions between the absorbing states occur with equal probability. However employing the fast-variable elimination described in the introduction to this section, we ﬁnd that in ﬁnite populations, there is noise-induced (equivalently, drift-induced) selection along the line: A(z) = 0 , 3 2 1 z(1 − z)(1 + z) (1 + z ){1 + (2 − z)z[4 + z(6 + z)]} A (z) = , (61) 4N [1 + z(2 + 3z)] 1 z(1 + z) {1 + z [1 + z(6 − z(6 + (3 − z)z))]} B(z) = , [1 + z(2 + 3z)] (4) S where z is the value of 1 − 2x on theCM[seeEq. (60)]. We ﬁnd that A (z) is positive along the entire length of the CM (see Fig. 9). Since z = 0 corresponds to the absorbing state in which the population is entirely comprised of XX and XY genotypes, and z = 1 corresponds to the absorbing state in which the population is entirely comprised of X Yand 123 34 G. W. A. Constable, A. J. McKane Fig. 9 Left: Figure illustrating the neutral dynamics of the model of transitions from male to female heteroga- mety. Only the frequency of female genotypes is shown. The blue dashed line is the CM, deﬁned by Eq. (60), with male genotypes assumed here to be ﬁxed at their value on the CM. Grey arrows indicate deterministic trajectories starting at various points and eventually reaching the CM. Right: Form of noise-induced selection on theSS[seeEq. (61)] (Color ﬁgure online) YY genotypes, we can conclude that the dominant mutation X is selected for along the entire CM. Therefore, if a dominant X mutation occurs in a resident XX|XY population, it is more likely to invade and ﬁxate than a recessive X mutation in a resident X Y|YY population. The above picture, whereby dominant sex-determining mutations are more likely to invade and ﬁxate, is recapitulated for successive invasions. If the resident population is X Y|YY, a dominant Y mutation can occur, Y , yielding X Y males. Again ﬁve genotypes can emerge following random mating including X X (female) and Y Y (male). This dominant mutation is once again more likely to invade, creating an emergent directionality to evolution with substitution rates of neutral dominant mutations being ﬁve times higher than those of recessive mutations [61]. Intriguingly, the behaviour described here recapitulates the empirical observation that sex chromosomes typically evolve to include a cascade of inhibitory mutations [64], with more recent mutations at the top of the cascade. The “bottom up” hypothesis [60] therefore pictures the sequential invasion of sex-determining genes, each one of which represses the action of the previous mutation. While it can be argued that the drift-induced phenomenon described here is not solely responsible for this pattern, and other deterministic explanations have been suggested [59], it is nevertheless interesting to see the emergence of the pattern in a model with minimal assumptions. In [61], the above analysis is extended to include biologically relevant selection. In addition more complicated models of transitions between sex chromosomes are also explored via the same timescale-separation arguments. 5 Discussion Our aim in this paper has been to describe and review an approach to the systematic mod- elling and analysis of a class of models in population genetics. Readers with a background in theoretical physics, and statistical physics in particular, will have found many of the ideas and techniques familiar. These include: the care taken in distinguishing between the form of the models at the microscale, the mesoscale and the macroscale; the idea that the model 123 Exploiting Fast-Variables to Understand Population Dynamics… 35 should be formulated and unambiguously deﬁned at the microscale and that the correspond- ing mesoscopic and macroscopic models should be derived from the former by some form of systematic approximation procedure; the nature of these approximation procedures (the dif- fusion limit and fast mode elimination); the formalism of non-equilibrium statistical physics including master equations, FPEs and SDEs. While the techniques and philosophy that we use come from statistical physics, the moti- vation and the models are taken from population biology, and especially from population genetics. Due to the difﬁculty of performing analytical calculations with models that have many variables and parameters, researchers in these ﬁelds frequently create models for rather speciﬁc situations and to answer one very limited and deﬁnite question. These models may be restricted to one particular situation, but they at least have the possibility of being anal- ysed. On the other hand the danger with this way of doing things is that the subject becomes characterised by a series of models with conﬂicting assumptions and little overlap. We have taken a different path: we have not held back from setting-up general models, even though they typically contain very many variables and parameters, but have instead tried to system- atically approximate these models to obtain ones which are more tractable. In this way the variables that we focus on are those which are likely to be the most important in the medium to long term (the slow modes) and the parameters which are chosen as being important at late times are combinations of the initial parameters of the original model which would have been impossible to guess apriori. There are also some modellers who would argue that to study more complex models, and especially those which feature individuals with more than one or two attributes and which have many parameters, one should turn to agent based models. Here there is no attempt at mathematical analysis and a computer simulation is set up in which the agents (i.e. the individuals) are simply allowed to interact with each other according to a set of rules. These rules will have many parameters associated with them, which will not in general easily map onto the parameters of the type we deﬁned here when formulating the microscopic models. While agent based models might be useful in some situations, the difﬁculty with the whole approach is that one may have no idea what aspect of the model is responsible for a particular outcome that one is interested in. In fact, this is frequently put forward as a strength of the agent-based method: behaviour emerges without having to have been programed-in. We believe that our approach straddles those where the models are simple and tailored for a particular situation and those agent based models which allow for complex systems, but which may be unable to provide insight into the underlying reasons for particular outcomes. Of course, fast-variable elimination has a long history in theoretical biology and ecology. Most prominently, the form of the Michaelis-Menten function [2] and the Hollings type II functional response [16] are derived based on arguments that assume a separation of timescales. However, these techniques tend to be employed most frequently in deterministic settings. The ﬁeld of population genetics provides a notable exception to this rule, with the crucial role of genetic drift providing motivation for studies that employed fast-variable elimination in a stochastic setting. We have explicitly discussed the classical assumption of diploid populations at Hardy–Weinberg frequencies in Sect. 3.3,whereaseparationof timescales was assumed apriori. However, there also exists a range of studies which apply more sophisticated techniques (see, for example [18,19,32,45,46,58]). The relatively frequent occurrence of stochastic fast-variable elimination in population genetics is the result of two factors. The ﬁrst, as we have already mentioned, is that population genetics models are often inherently stochastic. The second is that a common assumption in many population genetic models (even when fast-variable elimination techniques are not 123 36 G. W. A. Constable, A. J. McKane required) is that selection is weak relative to the other processes. If then the question we wish to ask is which of two alleles (or genotypes) ﬁxes in a population, the stage is set for methods based on the elimination of fast-variables to be useful: by deﬁnition, if there is no selection, there must be a CM connecting states in which either allele is ﬁxed; if selective forces are weak, the CM is simply replaced by an SS; if there are only two alleles competing, regardless of the other variables, the SS is one-dimensional, allowing an analytic treatment of the resulting effective equation. Given the ubiquity of the above listed features in models of population genetics, it is there- fore perhaps surprising that there are in fact not many more studies employing fast-variable elimination. Perhaps some of the reasons for this stem from the apparent difﬁculty of dealing with these types of problem, with previous studies [58] relying on the perturbative techniques described by Gardiner [24]. While certainly thorough, this approach is perhaps lacking in intuition, taking place as it does within the context of the FPE and involving operators. In contrast, we believe that there is much to say for the approaches that we have outlined in this paper. As they are described within the context of SDEs (entirely equivalent to the FPE), and as the effective equation is constructed from physically meaningful quantities (such as the left and right eigenvectors), it is relatively straightforward to apply the methodology. We note that while other approaches to fast-variable elimination in the SDE setting are well practiced (see [11], which gives a more complete review of this area of the literature), these too have not cemented themselves as methods within the population genetics literature. The fact that, as described in Sect. 4.2, the techniques we have presented are yielding novel insights to models ﬁrst developed in the 1970s is testament to the untapped potential of the approach. In terms of applications of the method, we have tried to some extent to move away from using the Moran model, which has been a very popular starting point among some researchers over the last decade or two. Instead we have introduced competition between individuals to control population growth. There are many other types of interaction between individuals that could be included, and this is an interesting question for the future. But even when only competition between individuals is included, some interesting points emerge. For example, as discussed in Sect. 3.2.2, the model with competition can be mapped into a Moran model, but only if the selection is frequency dependent, which suggests that the traditional use of frequency independent selection may be too restrictive. In addition, the mapping gives a direct connection between the competition rates and the elements of the payoff matrix. Only relative payoffs of a certain type appear in the mapping between the two models, which again could be used to constrain the nature of the payoffs. This would be useful, since the lack of prescriptions available to guide the choice of payoffs is a signiﬁcant weakness in the application of game theory to the subject. Of course, our methods are only valid when selection is weak, and moving beyond this regime, such as when strong mutualisms are accounted for, is known to invalidate this mapping between the SLVC and models of game theory [4]. Extending such a study using the methods we have presented is also an interesting avenue for future studies. Perhaps one of the most surprising insights that fast-variable elimination has provided is the occurrence of noise-induced dynamics in otherwise deterministically neutral systems. These dynamics, described in Sect. 4, can be interpreted as drift-induced selection and can give rise to surprising behaviours, such as selection reversal. While historically such behaviour has been observed and characterised in simpler settings [27], it has been often dismissed as a small second order effect of little biological signiﬁcance [29]. This is because noise induced selection is inversely proportional to the size of the system, and therefore, it is argued, likely to be swamped by other processes. However, this misses the fact that when selection is weak 123 Exploiting Fast-Variables to Understand Population Dynamics… 37 (a frequent assumption in population genetics) or when a population has a small effective population size (such as when a system is spatially distributed) noise-induced terms can in fact dominate the dynamical behaviour. Moreover, utilisation of fast-variable elimination techniques has led to an increasing awareness that noise-induced selection appears in many distinct evolutionary systems in ecology, epidemiology and population genetics. We have already discussed some of these in detail in Sect. 4. However more examples are arising with increasing frequency. In [39], stochastic Lotka–Volterra models for multiple islands (similar to the models described in Sect. 3.2.1) were used to show that demographic stochasticity induces a selective bias for species that disperse at a higher rate between identical islands, where there is no selection deterministically. Intriguingly, this bias persists even when the islands are inhomogeneous and a deterministic selective pressure exists for slower dispersal [40]. In [37] meanwhile, a model of viral evolution was investigated, showing that noise-induced dynamics selected for viral strains with a fast infection and recovery time in SIS and SIR models. It is clear that the methods we have described in this paper hold a great degree of utility, whether one is seeking to consolidate existing disparate models, add more biological detail to well understood systems or uncover and characterise entirely new dynamics. Beyond drawing attention to this fact, we hope that as presented, this work will prove of value to readers with both biology and statistical physics backgrounds. For the biologists, this paper has aimed to show that these techniques are not as conceptually formidable as other studies may have suggested. For the statistical physicists, this paper has aimed to show that many interesting and relevant problems in biology still exist, as yet untapped, and amenable to analysis with their approaches. Acknowledgements GWAC thanks the Finnish Center for Excellence in Biological Interactions and the Leverhulme Early Career Fellowship provided by the Leverhulme Trust for funding. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 Interna- tional License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix A: Calculational Details and Background for Sect. 3.3 In this appendix we give analytical support for the system discussed in Sect. 3.3.Webegin by brieﬂy reviewing the classic result (which assumes a separation of timescales apriori) before moving on to the fast-variable elimination itself. Classic Approach The classic treatment begins by considering the dynamics of alleles, rather than the genotypes that carry them. The probability of an allele reproducing is propor- tional to its frequency in the population (which is assumed to be given by Hardy–Weinberg frequencies) weighted by the ﬁtness of the genotypes carrying the allele. So, if there is a (1) (A ) (1) frequency n /(2N ) of allele A , then at Hardy–Weinberg frequencies the frequency of (1) (1) (1) (A ) 2 (1) (2) A A genotypes is approximately [n /(2N )] and the frequency of A A geno- (1) (1) (A ) (A ) 2 (1) typesis2n (N − n )/(2N ) [see Eq. (37)]. Assuming that the A allele has a ﬁtness (1) (1) (1) (2) 1 + s when it occurs in an A A genotype, and a ﬁtness 1 +sh when it occurs in a A A 123 38 G. W. A. Constable, A. J. McKane genotype, the transitions rates can be written (1) Weighted prob. A birth & '( ) ⎡ ⎤ (1) (1) (1) (A ) (A ) (A ) (1) (1) n n n (A ) (A ) ⎣ ⎦ T (n + 1|n ) = (1 + s) + (1 + sh) 1 − 2N 2N 2N (2) neutr al pr ob. A death & '( ) (1) (A ) × 1 − , 2N ⎡ ⎤ (1) (1) (1) (A ) (A ) (A ) (1) (1) n n n (A ) (A ) ⎣ ⎦ T (n − 1|n ) = (1 + sh) 1 − + 1 − 2N 2N 2N ( )& ' (2) Weighted prob. A birth (1) (A ) × . 2N ( )& ' (1) neutr al pr ob. A death (1) (2) (1) Note that since only “half” of the A A genotype carries the A allele, an extra factor (1) (2) (1) (A ) (A ) (A ) of a half arises in terms that involve both n and n = 2N − n . Applying the (1) (1) (A ) diffusion approximation with x = n /2N and assuming s is small one obtains Eq. (38). Note that since the diffusion approximation is here conducted with a large parameter 2N,an additional factor 1/2 occurs before the diffusion term in the FPE formulation [see Eq. (6)] and the noise term in the SDE formulation [see Eq. (7)]. Mechanistic Approach We want to build a similar model, but more mechanistically using (1) (1) (1) (2) (1) (2) the genotype frequencies. We will denote their numbers n (A A ), n (A A )and (3) (2) (2) (1) (2) (3) n (A A ), so that the ﬁxed size of the population is n + n + n = N. (Note (1) (1) (1) (2) (A A ) (1) (A A ) (2) the slight difference in notation with the main text; n ≡ n , n ≡ n and (2) (2) (A A ) (3) (αβ) n ≡ n ). We now denote by W the ﬁtness of genotype α reproducing with (αβ) genotype β, and assume W is symmetric. The transition rates may now be constructed by accounting for the probability of each different pairing and death. For instance, the probability (1) that n increases by 1 only is the probability that any pairing occurs (given by the product of frequencies of the genotypes in the pairing) multiplied by the probability that pairing gives (1) (1) rise to an A A genotype by Mendelian inheritance, multiplied by the probability that an (2) (2) (2) (2) A A genotype dies (given by the frequency of A A in the population). Accounting for these terms for each possible transition gives us the following transition rates: 1 1 (1) (2) (1) (2) (11) (1) 2 (12) (1) (2) T (n + 1, n |n , n ) = W (n ) + 2 W n n N 2 (22) (2) 2 (1) (2) + W (n ) (N − n − n ), 1 1 1 (1) (2) (1) (2) (22) (2) 2 (23) (2) (1) (2) T (n − 1, n |n , n ) = W (n ) + 2 W n (N − n − n ) N 4 2 (33) (1) (2) 2 (1) + W (N − n − n ) n , 123 Exploiting Fast-Variables to Understand Population Dynamics… 39 1 1 (1) (2) (1) (2) (12) (1) (2) (13) (1) (1) (2) T (n , n + 1|n , n ) = 2 W n n + 2W n (N − n − n ) 1 1 (22) (2) 2 (23) (2) (1) (2) + W (n ) + 2 W n (N − n − n ) 2 2 (1) (2) N − n − n , 1 1 1 (1) (2) (1) (2) (22) (2) 2 (23) (2) (1) (2) T (n , n − 1|n , n ) = W (n ) + 2 W n (N − n − n ) 4 2 (33) (1) (2) 2 (2) + W (N − n − n ) n , 1 1 (1) (2) (1) (2) (11) (1) 2 (12) (1) (2) T (n + 1, n − 1|n , n ) = W (n ) + 2 W n n (22) (2) 2 (2) + W (n ) n , 1 1 (1) (2) (1) (2) (12) (1) (2) (13) (1) (1) (2) T (n − 1, n + 1|n , n ) = 2 W n n + 2W n (N − n − n ) 1 1 (22) (2) 2 (23) (2) (1) (2) (1) + W (n ) + 2 W n (N − n − n ) n . 2 2 (αβ) (αβ) We now set W = 1 + sα and assume s is small to obtain a diffusion approximation (1) (1) (2) (2) [see Eqs. (6)and (7)] with y = n /N and y = n /N; (2) (1) (1) (1) (2) (1) (33) (1) A ( y) =−y (1 − y ) + y (y + ) + s −α y (11) (13) (33) (1) 2 (12) (23) (33) (1) (2) + (α − 2α + 2α )(y ) + (α − 2α + 2α )y y (22) (2) 2 (1) 3 (11) (13) (33) (1) 2 (2) + α (y ) − (y ) (α − 2α + α ) − (y ) y (12) (13) (23) (33) (1) (2) 2 (22) (23) (33) × (2α − 2α − 2α + 2α ) − y (y ) (α − 2α + α ) , (2) (2) (1) (1) (2) (1) (13) (1) A ( y) = 2 y (1 − y ) − y (y + ) + s 2α y (2) (23) (33) (13) (1) 2 (1) (2) (12) (13) (23) (33) + y (α − α ) − 2α (y ) + y y (α − 4α − α + 2α ) (2) 2 (22) (23) (33) (1) 2 (2) (11) (13) (33) + (y ) (α − 6α + 4α ) − (y ) y (α − 2α + α ) (1) (2) 2 (12) (13) (23) (33) (2) 3 (22) (23) (33) − 2y (y ) (α − α − α + α ) − (y ) (α − 2α + α ) , (62) and (11) (1) (1) (2) 2 (1) 3 (1) 2 (2) B ( y) = y (1 + y ) + (y ) − 2(y ) − 2(y ) y (1) (2) (2) +y y 1 − y , (22) (1) (2) (1) 2 (1) (2) (2) 2 (2) 3 B ( y) = 2(y + y ) − 2(y ) − 6y y − (y ) + (y ) (1) 2 (2) (1) (2) 2 +4(y ) y + 4y (y ) , 1 1 (12) (1) 2 (1) (2) (2) 3 (1) 3 (1) 2 (2) (1) (2) 2 B ( y) =−2(y ) − y y − (y ) + 2(y ) + (y ) y − y (y ) , 4 2 123 40 G. W. A. Constable, A. J. McKane (21) (12) B ( y) = B ( y). (63) We note however that unlike in the classic approach (where the noise correlation matrix −1 −1 has a prefactor N /2) this noise term has the standard pre-factor of N . Our end goal is to compare the mechanistic description in terms of genotypes given by Eqs. (62)and (63) to that proposed by Eq. (38). To this end the following calculation will involve two steps. The ﬁrst will be removing the fast degrees of freedom in the model (see Fig. 6) as described in the main text. The second will be transforming from the genotype (1) (2) (1) frequencies of the mechanistic model (y and y ) into the allele frequencies (x )in which the classic approach is couched. We note that in both these steps there are some additional mathematical complications that arise that, while formally important, are not crucial to understanding the basic principles of the method or insights developed. In the ﬁrst step, removing the fast-variables, a noise-induced selection term arises. Such terms are discussed in Sect. 4. In the second step, the non-linear transformation we employ makes use of the Ito¯ calculus necessary, giving rise to extra terms in the transformed equa- tion [53]. However in the problem at hand we will see that these two effects cancel precisely. Readers who only wish to understand the principles of the problem and physical intuitions developed are therefore free in this section to bypass references to noise-induced dynamics and Ito¯ calculus. In order to employ the fast-variable elimination procedure, we begin with Eq. (62). From this we can determine that a CM exists, given by Eq. (37). We wish to remove the fast variables as described in the Sects. 2 and 4. To begin, we must identify the left and right eigenvectors {1} {1} of the Jacobian evaluated on the CM that correspond to the zero eigenvalue, u and v . These are given by √ 1 2 z {1} {1} u = z , v = √ , (64) 1 1 − 2 z where z is y on the CM. We can now use these, along with Eqs. (16)–(17)and Eq.(52)in the main text to calculate the effective SDE. We ﬁnd (23) (33) (12) (13) (22) (23) (33) A(z) = 2sz(1 − z)(α − α + 3z(α − α − 2α + 3α − α ) (13) (22) (23) (33) + z(α + 2α − 6α + 3α ) 3/2 (11) (12) (13) (22) (23) (33) +z (α − 4α + 2α + 4α − 4α + α )) , (65) 3/2 2 B(z) = 4z − 4z , (66) from our standard analysis (see Sect. 2). From our more involved analysis given in Sect. 4 we also ﬁnd a noise-induced term A (z) = z − z . (67) We withhold attaching any biological meaning to this for the time being, however we can note that it arises entirely as a result of the curvature of the CM since the trajectories to the CM are not divergent (see Sect. 4). (1) (1) We now want to make a change of variables from z into x = z,where x is the (1) (1) (1) 2 frequency of A alleles in the population (recall from Eq. (37)that y = (x ) on the (1) CM and that z is our variable y on the CM). This is a non-linear transformation, and we must therefore employ Ito¯ calculus [53]. This yields a reduced SDE with , - , - (1) 2 (1) dx 1 d x (1) (1) 2 S (1) 2 (1) 2 ˜ ¯ ¯ ¯ A(x ) = A[(x ) ]+ A [(x ) ] + B[(x ) ] , (68) dz 2N dz 123 Exploiting Fast-Variables to Understand Population Dynamics… 41 (1) dx (1) (1) 2 ˜ ¯ B(x ) = B[(x ) ] . (69) dz (1) (1) Here the term in A(x ) involving a second derivative of x is that which arises from a proper implementation of Ito¯ calculus. However, after some algebra we ﬁnd that the term S (1) 2 (1) ¯ ˜ A [(x ) ] cancels this Ito¯ term, and so we ﬁnd the above expression for A(x ) simpliﬁes to (1) dx (1) (1) 2 ˜ ¯ A(x ) = A[(x ) ] , (70) dz (1) which is written in full in Eq. (40), while B(x ) is givenbyEq. (41). Note that in this particular case, these expressions are those that we would have obtained if we had naively ignored both noise-induced selection and Ito¯ calculus. Finally, we draw attention to the fact that once we account for the extra factor 1/2 that occurs in front of the noise correlation term (1) B(x ) in the classic approach (which we recall comes from the diffusion approximation −1 resulting in a coefﬁcient of N /2 in the FPE), the noise terms are identical in both the effective and classic equations, Eqs. (41)and (38). (αβ) With the effective SDE now in hand, we are free to choose which values of W = (αβ) 1 + sα we wish. We will look at some very simple parameterisations. First we consider (αβ) a multiplicative ﬁtness function where the ﬁtness W of any pairing of genotypes α and β is product of the ﬁtness of the genotypes; (αβ) (α) (β) W = W W (α) (β) = (1 + sα )(1 + sα ) (α) (β) 2 = 1 + s(α + α ) + O(s ). Substituting this into Eq. (40), we ﬁnd that this is exactly the same form as that posited in (1) (2) (3) classic population genetics [see Eq. (38)] if α = 1, α = h and α = 0. We next consider a system, similarly parameterised, but where the genotype pairings have the average ﬁtness of the two genotypes; (α) (β) W + W (αβ) W = (α) (β) (α + α ) = 1 + s . (1) (2) (3) Again we get a form similar to that in Eq. (38)if α = 1, α = h and α = 0, except now selection acts half as quickly. The case of additive ﬁtness can be tackled in precisely the (αβ) (α) (β) (α) (β) same fashion by setting W = W + W = 2 + s(α + α ), and then rescaling by a factor two. References 1. Blythe, R.A., McKane, A.J.: Stochastic models of evolution in genetics, ecology and linguistics. J. Stat. Mech. 2007, P07018 (2007) 2. Briggs, G.E., Haldane, J.B.S.: A note on the kinetics of enzyme action. Biochem. J. 19, 338–339 (1925) 3. Bull, J.J., Charnov, E.L.: Changes in the heterogametic mechanism of sex determination. Heredity 39, 1–14 (1977) 4. Chotibut, T., Nelson, D.R.: Evolutionary dynamics with ﬂuctuating population sizes and strong mutualism. Phys.Rev.E 92, 022718 (2015) 123 42 G. W. A. Constable, A. J. McKane 5. Chotibut, T., Nelson, D.R.: Population genetics with ﬂuctuating population sizes. J. Stat. Phys. 167, 777–791 (2017) 6. Constable, G.W.A., McKane, A.J.: Fast-mode elimination in stochastic metapopulation models. Phys. Rev. E 89, 032141 (2014) 7. Constable, G.W.A., McKane, A.J.: Population genetics on islands connected by an arbitrary network: an analytic approach. J. Theor. Biol. 358, 149–165 (2014) 8. Constable, G.W.A., McKane, A.J.: Models of genetic drift as limiting forms of the Lotka-Volterra com- petition model. Phys. Rev. Lett. 114, 038101 (2015) 9. Constable, G.W.A., McKane, A.J.: Stationary solutions for metapopulation Moran models with mutation and selection. Phys. Rev. E 91, 032711 (2015) 10. Constable, G.W.A., McKane, A.J.: A mapping of the stochastic Lotka-Volterra model to models of population genetics and game theory. Phys. Rev. E 96, 022416 (2017) 11. Constable, G.W.A., McKane, A.J., Rogers, T.: Stochastic dynamics on slow manifolds. J. Phys. A: Math. Theor. 46, 295002 (2013) 12. Constable, G.W.A., Rogers, T., McKane, A.J., Tarnita, C.E.: Demographic noise can reverse the direction of deterministic selection. Proc. Nat. Acad. Sci. USA 113, E4745–E4754 (2016) 13. Coullet, P.H., Spiegel, E.A.: Amplitude equations for systems with competing instabilities. SIAM J. Appl. Math. 43, 776–821 (1983) 14. Cox, J.C., Ingersoll, J.E., Ross, S.A.: A theory of the term structure of interest rates. Econometrica 53, 385–407 (1985) 15. Crow, J.F., Kimura, M.: An Introduction to Population Genetics Theory. The Blackburn Press, New Jersey (1970) 16. Dawes, J.H.P., Souza, M.O.: A derivation of Holling’s type I, II and III functional responses in predator prey systems. J. Theor. Biol. 327, 11–12 (2013) 17. Dirac, P.A.M.: The Principles of Quantum Mechanics, 4th edn. Clarendon Press, Oxford (1958) 18. Ethier, S.N., Nagylaki, T.: Diffusion approximations of Markov chains with two time scales and applica- tions to population genetics. Adv. Appl. Probab. 12, 14–49 (1980) 19. Ethier, S.N., Nagylaki, T.: Diffusion approximations of Markov chains with two time scales and applica- tions to population genetics, II. Adv. Appl. Probab. 20, 525–545 (1988) 20. Ewens, W.J.: Mathematical Population Genetics, 2nd edn. Springer, Berlin (2004) 21. Fisher, R.A.: On the dominance ratio. Proc. R. Soc. Edinb. 42, 321–341 (1922) 22. Fisher, R.A.: The Genetical Theory of Natural Selection. Clarendon Press, Oxford (1930) 23. Fox, R.F., Uhlenbeck, G.E.: Contributions to non-equilibrium thermodynamics. I. Theory of hydrody- namical ﬂuctuations. Phys. Fluids 13, 1893–1902 (1970) 24. Gardiner, C.W.: Handbook of Stochastic Methods. Springer, Berlin (2009) 25. Gillespie, D.T.: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 22, 403–434 (1976) 26. Gillespie, D.T.: Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 81, 2340–2361 (1977) 27. Gillespie, J.H.: Natural selection for within-generation variance in offspring number. Genetics 76, 601– 606 (1974) 28. Hallatschek, O.: Noise driven evolutionary waves. PLoS Comput. Biol. 7, e1002005 (2011) 29. Hansen, T.F.: On the deﬁnition and measurement of ﬁtness in ﬁnite populations. J. Theor. Biol. 419, 36–43 (2017) 30. Hanski, I.: Metapopulation Ecology. Oxford University Press, Oxford (1999) 31. Hardy, G.H.: Mendelian proportions in a mixed population. Science 28, 49–50 (1908) 32. Hossjer, O., Tyvand, P.A.: A monoecious and diploid Moran model of random mating. J. Theor. Biol. 394, 182–196 (2016) 33. Houchmandzadeh, B.: Fluctuation driven ﬁxation of cooperative behavior. Biosystems 127, 60–66 (2014) 34. Houchmandzadeh, B., Vallade, M.: The ﬁxation probability of a beneﬁcial mutation in a geographically structured population. New J. Phys. 13, 073020 (2011) 35. Houchmandzadeh, B., Vallade, M.: Selection for altruism through random drift in variable size popula- tions. BMC Evol. Biol. 12, 61 (2012) 36. Kimura, M., Weiss, G.H.: The stepping stone model of population structure and the decrease in genetic correlation with distance. Genetics 49, 561–576 (1964) 37. Kogan, O., Khasin, M., Meerson, B., Schneider, D., Myers, C.R.: Two-strain competition in quasi-neutral stochastic disease dynamics. Phys. Rev. E 90, 042149 (2014) 38. Lin, Y.T., Kim, H., Doering, C.R.: Features of fast living: on the weak selection for longevity in degenerate birth-death processes. J. Stat. Phys. 148, 646–662 (2012) 123 Exploiting Fast-Variables to Understand Population Dynamics… 43 39. Lin, Y.T., Kim, H., Doering, C.R.: Demographic stochasticity and evolution of dispersion I. Spatially homogeneous environments. J. Math. Biol. 70, 647–678 (2015) 40. Lin, Y.T., Kim, H., Doering, C.R.: Demographic stochasticity and evolution of dispersion II. Spatially inhomogeneous environments. J. Math. Biol. 70, 679–707 (2015) 41. Maddamsetti, R., Lenski, R.E., Barrick, J.E.: Adaptation, clonal interference, and frequency-dependent interactions in a long-term evolution experiment with escherichia coli. Genetics 200, 619–631 (2015) 42. Maruyama, T.: On the ﬁxation probability of mutant genes in a subdivided population. Genet. Res. Camb. 15, 221–225 (1969) 43. McKane, A.J., Biancalani, T., Rogers, T.: Stochastic pattern formation and spontaneous polarisation: the linear noise approximation and beyond. Bull. Math. Biol. 76, 895–921 (2014) 44. Moran, P.A.P.: Random processes in genetics. Math. Proc. Camb. Phil. Soc. 54, 60–71 (1957) 45. Nagylaki, T.: The strong migration limit in geographically structured populations. J. Math. Biol. 9, 101– 114 (1980) 46. Newberry, M.G., McCandlish, D.M., Plotkin, J.B.: Assortative mating can impede or facilitate ﬁxation of underdominant alleles. Theor. Popul. Biol. 112, 14–21 (2016) 47. Nowak, M.A.: Evolutionary Dynamics: Exploring the Equations of Life. Harvard University Press, Cam- bridge (2006) 48. Parra-Rojas, C., House, T., McKane, A.J.: Stochastic epidemic dynamics on extremely heterogeneous networks. Phys. Rev. E 94, 062408 (2016) 49. Parra-Rojas, C., McKane, A.J.: Reduction of a metapopulation genetic model to an effective one island model. Pre-print arXiv:1707.07145 (2017) 50. Parsons, T., Quince, C.: Fixation in haploid populations exhibiting density dependence I: the non-neutral case. Theor. Popul. Biol. 72, 121–135 (2007) 51. Parsons, T.L., Rogers, T.: Dimension reduction for stochastic dynamical systems forced onto a manifold by large drift. J. Phys. A: Math. Theor. 50, 415601 (2017) 52. Reichl, L.E.: A Modern Course in Statistical Physics. Wiley VCH, New York (1998) 53. Risken, H.: The Fokker-Planck Equation. Springer, Berlin (1989) 54. Roberts, A.J.: Appropriate initial conditions for asymptotic descriptions of the long term evolution of dynamical systems. J. Aust. Math. Soc. Ser. B 31, 48–75 (1989) 55. Roberts, A.J.: Model Emergent Dynamics in Complex Systems. SIAM, Philadelphia (2015) 56. Roughgarden, J.: Theory of Population Genetics and Evolutionary Ecology: An Introduction. Macmillan, New York (1979) 57. Rousset, F.: Genetic Structure and Selection in Subdivided Populations. Princeton University Press, Oxford (2004) 58. Stephan, W., Charlesworth, B., McVean, G.: The effect of background selection at a single locus on weakly selected, partially linked variants. Genet. Res., Camb. 73, 133–146 (1999) 59. van Doorn, G.S., Kirkpatrick, M.: Transitions between male and female heterogamety caused by sex- antagonistic selection. Genetics 186, 629–645 (2010) 60. van Doorn, G.S.: Patterns and mechanisms of evolutionary transitions between genetic sex-determining systems. Cold Spring Harb. Perspect. Biol. 6, a017681 (2014) 61. Veller, C., Muralidhar, P., Constable, G.W.A., Nowak, M.A.: Drift-induced selection between male and female heterogamety. Genetics 207, 711–727 (2017) 62. Watterson, G.A.: The application of diffusion theory to two population genetic models of Moran. J. Appl. Probab. 1, 233–246 (1964) 63. Weinberg, W.: On the detection of heredity in man (in German). Naturk. Wurttemb. 64, 368–382 (1908) 64. Wilkins, A.S.: Moving up the hierarchy: a hypothesis on the evolution of a genetic sex determination pathway. Bioessays 17, 71–77 (1995) 65. Wright, S.: Evolution in Mendelian populations. Genetics 16, 97–159 (1931)
Journal of Statistical Physics – Springer Journals
Published: Nov 1, 2017
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.