Probability of Identity by Descent in Metapopulations

GENERALITIES AND SCOPE

Our purpose is to model and study genetic diversity at individual, colony, and total population levels. The distribution of genetic diversity in the population is affected by five agencies: (1) mutation, (2) reproduction, (3) extinction and subsequent replacement of colonies, (4) migration between colonies and/or from an external source, and (5) population size.

We assume that mutation follows an infinite-allele model (Kimura and Crow 1964), which assumes that every mutation creates an allele that is not currently existing in the population. Under this model, two alleles that are identical in state are also identical by descent. Two types of reproduction systems are considered: asexual reproduction and sexual reproduction. Sexual reproduction includes both random and mixed mating (part of the population mating at random, the other part selfing). An extinct colony can be replaced by a given number of individuals originating from one of the existing colonies (propagule model; Slatkin 1977) or by a given number of individuals randomly sampled over the colonies (migrant pool model; Slatkin 1977). In either case, the initial size of the population is restored instantly. Alternatively, colonies can be gradually recolonized, allowing varying population sizes. Finally, colonies can exchange individuals (seeds) or single gametes (pollen).

We consider two models of basically different character, a fixed-size model where the colonies are equally large and of given size, and a variable-size model in which the sizes of the subpopulations are given by birth-and-death processes. We introduce next the various parameters and the precise dynamics in both models.

Model parameters: In both models there are n colonies. For the fixed-size model we let N denote the number of diploid individuals in each colony and k be the number of colonists. The distribution of genetic material evolves according to the nonnegative Poisson intensities

• u, mutation of an allele into a new type per gene,

• 1/2, reproduction/compensatory death per individual,

• λ, extinction/recolonization per colony,

• m_s, seed migration between colonies per individual.

In addition, two probabilities enter as parameters, namely,

• s, selfing probability,

• m_p, pollen migration probability between colonies during mating.

The corresponding set of parameters for the variable-size model are the Poisson intensities

• u, mutation of an allele into a new type per gene,

• 1/2, reproduction per individual,

• δ, death per individual,

• β, immigration of a new individual from the mainland per colony,

• m_s, seed migration between colonies per individual,

and the selfing probability s, 0 ≤ s ≤ 1. The reproduction intensity of 1/2 was retained to conform with the usual choice of normalization. In that case, with probability h/2 + o(h) one individual, and hence with probability h + o(h) one of two alleles in different individuals considered at time t + h, is the result of reproduction during (t, t + h].

Mutation acts in the same way in both models. It is the basic cause of reduction in identity by state of the genes in the population. The other mechanisms require more precise descriptions.

Reproduction: Obviously the main difference between fixed and variable size is that in the first case the Moran model mechanism of automatic compensating deaths following the formation of every new individual applies (Moran 1958), whereas in the latter a birth event more realistically increases the size of the population by one unit. In the variable-size case the assumption δ> 1/2 will guarantee that the population does not grow without bound and the assumption β > 0 will prevent permanent extinction.

To define reproduction in the fixed-size model we distinguish first of all whether reproduction occurs by random mating or asexually, and treat these cases separately. In both cases we assume the Moran model mechanism of automatic compensating deaths.

The asexual case is simpler. Here, in the event that the Poisson intensity signals that an individual is going to reproduce, a genetic copy of this diploid individual is produced followed by the death of one of the N original individuals. It is worth noting that the killed individual may thus be the one which itself triggered the production of the new.

The more complicated sexual reproduction case also involves the selfing parameter s as follows. Consider one of the N individuals and suppose its Poisson clock goes off. This we understand as an instruction to first make a copy of one of its own genes that is to form the first gene in the pair of genes representing the new individual. To find the second gene in the new pair there are two possibilities. Namely, with probability s the reproducing individual selects the second gene again to be a copy of one of its own, and with probability 1 – s it selects the second gene to be a copy of any one of the 2N existing genes sampled uniformly in the subpopulation. The actual selfing probability is therefore s + (1 – s)/ N rather than s. Again the formation in this manner of the new individual is instantly followed by a compensatory death.

In the variable-size model, for simplicity we restrict the analysis to sexual reproduction only.

Extinction/recolonization: As already stated we cover the two schemes introduced by Slatkin (1977) within the fixed-size model. The first mechanism is the so-called propagule recolonization, which means that upon extinction of a colony, k individuals are picked randomly in one of the n existing colonies to replace it. The second is the migrant pool. The difference is that now the replacement colony is put together by sampling k individuals uniformly over the whole population. Cases intermediate between those two extreme models are also taken care of by considering the probability of common origin of two colonists (Whitlock and McCauley 1990). In all cases, the k individuals then reproduce to restore the size of the population.

For the variable-size model one can say that deaths and immigration replace extinction/recolonization. Even if colonies may go extinct from time to time, influx from an external source is guaranteed. The effect of β is that external individuals enter into each colony independently, and an additional assumption will then be made regarding the genetic relationship between existing and immigrating gametes.

Migration: To each individual is associated the seed migration intensity m_s, which, when its Poisson clock rings, forces this individual to select randomly and uniformly one of the existing individuals in the population and then exchange position with the particular individual that was chosen. For example, in the fixed-size model this yields a probability 1/nN that nothing at all happens. Pollen migration, on the other hand, is a refinement of the random mating mechanism. We let m_p denote the probability that one of the genes composing a newly created individual due to pollen migration should be considered as being selected randomly in the full existing population and not only within its own colony.

Now we introduce more formally the IBD probabilities that we wish to study in this work. Pick two different alleles, A and B, randomly in the population. We consider the probabilities θ_k, k = 1, 2, 3, that A and B are identical by descent given their hierarchical “distance” k,

• θ₁, if A and B belong to the same individual,

• θ₂, if A and B belong to different individuals in the same colony,

• θ₃,if A and B are taken from different colonies in the population.

FIXED-POPULATION SIZE MODEL

Mutation and reproduction: We consider the IBD probabilities θ_k(t), t ≥ 0, k = 1, 2, 3, as functions of time, and investigate to what extent they change during a small interval of time (t, t + h], h > 0. Pick two alleles at hierarchical distance k from each other and observe them at time t + h.

First suppose that only mutation is in effect. Because the probability of a mutation in (t, t + h] of one of the selected genes is proportional to 2uh, θk(t+h)=2uh⋅0+(1−2uh)θk(t)+o(h); hence, θk(t+h)−θk(t)h=−2uθk(t)+o(h)h, and letting h tend to zero, θk′(t)=−2uθk(t). The solution θk(t)=θk(0)e−2ut,t≥0 (2) ultimately tends to zero as a result of the assumed infinite-allele model.

Now we include reproduction, which stabilizes the genetic structure in the population and counteracts the diversifying effect of mutation through genetic drift toward fixation. The balance of mutation and genetic drift then leads to stationary probabilities θ₁ and θ₂. However, there is no genetic drift over colony borders, and hence for k = 3 the only solution is still θ₃ = 0.

Asexual reproduction: Assume that all reproduction is asexual. Then the genetic relationship between two genes in a single individual does not change in case that individual was just formed as a copy of a previously existing individual. Mutation, on the other hand, will affect their interrelationship, forcing the IBD probability to vanish asymptotically as in (2). Hence the only stationary solution for θ₁(t) is θ₁ = 0.

For θ₂(t) note that with probability h + o(h) one of the chosen genes belongs to an individual that during (t, t + h] was produced by asexual reproduction. Then one of the two individuals involved, call it A, existed at time t and the other one, call it B, has been replaced by a copy of a randomly picked individual that existed at time t. With probability 1/N, B happens to be a copy of A in which case the two chosen genes must be considered belonging to the same individual. In this last case we have selected one gene in A and one gene in the copy of A called B. But the probability that two different genes in the same individual are IBD at t is θ₁(t), and thus our two genes are IBD with probability 1/2 + θ₁(t)/2. With probability 1 – 1/N their IBD probability remains at θ₂(t). Hence, θ2(t+h)=[h(1N(12+12θ1(t))+(1−1N)θ2(t))+(1−h)θ2(t)]⋅(1−2uh)+2uh⋅0+o(h), so that θ2(t+h)−θ2(t)=−2uhθ2(t)+h1N(12+12θ1(t)−θ2(t))+o(h) (3) and thus θ2′(t)=−2uθ2(t)+1N(12+12θ1(t)−θ2(t)). The stationary solution for θ₂(t) is obtained by taking θ′₂(t) = θ₁(t) = 0 above; thus θ2=12(1+2uN), which is half the value of the stationary solution of θ₂(t) in a haploid model.

Sexual reproduction: We start by deriving an equation for θ₂(t). Consider at time t + h two alleles, a and b, from two different individuals in the same colony and suppose one of them, a say, belongs to an individual that was created sexually in the same time interval elapsed since time t. With probability 1/N, allele a is copied from the individual to which b belongs, in which case they are identical by descent with probability 1/2 + θ₁(t)/2. With probability 1 – 1/N, the gene a was copied elsewhere; hence their IBD probability remains θ₂(t). Some reflection shows that these same relations hold whether or not the individual decided to self. Hence, θ2(t+h)=[h(1N(12+12θ1(t))+(1−1N)θ2(t))+(1−h)θ2(t)]⋅(1−2uh)+2uh⋅0+o(h). This is the same equation as the one obtained for the asexual case, thus again leading to θ2′(t)=−2uθ2(t)+1N(12+12θ1(t)−θ2(t)). (4) The corresponding equation for θ₁(t), however, is different. In fact, by splitting selfing from nonselfing, weighted by probabilities s and 1 – s, we find θ1(t+h)=[h(1−s)2(1N(12+12θ1(t))+(1−1N)θ2(t))+hs2(12+12θ1(t))+(1−h2)θ1(t)]⋅(1−2uh)+2uh⋅0+o(h), leading to θ1′(t)=−12(1+4u)θ1(t)+12θ2(t)+12(1−sN+s)(12+12θ1(t)−θ2(t)). (5) Relation (4) shows that at equilibrium θ2=1+θ12(1+2uN). (6) By inserting this into (5) taking θ′₁(t) = 0 we obtain the stationary solutions for mutation/genetic drift balance. They can be written θ1=(1+2uN⋅s(1−1∕N)−1∕N1+4u)/(1+2uN(2−s(1−1∕N)−1∕N1+4u)) (7) θ2=1∕(1+2uN(2−s(1−1∕N)−1∕N1+4u)). (8) Note that the case s=1∕N1−1∕N=1N−1, which corresponds to an actual selfing rate s + (1 – s)/N = 2/N, leads to the traditional expressions θ1=θ2=11+4uN.

Turning to Wright's indices, we obtain for the asexual case described previously, FIS=−11+4uN, and for the sexual system (7, 8), FIS=s(1−1∕N)−1∕N2(1+4u)−s(1−1∕N)+1∕N. (9) Negative values of F_IS, which occur for s < 1/(N – 1) in the sexual case, indicate that two genes in different individuals are more likely to be IBD than two genes in the same individual.

Adding migration and extinction: The next step toward enriching the applicability of the model and taking into account the colony structure that until now did not enter the derivations, is to add migration and extinction. Seed migration: To include seed migration we must consider the equations for θ₂(t) and θ₃(t), noting that migration of individuals does not affect θ₁(t). Consider first a pair of genes selected in the same colony. Then the intensity of a migration event in (t, t + h], which corresponds to selecting the two genes, in fact, from different colonies is 2m_s(1 – 1/n). Hence θ2(t+h)=2ms(1−1n)hθ3(t)+(1−2uh−2ms(1−1n)h)⋅[h1N(12+12θ1(t)−θ2(t))+θ2(t)]+o(h). (10) Similarly, if two genes are selected in different colonies at time t + h the intensity is 2m_sh/n such that they are considered to be stemming from the same colony at time t due to migration during (t, t + h]. Hence θ3(t+h)=θ3(t)−2(nu+ms)1nhθ3(t)=2ms1nhθ2(t)+o(h). Together with the previously derived Equation 5 for θ₁(t), these relations yield in the limit h → 0 the coupled system of equations θ1′(t)=−12(1+4u)θ1(t)+12θ2(t)+12(1−sN+s)(12+12θ1(t)−θ2(t)) (11) θ2′(t)=−2(u+ms(1−1n))θ2(t)+2ms(1−1n)θ3(t)+1N(12+12θ1(t)−θ2(t)) (12) θ3′(t)=−2(nu+ms)1nθ3(t)+2ms1nθ2(t). (13) By solving the system of Equations 11, 12, and 13 with vanishing left-hand terms for its stationary solution (θ₁, θ₂, θ₃) in steady state, we find the equilibrium relation θ2=1+θ12(1+2uαN), (14) which replaces (6), and then θ1=1+2uα(s(N−1)+1)(1+4u)(1+4uαN)+4u−2uα(s(N−1)+1)θ2=1+4u(1+4u)(1+4uαN)+4u−2uα(s(N−1)+1)θ3=msms+nuθ2, (15) where α=nms+nums+nu.

Pollen migration: Introduce the total migration parameter m=ms+12mp. To include migration of pollen, we start again with θ₁(t). During the mating phase the probability is m_p that one gene ends up in the new individual from a migrating pollen; hence the probability is m_p(1 – 1/n) that it originated in another colony. Therefore the term m_p(1 – 1/n)θ₃(t) enters the differential equation for θ₁(t), and the corresponding reduction occurs for the other, random mating terms, more exactly θ1′(t)=−12(1+4u)θ1(t)+12(1−mp(1−1n))θ2(t)+12(1−sN+s)(1−mp(1−1n))(12+12θ1(t)−θ2(t))+12mp(1−1n)θ3(t). (16) The differential equation for θ₂(t) modifies similarly. Exploiting the notation m introduced above we see that Equation 12 changes into θ2′(t)=−2(u+m(1−1n))θ2(t)+2m(1−1n)θ3(t)+1N(1−mp(1−1n))(12+12θ1(t)−θ2(t)). (17) Then we need to take into account that, because of pollen migration, reproductions in (t, t + h] now will affect θ₃. With probability m_p/n the selected genes belong to the same colony. This leads to θ3′(t)=−2(nu+m)1nθ3(t)+2m1nθ2(t)+mpn1N(12+12θ1(t)−θ2(t)). (18) In steady state the relation between θ₂ and θ₃ is found from (17) and (18), θ3=m+mpum+nu(1−mp(1−1∕n))θ2, (19) and the analog of system (15) is obtained similarly. In particular, it is seen that Equation 14 still holds in this case with α replaced by α′=nm+num+nu(1−mp(1−1∕n)).

To complete the analysis of the migration model, it remains to do the algebraic manipulations leading to the following expressions for Wright's indices: FIS=[2(u+m)g(s)+[(4ms+mp)(1−mp)⋅(1−1∕n)∕(1+4u)](1∕N)]∕[2(u+m)(2−g(s))−[(4ms+mp)(1−mp)(1−1∕n)∕(1+4u)](1∕N)] and FST=[1−mp]∕[(1−mp)(1−(4ms+mp)(1−1∕n)∕(1+4u))+2(u+m)N(2−g(s))], (20) where g(s)=1−mp(1−1∕n)1+4u(s(1−1N)−1N). In particular, for m_p = 0, FST=1∕[1−4ms(1−1∕n)(1+4u)+2(u+ms)N∕(2−[s(1−1∕N)−1∕N](1+4u))].

Extinction/recolonization: The recolonization event consists of two steps that are supposed to occur instantly in case of extinction. First k,1 ≤ k ≤ N, founder individuals, or colonists, are selected from a single colony (propagule model) or randomly over the n colonies (migrant pool model). The N – k additional individuals are produced by sexual reproduction from the k colonists. To avoid unnecessarily complicated expressions we restrict to the case m_p = 0, ignoring pollen migration.

The simplest case is the propagule model with k = N colonizers, in which an existing colony is copied and reinstalled in full to account for recolonization in the event of extinction. It is clear that extinction and the consecutive reinstallment of a colony have no effect upon θ₁ or θ₂. Hence we derive the stationary solution for θ₃(t), taking into consideration the possibility that one of the two colonies representing the pair of genes selected for the determination of θ₃(t + h) goes extinct during (t, t + h], the probability of which is 2λh. The probability is 1/n that the particular replacement colony is chosen to which the other gene of the given pair belongs, and the probability is 1/nN that even the same individual is chosen to which the other gene of the given pair belongs. Thus θ3(t+h)=θ3(t)−2(nu+λ+ms)1nhθ3(t)+2(λ+ms)1nhθ2(t)+2λ1n1Nh(12+12θ1(t)−θ2(t))+o(h) and so, replacing Equation 13, θ3′(t)=−2(nu+λ+ms)1nθ3(t)+2(λ+ms)1nθ2(t)+2λnN(12+12θ1(t)−θ2(t)). (21) The modified equation for θ₃(t) results in the equilibrium probabilities θ1=1+2uα(s(N−1)+1)(1+4u)(1+4uαN)+4u−2uα(s(N−1)+1)θ2=1+4u(1+4u)(1+4uαN)+4u−2uα(s(N−1)+1)θ3=λ+ms+2λu+2λms(1−1∕n)λ+ms+nu+2λms(1−1∕n)θ2, (22) where α=λ+nms+nuλ+ms+nu+2λms(1−1∕n). Moreover, FIS=[(λn+u+ms)(s(1−1N)−1N)+2ms(1−1n)⋅(1−2λn)1N]/[(λn+u+ms)(2(1+4u)−s(1−1N)+1N)−2ms(1−1n)(1−2λn)1N] and FST=(1−2λn)/[(1−2λn)(1−4ms(1−1∕n)1+4u)+2(λn+u+ms)⋅N(2−s(1−1∕N)−1∕N1+4u)]. (23) We observe that the case λ > n/2 is a nonbiological regime where the effect of recolonization is so strong that it is more likely for two genes to be identical by descent if they are picked from separate colonies rather than from the same. Restricting to λ ≤ n/2, we note that F_ST is decreasing as a function of λ, u, and m_s. The effect of increasing the number of colonies n is to reduce the effect of extinctions and enhance the effect of migration on F_ST.

Next we turn to the general propagule and migrant pool schemata with k colonizers. The Equations 11 for θ₁(t) and 21 for θ₃(t) with k = N remain unchanged, whereas Equation 12 for θ₂(t) changes. To find the modified equation for θ₂(t), suppose we pick at time t + h two alleles from two different individuals in a given colony. With probability λh + o(h), the colony suffers extinction and is recolonized by k founders, all this within (t, t + h]. Let ak=P(both alleles originate from the same colony-founder), noting a₁ = 1, a_N = 0. If the two alleles represent two different founders of their colony we make a further distinction by introducing ψ=P(two different colony founders originate from the same colony), so that in the propagule model we have ψ = 1 and in the migrant pool model ψ = 1/n.

Thus, temporarily switching off mutation, reproduction, and migration and considering only extinction/recolonization, θ2(t+h)=λhak(12+12θ1(t))+λh(1−ak)ψ(1N(12+12θ1(t))+(1−1N)θ2(t))+λh(1−ak)(1−ψ)θ3(t)+(1−λh)θ2(t)+o(h), and therefore θ2′(t)=−λ(1−ak)(1−ψ)θ2(t)+λ(1−ak)(1−ψ)θ3(t)+λ[ak+(1−ak)ψ1N](12+12θ1(t)−θ2(t)). In combination with (12) we find the new equation θ2′(t)=−[2u+(2ms(1−1n)+λ(1−ak)(1−ψ))]θ2(t)+(2ms(1−1n)+λ(1−ak)(1−ψ))θ3(t)+(λ[ak+(1−ak)ψ1N]+1N)(12+12θ1(t)−θ2(t)). (24)

The corresponding steady-state solutions (22) change into θ1=1+2uαk(s(N−1)+1)(1+4u)(1+4uαkN)+4u−2uαk(s(N−1)+1),θ2=1+4u(1+4u)(1+4uαkN)+4u−2uαk(s(N−1)+1),θ3=(λ+ms)B+2u+2ms(1−1∕n)+λ(1−ak)(1−ψ)(λ+ms+nu)B+2ms(1−1∕n)+λ(1−ak)(1−ψ)θ2, (25) where B=1λ+akN+(1−ak)ψ and αk=(λ+nms+nu)∕λ+n(1−ak)(1−ψ)(λ+ms+nu)B+2m(1−1∕n)+λ(1−ak)(1−ψ).

It remains to calculate the probabilities α_k, a task that can be rephrased as follows. Suppose we start with k individuals marked 1,..., k, and select from them by random sampling with replacement 2(N – k) allele copies, forming N – k new individuals. This scheme yields a colony of 2N genes classified in k categories that are represented by at least one gene each. Now select randomly with equal probabilities two individuals and in each of them a gene. We want to find the probability a_k that both genes belong to the same category (originate from the same founder). However, as long as the two chosen individuals do not both belong to the original group of size k, the probability is 1/k that they are marked with the same number. Hence, ak+(1−(k2)(N2))1k=1k−k−1N(N−1). A simpler scheme is to assume that a new colony consisting of N individuals is instantly sampled from the k founding individuals. Then ak=1∕k. (26)

The exact expressions for Wright's indices are somewhat unwieldy in this case, so we defer the further discussion to the final section where we apply some approximations and state a result.

IBD PROBABILITY IN VARYING POPULATION SIZE

An obvious drawback with the previous model is the restriction to fixed colony sizes N throughout the population. In population genetics modeling it has been natural to apply simple birth-death-immigration (BDI) models to widen the scope toward varying population size. In particular, there is a rich probabilistic structure in haploid populations based on the demographic decomposition into the number of families of allele types with a given number of representatives in the population (Tavaré 1989). Rannala and Hartigan (1995) explore some of these ideas for computing IBD probabilities in a metapopulation run by birth, death, and immigration. Rannala (1996) considers the corresponding sampling distribution of alleles under the same model and derives an expression for F_ST, which then coincides with θ₂ in our notation.

We take a different approach and extend the analysis of Rannala and Hartigan to include, under equilibrium conditions, selfing and migration between colonies. Our study may also be seen as a further alternative to the extinction/recolonization scheme in the fixed-size model.

A birth-death-immigration model: We start afresh from the situation that led to Equation 15, keeping seed migration but replacing the size-preserving Moran mechanism by death and immigration, as explained in generalities and scope.

Within this setting, δ introduces the potential for local population extinctions and β that for local recolonizations. Out of several possible schemes regarding the probability of IBD between a newly immigrated gene and an existent one, we assume that all immigrants are taken from a pool of individuals such that P(apair with one immigrant gene and one existing gene isIBD)=f≥0.

We assume that only seed migration is in effect. Hence individuals are exchanged between colonies without changing the population size. Let Ntj=the number of individuals at timetin colony nrj,1≤j≤n, and Ntpop=∑j=1nNtj=the total population size process. It is clear that N^pop_t is a continuous-time Markov birth-and-death process and that the components in the vector of colony sizes (N¹_t,..., Nⁿ_t) are dependent random processes.

The conditional jump probabilities at a time t, such that N^pop_t = r, are given by P(Nt+hpop−Ntpop=1∣Ntpop=r)=(nβ+r∕2)h+o(h),P(Nt+hpop−Ntpop=−1∣Ntpop=r)=δrh+o(h),h→0. This is a standard birth-and-death process, known as a BDI process with birth intensity ½, death intensity δ, and immigration intensity nβ. For each δ > ½, it has a negative binomial distribution as its stationary steady-state distribution. We let N^pop denote a random variable having this distribution. From properties of the BDI process, for large t, P(Ntpop=r)≈P(Npop=r)=(2nβ+r−1r)(1−12δ)2nβ(12δ)r. (27) Moreover, expectation and variance are given by E(Npop)=2nβ∕(2δ−1),V(Npop)=4nβδ∕(2δ−1)2.

Next we study the sizes of the single colonies. We refer to their stationary distributions by means of random variables (N¹,..., Nⁿ). By symmetry we immediately have E(Nj)=2β∕(2δ−1),1≤j≤n. Furthermore, in the special case m_s = 0 the n components N^j are obviously independent and thus each N^j possesses a negative binomial distribution with parameters 2β and 2δ. In fact, their sum is then also negative binomial with parameters 2nβ and 2δ, in agreement with (27). However, we do want to consider the case m_s > 0 where the colony sizes are no more independent.

To this end, we note that in the present framework migration corresponds to the death of an individual in one colony followed instantly by immigration into a randomly chosen colony. Hence, focusing on a specific colony size, N^j_t say, the total external immigration intensity at time t is given by β+ms1n∑i≠jNi≈β+ms(1−1n)E(Nj), because a migrating gene from any other colony will end up in colony nr j with probability 1/n. The approximation should be valid at least for large values of n in view of the law of large numbers. Hence we introduce for each component N^j_t a local immigration intensity. β∼=β+ms(1−1n)E(Nj)=β+2βms(1−1∕n)2δ−1, (28) and, similarly, a local death intensity by δ∼:=δ+ms(1−1∕n)>1∕2. Indeed, a migrating gene from colony nr j is going to enter a different colony with probability 1 – 1/n.

As an approximation for the colony size we have thus found a slightly different BDI process. Its conditional jump probabilities given that N^j_t = r are P(Nt+hj−Ntj=1∣Ntj=r)=(β∼+r∕2)h+o(h),P(Nt+hj−Ntj=−1∣Ntj=r)=δ∼rh+o(h),h→0. Consistent with our previous results the corresponding average colony size at stationarity is again E(Nj)=2β∼2δ∼−1=2β2δ−1.

It is interesting to compare the variance for the components N^j with the variance of their sum N^pop in the dependent case m > 0. We obtain V(Nj)=4nβ∼δ∼∕(2δ∼−1)2=2β2δ−12δ+2ms(1−1∕n)2δ−1+2ms(1−1∕n). It is a simple observation that V(N^j) is maximal for m_s = 0 [where of course V(N^pop) = nV(N¹)] and decreases toward its minimum value, coinciding with the mean E(N^j), as the migration rate m is allowed to increase. In this sense migration is seen to stabilize the system and counteract fluctuations in the colony sizes.

Probabilities of IBD in the varying-population-size model: Because in the present model with positive probability some colonies or even the whole population may be extinct from time to time, the definition of the IBD probabilities θ_k is not as straightforward as in the previous constant-size model. To define θ₂ there must be at least one colony with two individuals, and to define θ₃ at least two nonempty colonies. It is reasonable to argue that otherwise corresponding weights should be assigned to let θ₃ = 0 and θ₂ = 0, respectively. However, we have, for example, P(Npop=0)=(1−12δ)2nβ≈e−nβ∕δ, (29) which suggests that we simply ignore this effect. Hence to derive a balanced differential equation for the IBD probability within a colony we assume that at some time t a colony j has size Ntj=r , where r ≥ 3, and we let θ₂(t) refer to the probability that two genes from different individuals in this colony are identical by descent.

Now we go over the various events that may occur in (t, t + h]. First of all, with intensity δ∼r there was a death leading to r – 1 ≥ 2 individuals in the selected colony at time t + h. However, the removal of r – 1 such potential pairings does not change the probability of IBD for the remaining pairs (cf. Appendix in Rannala and Hartigan 1995).

Next, with intensity β∼ an individual immigrated between t and t + h, increasing the colony size from r to r + 1. With probability (r1)(r+12)=2r+1 a gene from this newly immigrated individual is selected to the pairing needed to determine θ₂. If the colony immigrant is an external immigrant to the population this results in the probability f for the pair to be identical by descent. If the colony immigrant migrated from another colony the gene pair is related according to θ₃. Otherwise, if the immigrant was not chosen for the pair, the probability for IBD remains the same as before immigration.

Finally, with intensity r/2 a birth occurs in (t, t + h] such that again the size changes from r to r + 1. As a pair of genes A and B are selected at time t + h, the probability is 2/(r + 1) that one of them belongs to the newly reproduced individual. Let us suppose that A does. Then we have two possibilities. With probability 1/r gene A was produced by random sampling from the individual that gene B belongs to, in which case the probability of IBD is 1/2 + θ₁(t)/2. With probability 1 – 1/r the relationship between A and B is given simply by θ₂(t).

Summing up, θ2(t+h)=h2r+1(βf+2msβ2δ−1(1−1n)θ3(t))+β∼h(1−2r+1)θ2(t)+12rh{2r+1[1r(12+12θ1(t))+(1−1r)θ2(t)]+(1−2r+1)θ2(t)}+(1−2uh−β∼h−rh∕2)θ2(t)+o(h). Hence, θ2′(t)=2βr+1f+4βms(1−1∕n)(r+1)(2δ−1)θ3(t)+1r+1(12+12θ1(t)−θ2(t))−2(u+β∼r+1)θ2(t). By a similar argument θ1′(t)=12rr+1[(s+1−sr)(12+12θ1(t)−θ2(t))+θ2(t)]−12(rr+1+4u)θ1(t).

Turning to θ₃ suppose two colonies i and j have been chosen. We must determine the effect of immigration into any one of them during (t, t + h]. Suppose colony i of size r is subject to external immigration (or colony j of size q). Then with probability 1/(r + 1) the immigrant belongs to the select pair resulting in IBD probability f. Similarly, if the immigrant is a migrant from j, then again with probability 1/(r + 1) the IBD probability of the select pair changes, this time into θ₂. The intensities are for the first event β and for the second approximately mE(N^j)/n. Hence θ3′(t)=(1r+1+1q+1)(fβ+1nmsE(Nj)θ2(t))−2uθ3(t)−(1r+1+1q+1)(1nmsE(Nj)+β)θ3(t).

At this point we observe that the steady-state probabilities corresponding to the balanced equations we have derived will depend on r and q unless u = 0. In the present model the external immigration parameter β has a similar effect and purpose as mutation, namely, to guarantee fresh import of genetic material and enhance genetic diversity. The difference is that mutation acts only on the individual level while immigration occurs in colonies regardless of their sizes. We therefore make the approximation in the sequel to ignore mutation and simply restrict the case to u = 0. Moreover, because θ′₁(t) = 0 implies (s+1−sr)(12+12θ1−θ2)+θ2−θ1=0, and thus θ1=1−(1−s)(1−1∕r)+2(1−s)(1−1∕r)θ21+(1−s)(1−1∕r)≈s+2(1−s)θ22−s for large r, which is again the consistency relation (34), we are led to study the balanced equations θ∼1=s+2(1−s)θ∼22−s,2βf+2ms(1−1n)E(Nj)θ∼3+12+12θ∼1−θ∼2=2β∼θ∼2,2βf+2msE(Nj)nθ∼2=2(msE(Nj)∕n+β)θ∼3. The solution is given by θ∼1=1+2sαβ+4(2−s)αβf1+2(2−s)αβ,θ∼2=1+2(2−s)αβf1+2(2−s)αβ,θ∼3=(2δ−1)f+2m2θ∼2∕n2ms∕n+2δ−1, (30) where now α=β+msE(Nj)β+msE(Nj)∕n=2ms+2δ−12ms∕n+2δ−1.

In conclusion we have found, under the simplifying assumptions of mutation being negligible and the typical size of a colony being large, that the conditional IBD probabilities at equilibrium for given colony sizes are expressed by the relations (30). Again relying on, e.g., (29), we let expression (30) represent the unconditional IBD probabilities for the variable-size model. Thus Wright's genetic indices for the BDI model turn out to be FIS=s2−s and FST=1∕(1+2β(2ms+2δ−12δ−1)(2−s)). (31) Interestingly, the last quantity depends on neither f nor n.