Reduction of a metapopulation genetic model to an effective one-island model^(a)

The subject of population genetics holds a particular fascination for statistical physicists because of the many analogies it has with various models in non-equilibrium statistical mechanics [1,2]. Much of the formalism used by physicists in the study of non-equilibrium systems derives from viewing these as stochastic processes, and is directly applicable to the investigation of models of population genetics [3,4]. The concepts that are frequently of interest there, such as the probability that a particular allele fixes, and the mean time to fixation, are also the focus of attention in many physical systems out of equilibrium [5,6].

As the genetic models have become increasingly complicated, incorporating spatial structure, sexual reproduction or several gene loci, the methods of solution previously employed are no longer efficacious. The purpose of this article is to describe a systematic method for reducing the full models to effective models, which still provide good predictions for quantities relating to the fixation of alleles, but which are simple enough to be analysed mathematically. This method has previously been applied to several models of population genetics; here we apply it to a model not previously considered. In this way we hope that the article has the dual function of serving as a concise review of the approach, but also providing some original results.

The specific model we will discuss will have a spatial aspect: several subpopulations in distinct regions, with individuals able to migrate to one region from another. It will therefore have many parameters: birth, death and competition rates which differ between alleles and between regions; the regions, in turn, vary in size (in the sense that they can sustain different numbers of individuals), and the migration rates between them are also variable. We are therefore confronted with the difficulty of analysing a rather complex model, as discussed above. This is managed by making two approximations, which we will show give an excellent agreement with results found by simulating the original model.

The first is the standard diffusion approximation [7], which in the language of statistical physics consists of moving from the microscopic description in terms of individuals to a mesoscopic description in terms the fraction of the population in the various regions that is of one type or the other. The second approximation is the neglect of degrees of freedom that decay rapidly on time scales that are of interest to us. This approximation also has a long history, and is known variously as adiabatic elimination [8], fast variable elimination [9], centre manifold (CM) theory [10], among others. In the present application it will turn out that all degrees of freedom but one, decay away relatively quickly, leaving an effective theory which is sufficiently simple to be analytically tractable.

The modelling procedure that we will adopt will include the effects of migration, selection and genetic drift, but the processes of birth and death will be taken to be distinct, unlike the conventional approach in population genetics where birth and death are coupled in order to keep the population size fixed [11–13]. Instead, a competition between the individuals in the system will be introduced that will have the effect of keeping the population fixed on average, but with ever present fluctuations about this average. In this way the basic elements of the model will more closely resemble an ecological model with the processes of birth, death and competition, but where the different species are identified by the fact that they carry different alleles. We will only examine the case of a single gene in haploid individuals that can only have two variants; we will refer to the alleles as type 1 and type 2. The method can be extended to diploid and multiallelic individuals, but here we prefer to focus on the effects of spatial structure, selection due to varying birth, death and competition rates between the species, and genetic drift due to stochastic effects resulting from the finite number of individuals present in the system.

We seek to make the model as generic as possible, and so we will construct it at the fundamental level of individuals undergoing the processes of birth, death, competition and migration. The simplest choices for these processes lead to a Lotka-Volterra competition model [14], and since the model will be stochastic, we will refer to it as a stochastic Lotka-Volterra competition (SLVC) model. The spatial structure will be introduced by asking that the population is divided into ${\cal D}$ subpopulations in distinct regions. In population genetics these might be referred to as demes or islands; here we will use the terminology of islands, following the practice in ecology. Similarly we will refer to the population as a metapopulation [15], since it will have the structure of a network where the nodes are islands, with different sizes and with varying link strength (level of migration) between them.

As we have stressed above, we believe it is important to begin at the level of discrete individuals and the interactions between them. As also mentioned, in common with most authors, we make the diffusion approximation [7] in order to make progress in analysing the model. Within this approximation the variables are the number density of individuals of type α on island i, denoted by $x^{(\alpha)}_i$. The parameters of the model are both local (the birth and death rates of these individuals, respectively $b^{(\alpha)}_i$ and $d^{(\alpha)}_i$, and the competition between types α and β on island i denoted by $c^{(\alpha \beta)}_i$) and non-local (the rate $\mu_{ij}$ at which an individual from island j will migrate to island i). The specification of the model and the application of the diffision approximation is by now standard [7,16], and is discussed in detail for this particular model in sect. 2 of the supplementary material Supplementarymaterial.pdf (SM). Our interest here is in the second approximation discussed in the introduction, which can be made after this first approximation has been carried out, and therefore our starting point will be the stochastic differential equation which is the outcome of the analysis described in sect. 2 of the SM.

To simplify the form of the stochastic differential equation it is useful to introduce an index I that runs from 1 to $2{\cal D}$, so that I = i if the allele labelled is 1 and if the island being considered is i, and $I={\cal D}+i$ if the allele labelled is 2 and if the island being considered is i. The state of the system is denoted by the vector $\bm{x} = (x^{(1)}_1,x^{(2)}_1,\ldots,x^{({\cal D})}_1,x^{(1)}_2,x^{(2)}_2,\ldots,x^{({\cal D})}_2)$. As discussed in sect. 2 of the SM, the model also contains a set of ${\cal D}$ parameters, V_i, which denote the potential capacity of island i, both in terms of environmental factors required to sustain a population and the size of the island. Within the diffusion approximation we set $V_i = \beta_i V$, where $\beta_i$ is a number of order one that characterises the capacity of each island compared to the others, and where V is the typical carrying capacity of an island, which is the central parameter which controls the diffusion approximation. After these definitions, we may now write the stochastic differential equation (defined in the sense of Itō [5]) derived in the SM in the form

$$\begin{equation} \frac{\mathrm{d}x_I}{\mathrm{d}\tau} = A_I(\bm{x}) + \frac{1}{\sqrt{V}} \eta_I(\tau), \end{equation} \tag{ 1 }$$

where $\tau = t/V$ is a rescaled time and $\eta_I(\tau)$ is a Gaussian white noise with zero mean and with a correlator

$$\begin{equation} \langle \eta_I(\tau) \eta_J(\tau')\rangle = B_{I J}(\bm{x}) \delta( \tau - \tau'). \end{equation} \tag{ 2 }$$

The functions $A_I(\bm{x})$ and $B_{I J}(\bm{x})$ which specify the model, are derived in sect. 2 of the SM, beginning from the microscopic description given by eqs. (SM1)–(SM5). They are given by

$$\begin{eqnarray} \hspace{-1pc}A^{(\alpha)}_i(\bm{x}) &= \frac{1}{\beta_i} \Biggl\{ \left( b^{(\alpha)}_i - d^{(\alpha)}_i\right)x^{(\alpha)}_i - \sum^{2}_{\beta = 1} c^{(\alpha \beta)}_i x^{(\alpha)}_i x^{(\beta)}_i \nonumber \\ \hspace{-1pc}&\quad +\,{\cal M}^{(\alpha,-)}_i \Biggr\}, \ \ \ \ i=1,\ldots,{\cal D}, \ \ \ \alpha=1,2, \end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray} \hspace{-1pc}B^{(\alpha \alpha)}_{ii}(\bm{x}) &= \frac{1}{\beta^2_i} \Biggl\{ \left( b^{(\alpha)}_i + d^{(\alpha)}_i\right)x^{(\alpha)}_i + \sum^{2}_{\beta = 1} c^{(\alpha \beta)}_i x^{(\alpha)}_i x^{(\beta)}_i \nonumber \\ \hspace{-1pc}&\quad +\,{\cal M}^{(\alpha,+)}_i \Biggr\}, \ \ \ \ i=1,\ldots,{\cal D}, \ \ \ \alpha=1,2, \end{eqnarray} \tag{ 4 }$$

where the non-local contributions due to migration, ${\cal M}^{(\alpha,\pm)}_i$, are given by

$$\begin{equation} {\cal M}^{(\alpha,\pm)}_i = \sum_{j \neq i}\left[ \mu_{ij} x^{(\alpha)}_j \pm \mu_{ji} x^{(\alpha)}_i \right]. \end{equation} \tag{ 5 }$$

In addition,

$$\begin{equation} B^{(\alpha \alpha)}_{ij}(\bm{x}) = - \frac{1}{\beta_i \beta_j}\,\left[ \mu_{ij} x^{(\alpha)}_j + \mu_{ji} x^{(\alpha)}_i\right], \ \ (i \neq j), \end{equation} \tag{ 6 }$$

and $B^{(12)}_{ij}=B^{(21)}_{ij}=0$, for all i, j.

While the transition rates which define the model at the level of individuals (given by eq. (SM1)) are rather transparent, and can be written down from the model description, the forms of the equivalent mesoscopic quantities $A_I(\bm{x})$ and $B_{IJ}(\bm{x})$ given above are rather less clear. The $A_I(\bm{x})$, from which the deterministic dynamics follow, has some familiar elements, namely the first two terms in the curly brackets which are the usual Lotka-Volterra interaction terms. So although analytic progress is helped by making the diffusion approximation, the fact that the functions given by eqs. (3) and (4) are still very complex, means that further approximations are required. We will now show that the elimination of fast modes is an approximation which can be justified biologically, and yields a simplified model which retains the power to make accurate predictions for quantities such as probabilities of fixation and mean fixation times.

In this second approximation, the mesoscopic model with $2{\cal D}$ degrees of freedom may be reduced to one with only a single degree of freedom. This reduced model can essentially be thought of as one with no spatial structure, but defined by a set of effective parameters, which encapsulate those of the full model. Later we will compare the result of calculations from the reduced model to numerical simulations of the original.

The method is based on the observation that the dynamics of the full model consists of two stages. The first consists of a relatively rapid decay from the initial state to the vicinity of a CM (if selection is absent) or a slow subspace (SS) (if selection is present). It then enters the second stage where it wanders stochastically on or near the CM (and also weakly deterministically on a SS if weak selection is present) until fixation of one or other of the alleles; this is shown in fig. 1 for a neutral system with ${\cal D}=5$ islands. This is the heart of the time-scale separation: the rate of migration, which controls the collapse onto the SS, is much greater than the rate of genetic drift, which eventually leads to global fixation. Time-scale separation arguments have also been used on similar models elsewhere [17,18]. In the dynamics of the first stage, stochastic effects play very little role; there is what is in essence a deterministic collapse onto the CM (or SS). We will therefore study this first stage of the processs deterministically, beginning with the case of no selection, where a true CM exists.

Zoom In Zoom Out Reset image size

Neutral model

In SLVC models, selection is introduced through the parameters $b^{(\alpha)}_i, d^{(\alpha)}_i$, and $c^{(\alpha \beta)}_i$, which if made to vary with α and β, give a selective advantage to those individuals carrying either allele α or allele β. Therefore to have no selection we set $b^{(\alpha)}_i = b^{(0)}_i$, $d^{(\alpha)}_i = d^{(0)}_i$, and $c^{(\alpha \beta)}_i = c^{(0)}_i$ for all $\alpha, \beta=1,2$. Substituting this into the deterministic equation $\mathrm{d}x_I/\mathrm{d}\tau = A_I(\bm{x})$, obtained by taking the $V \to \infty$ limit of eq. (1), yields

$$\begin{eqnarray} \frac{\mathrm{d}x^{(\alpha)}_i}{\mathrm{d}\tau} &= \frac{1}{\beta_i} \Biggl\{\left( b^{(0)}_i - d^{(0)}_i\right)x^{(\alpha)}_i - c^{(0)}_i x^{(\alpha)}_i \sum^{2}_{\beta = 1} x^{(\beta)}_i \nonumber\\ &\quad +\,{\cal M}^{(\alpha,-)}_i \Biggr\}, \ \ \ \ i=1,\ldots,{\cal D}, \ \ \ \alpha=1,2. \end{eqnarray} \tag{ 7 }$$

To achieve the maximum reduction, we are searching for a low-dimensional CM. In this case we can find one which is one dimensional, by seeking fixed points of eq. (7) that are independent of i, that is, solutions of

$$\begin{equation} x^{(\alpha)} \left[\left( b^{(0)}_i + q_i - d^{(0)}_i\right) - c^{(0)}_i \sum^{2}_{\beta = 1} x^{(\beta)} \right] = 0, \ \ \ \alpha=1,2, \end{equation} \tag{ 8 }$$

where

$$\begin{equation} q_i \equiv \sum_{j \neq i}[\mu_{ij} - \mu_{ji}]. \end{equation} \tag{ 9 }$$

The only solution of eq. (8), apart from the trivial solution $x^{(1)}=x^{(2)}=0$, is

$$\begin{equation} x^{(1)} + x^{(2)} = \frac{(b^{(0)}_i + q_i - d^{(0)}_i)}{c^{(0)}_i}, \end{equation} \tag{ 10 }$$

which, for consistency, requires that $( b^{(0)}_i + q_i - d^{(0)}_i) = \kappa c^{(0)}_i$ for all i, where κ is a constant. This condition should perhaps not be surprising, since we are reducing the model from one with $2{\cal D}$ degrees of freedom to one with only one degree of freedom (x⁽¹⁾, with x⁽²⁾ determined from eq. (10)). Therefore each island has in some sense to be neutral in order to obtain a neutral one-island model. Later, when we introduce selection, we will be able to move away from this assumption.

Equation (10) defines the one-dimensional CM, which we show for a two-island system in the phase diagram of fig. 2 —fig. SM1 (see SM) further shows that the solution is island independent. Before proceeding any further, we scale the original variables of the system, in order to make the analysis more transparent. To do this, we define variables

$$\begin{equation} y^{(\alpha)}_i = \frac{c^{(0)}_i x^{(\alpha)}_i}{(b^{(0)}_i + q_i - d^{(0)}_i)} = \kappa^{-1} x^{(\alpha)}_i, \end{equation} \tag{ 11 }$$

with $i=1,\ldots,{\cal D}$ and $\alpha=1,2$. Then repeating the analysis of this section, but in the $y^{(\alpha)}_i$ variables, rather than in the $x^{(\alpha)}_i$, we find a CM where $y^{(\alpha)}_i = y^{(\alpha)}$ for all i and $\alpha=1,2$, with

$$\begin{equation} y^{(1)} + y^{(2)} = 1. \end{equation} \tag{ 12 }$$

We will choose the CM to be parameterised by y⁽¹⁾ which we will denote by z, the only variable of the reduced system. Then $y^{(2)} = 1 - z$.

Zoom In Zoom Out Reset image size

The more complete analysis carried out in sect. 3 of the SM, involves finding the eigenvalues and eigenvectors of the Jacobian on the CM. The decay rates of the modes associated with the various eigenvectors are proportional to the inverse of the corresponding eigenvalues. In the SM the $2{\cal D} - 1$ "fast" modes are identified; the single slow mode —which is actually static when there is no selection, since it has eigenvalue zero— corresponds to the CM. For the purposes of the general overview presented here, the fast modes simply take the system from its initial condition (IC) to the CM, the initial point of contact being referred to as the initial condition on the CM (CMIC).

As discussed earlier in this section, we assume that in this first part of the dynamics —the decay from the initial condition, $\bm{y}^{\text{IC}}$, to the CM— the deterministic dynamics completely dominates the stochastic dynamics. In effect, this means that it is assumed that the stochastic system still reaches the CM at the point $z^{\text{CMIC}}$ found from the deterministic neutral dynamics, and that this can be used as an initial condition for the second stage of the dynamics, which takes place entirely on the CM. This assumption will be examined in the numerical simulations which are discussed later and in the SM.

Model with selection

To go on to analyse the non-neutral case we write the birth, death and competition parameters as follows:

$$\begin{eqnarray} & b^{(\alpha)}_i = b^{(0)}_i\left(1 + \epsilon \hat{b}^{(\alpha)}_{i}\right); \ \ d^{(\alpha)}_i = d^{(0)}_i \left(1 + \epsilon \hat{d}^{(\alpha)}_{i}\right);\nonumber \\ & c^{(\alpha \beta)}_{i} = c^{(0)}_i \left(1 + \epsilon \hat{c}^{(\alpha \beta)}_i\right). \end{eqnarray} \tag{ 13 }$$

Here is the selection strength. As described in sect. 4B of the SM, we assume that and $V^{-1}$ are of the same order, and therefore keep order terms in $A_{I}(\bm{y})$, but only order one terms in $B_{IJ}(\bm{y})$. The noise correlator will then correspond to the one obtained from the neutral theory (see sect. 4A of the SM).

In order to find $A_{I}(\bm{y})$ to first order in , we write the coordinates on the SS as

$$\begin{equation} y^{(1)}_i=z + \epsilon Y^{(1)}_i + {\cal O}(\epsilon^2), \ \ \ y^{(2)}_i=(1-z) + \epsilon Y^{(2)}_i + {\cal O}(\epsilon^2), \end{equation} \tag{ 14 }$$

where $Y^{(1)}_i$ and $Y^{(2)}_i$ are to be determined. Substituting these coordinates into the expressions for $A^{(1)}_{i}(\bm{y})$ and $A^{(2)}_{i}(\bm{y})$ (see eq. (3)), but restricted to the SS, together with some further analysis, gives eq. (SM29) for the equation of the SS.

So far we have identified the one-dimensional subspace that the system collapses onto (the SS) and have identified the variable which moves the system along this subspace (z). The subspace itself was found by starting from eq. (14) and asking that $A_I(\bm{y})$ only had components along the subspace. We can also ask that the noise only acts along the SS; technically this is best achieved through the construction of a projection operator which in effect projects the stochastic differential equation (1) onto a one-dimensional stochastic differential equation consisting of an effective deterministic function $\bar{A}(z)$, with the noise having an effective correlator $\bar{B}(z)$.

The details of how this projection is carried out are given in the SM where it is shown (see sect. 4) that we arrive at the following form for the stochastic differential equation describing the stochastic dynamics after the fast-mode elimination:

$$\begin{equation} \frac{\mathrm{d}z}{\mathrm{d}\tau} = \bar{A}(z) + \frac{1}{\sqrt{V}} \zeta(\tau), \end{equation} \tag{ 15 }$$

where $\zeta(\tau)$ is a Gaussian noise with zero mean and correlator

$$\begin{equation} \langle \zeta(\tau) \zeta(\tau')\rangle = \bar{B}(z) \delta(\tau - \tau'). \end{equation} \tag{ 16 }$$

Here

$$\begin{equation} \bar{A}(z) = \epsilon z (1 - z)(a_1 + a_2 z) + {\cal O}(\epsilon^2), \end{equation} \tag{ 17 }$$

where

$$\begin{eqnarray} \hspace{-1.3pc}a_1 &= \sum^{\cal D}_{i=1} \frac{u^{\{1\} }_i}{\beta_i}\,\left\{ \left[ \left( b^{(0)}_i \hat{b}^{(1)}_i - d^{(0)}_i \hat{d}^{(1)}_i \right) \right. \right.\nonumber \\ \hspace{-1.3pc}&\quad -\, \left. \left. \left( b^{(0)}_i \hat{b}^{(2)}_i - d^{(0)}_i \hat{d}^{(2)}_i \right) \right] + \kappa c^{(0)}_i \left[ \hat{c}^{(2 2)}_i - \hat{c}^{(1 2)}_i \right] \right\} \end{eqnarray} \tag{ 18 }$$

and

$$\begin{equation} a_2 = - \sum^{\cal D}_{i=1} \frac{\kappa u^{\{1\} }_i c^{(0)}_i \Gamma_i }{\beta_i}, \end{equation} \tag{ 19 }$$

and where we have defined $\Gamma_i \equiv \hat{c}^{(11)}_i - \hat{c}^{(12)}_i - \hat{c}^{(21)}_i + \hat{c}^{(22)}_i$. In addition, $\bm{u}^{\{1\} }$ is the eigenvector of the ${\cal D} \times {\cal D}$ matrix with off-diagonal elements $\mu_{ij}/\beta_i$ and diagonal elements $- \sum_{j \neq i} \mu_{ij}/\beta_i$, having eigenvalue zero.

In the same way, the reduced noise correlator is found to be (see eq. (SM38))

$$\begin{equation} \bar{B}(z) = 2bz(1 - z), \end{equation} \tag{ 20 }$$

where

$$\begin{equation} b = \kappa^{-1} \sum^{\cal D}_{i=1} \frac{[u^{\{ 1\} }_i]^{2}}{\beta^2_i} b^{(0)}_i. \end{equation} \tag{ 21 }$$

We see that the forms for $\bar{A}(z)$ and $\bar{B}(z)$ are similar to those that we might expect from a model with only one degree of freedom, but with the parameters of the model (a₁, a₂ and b) encapsulating some of the structure of the original $2{\cal D}$-degrees-of-freedom model. The reduced stochastic differential equation (15), together with the correlation function in eq. (16) and eqs. (17) and (20), completely describe the stochastic dynamics of the reduced system.

It is straightforward to check that the results obtained above agree with an earlier analysis carried out for a single island, i.e., ${\cal D}=1$ [19]. In the single-island reduction, a further simplication was made, which while not necessary, does simplify the analysis. This consisted in asking that the SS passes through the two points $\bm{y}=(1,0)$ and $\bm{y}=(0,1)$ [19]. The analogue in the present case is the requirement that when z = 1, $y^{(1)}_i=1$ and $y^{(2)}_i=0$, for all i. Similarly when z = 0, $y^{(1)}_i=0$ and $y^{(2)}_i=1$, for all i. If these conditions are not imposed, there is a stochastic drift along the SS until either of the axes is reached and fixation of one of the types is achieved. The imposition of the conditions reduces the number of parameters of the model and ensures that fixation occurs at z = 0 and z = 1. In sect. 4B of the SM we show that these conditions imply that

$$\begin{eqnarray} \left(b^{(0)}_i \hat{b}^{(1)}_i - d^{(0)}_i \hat{d}^{(1)}_i \right) &= \kappa c^{(0)}_i \hat{c}^{(1 1)}_i,\nonumber \\ \left(b^{(0)}_i \hat{b}^{(2)}_i - d^{(0)}_i \hat{d}^{(2)}_i \right) &= \kappa c^{(0)}_i \hat{c}^{(2 2)}_i, \end{eqnarray} \tag{ 22 }$$

where $i=1,\ldots,{\cal D}$. The substitution of the conditions in eq. (22) into eq. (18), leads to a form for eq. (17), at order , which is given by

$$\begin{eqnarray} \bar{A}(z) &= \epsilon z ( 1 - z ) \sum^{\cal D}_{i=1} \frac{\kappa c^{(0)}_i u^{\{1\} }_i}{\beta_i}\,\left[\phi^{(1)}_i - \Gamma_i z \right]\,, \end{eqnarray} \tag{ 23 }$$

where

$$\begin{equation} \phi^{(1)}_i \equiv \hat{c}^{(1 1)}_i - \hat{c}^{(1 2)}_i. \end{equation} \tag{ 24 }$$

This shows that all dependence on the birth and death parameters has been eliminated; the result for $\bar{A}(z)$ only depends on the competition parameters.

In the same way as was done in the general case, effective parameters, which contain information about the full model, can be introduced:

$$\begin{eqnarray} \Gamma_{\text{eff}} &\equiv \sum^{\cal D}_{i=1} \frac{\kappa c^{(0)}_i u^{\{1\} }_i}{\beta_i}\,\Gamma_i, \nonumber \\ \phi^{(1)}_{\text{eff}} &\equiv \sum^{\cal D}_{i=1} \frac{\kappa c^{(0)}_i u^{\{1\}}_i}{\beta_i}\,\phi^{(1)}_i. \end{eqnarray} \tag{ 25 }$$

This then yields

$$\begin{equation} \bar{A}(z) = \epsilon z (1 - z)\,\left[ \phi^{(1)}_{\text{eff}} - \Gamma_{\text{eff}} z + {\cal O}(\epsilon) \right], \end{equation} \tag{ 26 }$$

which has the same form as in the one-island case [19], but now with effective parameters. It should be stressed that the simplification leading to eq. (22) was simply made as a special case which leads to a simpler end result, which can be useful in checking the efficacy of the method; the more general form given by eqs. (17)–(19) should and can be used in general.

Figure 3 shows a phase diagram for a system with ${\cal D}=2$ islands and selection. The rather strong level of selection allows us to clearly appreciate the fact that a CM no longer exists, and the system collapses towards a curved SS instead; on the latter, both deterministic and stochastic dynamics take place. In the next section we will use the reduced model to make predictions, and test these through a numerical simulation of the original model.

Zoom In Zoom Out Reset image size

The purpose of this section is twofold. Firstly, to note that the one-degree-of-freedom model given in the previous section can be analysed mathematically, and to compare the predictions of this reduced model to simulations of the full model. Secondly, to use these results to investigate the quality of the approximations made to obtain the reduced model.

Although the form of the reduced model closely resembles those of one-dimensional stochastic models in population genetics [7], there is one significant difference. This is that $\bar{A}(z)$ is in general cubic in the variable z, rather than having a simple quadratic form such as $sz(1-z)$, where s is a selection coefficient. This difference implies that there is a possibility of an "internal" fixed point —one away from the boundaries at z = 0 and z = 1. One might naively expect that the presence of a stable fixed point would lead to a longer mean time to fixation and an unstable fixed point to a shorter mean time to fixation.

To investigate this, we use the form of $\bar{A}(z)$ given by eq. (26). There is the possibility of an internal fixed point at $z^* = \phi^{(1)}_{\text{eff}}/\Gamma_{\text{eff}}$ if $\Gamma_{\text{eff}} \neq 0$, but clearly we require $0 < z^* < 1$, for this to be an internal fixed point in a biologically relevant regime. If we introduce the quantity

$$\begin{equation} \phi^{(2)}_i \equiv \hat{c}^{(2 2)}_i - \hat{c}^{(2 1)}_i, \end{equation} \tag{ 27 }$$

in an analogous way to $\phi^{(1)}_i$, then we can easily show, as in the one-island case [19], that if $0 < z^* < 1$, then either $\phi^{(\alpha)}_{\text{eff}} > 0$ (for both $\alpha=1$ and $\alpha=2$) or $\phi^{(\alpha)}_{\text{eff}} < 0$ (again for both $\alpha=1$ and $\alpha=2$). We can also investigate the stability of the internal fixed point. A simple calculation shows that the internal fixed point is stable if $\Gamma_{\text{eff}} > 0$ and unstable if $\Gamma_{\text{eff}} < 0$. Since $\Gamma_{\text{eff}} = \phi^{(1)}_{\text{eff}} + \phi^{(2)}_{\text{eff}}$, an internal fixed point exists and is stable if both $\phi^{(\alpha)}_{\text{eff}}$ are positive —as shown in fig. 3— and it exists and is unstable if both $\phi^{(\alpha)}_{\text{eff}}$ are negative.

Two quantities which are of interest to calculate are the fixation probability of a given allele and the mean time to fixation of the system, given a set of initial allele frequencies. These are also useful to test the approximations that have been made to obtain the reduced model, since they are long-time properties in the sense that we expect fixation to occur after the system has reached the SS, and has moved along the SS to reach either z = 0 or z = 1.

To calculate the fixation probability and mean time to fixation, we revert to the formalism of Fokker-Planck equations. The details of the calculation are given in the SM (sect. 5); here we simply compare these results against simulations of the full system, shown in fig. 4 for ${\cal D}=2$ —and ${\cal D}=4$ in fig. SM2 (see SM). When there is no selection, we find that the agreement between theory and simulation is excellent. When selection is present, we also see that in spite of the relatively large values of the selection parameter explored, the calculation carried out to linear order in captures the behaviour of the full system extremely well. Furthermore, we corroborate the supposition that the existence of a stable (respectively, unstable) internal fixed point of the reduced system leads to larger (resp. smaller) values of the fixation time. In fig. 4, we present a version of the system with $\hat{c}^{(11)}_i,\hat{c}^{(22)}_i>0$ and $\phi_i^{(1)},\phi^{(2)}_i<0$ for all i, so that $\phi_{\text{eff}}^{(\alpha)} < 0$, $\alpha=1,2$, yielding an unstable fixed point. This is compared to a version with the signs of $\hat{c}^{(12)}_i$ and $\hat{c}^{(21)}_i$ reversed so that, all the other parameters being equal, in this case $\phi_i^{(1)},\phi^{(2)}_i>0$ for all i and the fixed point is stable. The difference between both scenarios is clearly seen. A stronger effect is observed for the case with ${\cal D}=4$ —see fig. SM2— which shows much longer times to fixation when a stable fixed point is present.

Zoom In Zoom Out Reset image size

In this paper we have investigated a model of metapopulation genetics and shown that, despite its relative complexity, it could be reduced to an effective model with only one degree of freedom. This model is amenable to mathematical analysis.

Our starting point differed from that used by many theoretical population geneticists in so far that we did not use the Wright-Fisher or Moran model in their original microscopic form or in their mesoscopic form obtained through the diffusion limit. Although these models are widely used, they have several disadvantages. We have already mentioned the artifically fixed population size, which is required because the models do not include competition between individuals which potentially leads to a rapid increase in population size. Another example, especially relevant in this paper, is the convoluted way in which the migration process is described in the Moran model.

In the SLVC model, individuals simply migrate at a certain rate, just as they are born, die or compete with each other at a certain rate. Therefore, in eq. (SM1), the transition rates for migration only depend on the population density of the relevant allele on the island from which the migration takes place, j. As a consequence it is linear in this density, but it changes the population size on both island j and on island i where the migrant moves to. By contrast, in the Moran model the transition rates depend on the population density of the relevant allele on both islands. It is quadratic in the densities, although cancellations mean that eventually it turns out to be linear, but still depending on the densities of the relevant allele on both j and i. In addition, the migration process only changes the make-up of the population on island i (by perhaps displacing a resident of that island), but does not change the make-up of the population on island j, since all that happens here is that an offspring of an individual migrates as soon as it is born. The process then, in the SLVC model, is clearly simpler and more intuitive. A disadvantage of the SLVC model is, of course, that it doubles the number of variables, as compared to the Moran model, but it can still be reduced to an effective one-variable model, just as in the case of the Moran model [20,21].

The method we have discussed in this paper can be extended to SLVC models with additional features. For instance, in addition to migration, selection and genetic drift, the process of mutation could be added, as has been done for the Moran model [22]. There are, however, many other effects that could be included: the individuals could be assumed to be diploid, or the effect of more than one loci could be included or other types of ecological interactions could be incorporated. There would then be many types of fast modes, but as long as there was a time-scale separation between these and a few slow modes, there would be the possibility of an effective model with just a few degrees of freedom which would encapsulate the essence of the full model. In this way it may be possible to gain quantitative insights into quite complex models.