Diffusion approximation of a multilocus model with assortative mating

To understand the effect of assortative mating on the genetic evolution of a population, we consider a finite population in which each individual has a type, determined by a sequence of n diallelic loci. We assume that the population evolves according to a Moran model with weak assortative mating, strong recombination and low mutation rates. With an appropriate rescaling of time, we obtain that the evolution of the genotypic frequencies in a large population can be approximated by the evolution of the product of the allelic frequencies at each locus, and the vector of the allelic frequencies is approximately governed by a diffusion. We present some features of the limiting diffusions (in particular their boundary behaviour and conditions under which the allelic frequencies at different loci evolve independently). If mutation rates are strictly positive then the limiting diffusion is reversible and, under some assumptions, the critical points of the stationary density can be characterised.


Introduction
The aim of this paper is to construct and analyse a diffusion approximation for a diallelic multilocus reproduction model with assortative mating, recombination and mutation. Our starting point is a variant of the Moran model. We suppose that the population is monoecious 1 , haploid 2 and of constant size N . This will be an overlapping generation model, but, in contrast to the usual Moran framework, we suppose that reproduction takes place at discrete times 1, 2, . . . . In each time step, a mating event occurs between two individuals I 1 and I 2 ; I 1 is replaced by an offspring, so that the size of the population is kept constant. The genotype of the offspring is obtained from those of I 1 and I 2 through a process of recombination followed by mutation which we make precise in §2. In the classical Moran model, the two individuals I 1 and I 2 are chosen at random from the population. Here, to study the effects of assortative mating, we assume that the first individual, I 1 , is still chosen at random, but the second individual, I 2 , is sampled with a probability that depends on its genotype and on the genotype of the first selected individual. The genotype of an individual is composed of a finite number, n, of loci with two alleles per locus denoted by 0 and 1. To characterise the assortative mating, we introduce a real parameter s i,j for every pair of genotypes (i, j). If I 1 has genotype i, then, in the draw of I 2 , an individual with genotype j has a probability proportional to 1 + 1 N s i,j of being selected.
Diffusion approximations for different selection-mutation models have been studied extensively in the one-locus case (see, for example, Ethier & Kurtz 1986, Chapter 10). The coefficients 1 + 1 N s i,j of our model play the same rôle as the (viability) selection coefficients in a Wright-Fisher model for a diploid population. Since they depend on the types of both parents they result in nonlinear (frequency dependent) selection (see §4.4). Ethier & Nagylaki (1989) study two-locus Wright-Fisher models for a panmictic 3 , monoecious, diploid population of constant size under various assumptions on selection and recombination. Depending on the strength of the linkage between the two loci, they obtain different types of diffusion approximation: limiting diffusions for gametic 4 frequencies if the recombination fraction multiplied by the population size tends to a constant as the size tends to +∞ (so-called tight linkage) 1 Introduction 3 and limiting diffusions for allelic frequencies if the recombination fraction multiplied by the population size tends to +∞ (so-called loose linkage). To our knowledge, this work has not been extended to more general multilocus models with recombination and assortative mating.
Nevertheless, there is a large body of work on multilocus genetic systems. Most theoretical investigations assume that the size of the population is infinite, so that the random genetic drift can be ignored; the evolution of genotypic frequencies is then described by recursive equations or by differential equations (see Christiansen 2000 and references within). A comparison between infinite and finite population models with random mating is presented in Baake & Herms (2008). A review of several simulation studies can be found in the introduction of Devaux & Lande (2008). Among these, the 'species formation model', introduced by Higgs & Derrida (1992) inspired our work. In their model, mating is only possible between individuals with sufficiently similar genotypes, so that from the point of view of reproduction the population is split into isolated subgroups. Their simulations display a succession of divisions and extinctions of subgroups. In this paper we generalise their assortative mating criterion to one defined through the family of parameters s i,j and provide a general theoretical treatment.
To give an overview of our results, we first consider a particular pattern of assortative mating. Let us assume that the frequency of matings between two individuals of types i and j depends only on the number of loci at which their allelic types differ (and not on the positions of those loci along the genome). We then have a model with n + 1 assortment parameters, denoted by s 0 , . . . , s n , obtained by setting s i,j = s k if the genotypes i, j are different at exactly k loci (regardless of their positions). This mating criterion will be called the Hamming criterion in what follows. A decreasing sequence s 0 ≥ s 1 ≥ . . . ≥ s n will describe a positive assortative mating (individuals mate preferentially with individuals that are similar). An increasing sequence s 0 ≤ s 1 ≤ . . . ≤ s n will describe a negative assortative mating (individuals mate preferentially with individuals that are dissimilar).
We establish a weak convergence of the Markov chain describing the genetic evolution of the population as its size tends to +∞, under a hypothesis on the recombination distribution that corresponds to loose linkage (during each reproduction event, recombination between any pair of loci occurs with a positive probability) and under the assumption that mutations occur independently at each locus with the same rates (at each locus, the rate of mutation of a type 0 allele to a type 1 allele is µ 0 N and the rate of mutation of a type 1 to a type 0 is µ 1 N ). In particular, while mutation and assortment parameters are rescaled with population size, recombination is not. As a result, we see a separation of timescales. Due to recombination, the genotypic frequencies rapidly converge to a product distribution which is characterised by its marginals, that is by the 0-allelic frequencies at each locus. We show that, at a slower rate, the allelic frequencies converge to a multidimensional diffusion, whose components are coupled only through an infinitesimal drift term (in the mathematical sense) arising from the assortative mating.
Let us describe some features of the limiting diffusion. If s 1 − s 0 = s 2 − s 1 = . . . = s n − s n−1 then the frequencies of the 0-allele at each locus evolve according to independent Wright-Fisher diffusions with mutation rates µ 0 and µ 1 and symmetric balancing selection with strength 1 2 (s 1 − s 0 ); that is they solve the following stochastic differential equation: where (W t ) t≥0 is a standard Brownian motion. In all other cases, the allelic frequencies at different loci no longer evolve independently. Instead the vector of 0-allelic frequencies (x t (1), . . . , x t (n)) is governed by the stochastic differential equation: where (W t (1)) t≥0 ,. . . , (W t (n)) t≥0 denote n independent standard Brownian motions and P (x (i) ) is a symmetric polynomial function of the n − 1 variables x(j)(1 − x(j)), j ∈ {1, . . . , n} \ {i} whose coefficients depend only on the parameters s 1 − s 0 , . . . , s n − s n−1 . More precisely, P (x (i) ) is an increasing function of each parameter s 1 − s 0 , . . . , s n − s n−1 (see Theorem 4.1 for an explicit formula of P (x (i) )). When the mutation rates µ 0 and µ 1 are strictly positive, the limiting diffusion has a reversible stationary measure, the density of which is explicit. When the two mutation rates are equal to µ > 0, we describe the properties of the critical points of the density of the stationary measure. In particular, we find sufficient conditions on µ and s 1 − s 0 , . . . , s n − s n−1 for the state where the frequencies of the two alleles are equal to 1/2 at each locus to be a global maximum and for the stationary measure to have 2 n modes. These sufficient conditions generalise the independent case. For example, when µ > 1/2 they imply the following results: 1. if s ℓ − s ℓ−1 ≥ −(8µ − 4) for every ℓ ∈ {1, . . . , n}, then (1/2, . . . , 1/2) is the only mode of the stationary measure; 2. if s n − s n−1 ≤ . . . ≤ s 1 − s 0 < −(8µ − 4) then the stationary measure has 2 n modes.
These results can be extended to other patterns of assortative mating. In fact, we need only make the following assumption on the parameters s i,j : the value of the assortment parameter s i,j between two genotypes i and j is assumed to be the same as the value of s j,i and to depend only on the loci at which i and j differ. In particular this implies that the value of s i,i is the same for every genotype i. This generalises the Hamming criterion and allows us to consider more realistic situations in which the influence on mating choice differs between loci (see §2.3). It transpires that, under these assumptions, the limiting diffusion does not depend on the whole family of assortment parameters, but only on one coefficient per subgroup of loci L. We denote this coefficient m L (s). It is the mean of the assortment parameters for pairs of genotypes that carry different alleles on each locus in L and identical alleles on all other loci.
The stochastic differential equation followed by the limiting diffusion can still be described by equation (1.1) if the symmetric polynomial term P (x (i) ) in the drift of the i-th coordinate is replaced by a non-symmetric polynomial term P i (x (i) ) in the coefficients of which the quantities The rest of the paper is organized as follows. In §2, we present our multilocus Moran model.
In §3, we describe the diffusion approximation for the one-locus model and compare it with a diffusion approximation for a population undergoing mutation and 'balancing selection'. We recall some well-known properties of this diffusion, in particular the boundary behaviour and the form of the stationary measure, for later comparison with the multilocus case. In §4, we state our main result concerning convergence to a diffusion approximation in the multilocus case (Theorem 4.1) and give two equivalent expressions for the limiting diffusion. We then compare with the two-locus diffusion approximation obtained in Ethier & Nagylaki (1989).
The proof of Theorem 4.1 is postponed until §7. In §5, we derive some general properties of the limiting diffusion. §6 is devoted to the study of the density of the stationary measure. An appendix collects some technical results used in the description of the limiting diffusion.

The discrete model
This section is devoted to a detailed presentation of the individual based model. The assumptions on assortative mating, recombination and mutation that we will require to establish a diffusion approximation for the allelic frequencies are discussed at the end of the section.

Description of the model
We consider a monoecious and haploid population of size N where the type of each individual is described by a sequence of n diallelic loci. For the sake of brevity, let the set of loci be identified with the set of integers 1 ; n := {1, . . . , n} and let the two alleles at each locus be labelled 0 and 1. The type of an individual is then identified by an n-tuple k := (k 1 , ..., k n ) ∈ {0, 1} n . Let A = {0, 1} n be the set of possible types. The proportion of individuals of type k at time t ∈ IN will be denoted by Z (N ) t (k) so that the composition of the population is described by the set Z At each unit of time the population evolves under the effect of assortative mating, recombination and mutation as follows.
Assortative mating: at each time t, two individuals are sampled from the population in such a way that: 1. the first individual has probability Z (N ) t (i) of being of type i; 2. given that the first individual chosen is of type i, the probability that the second individual is of type j is 1 + s where the assortment parameters {s (N ) i,j , i, j ∈ A} are fixed nonnegative real numbers 5 .
The population at time t+1 is obtained by replacing the first chosen individual with an offspring whose type is the result of the following process of recombination followed by mutation.
Recombination: for each subset L of 1 ; n , let r L denote the probability that the offspring inherits the genes of the first chosen parent at loci ℓ ∈ L and the genes of the second parent at loci ℓ ∈ L. The family of parameters {r L , L ⊂ 1 ; n } defines a probability distribution, called the recombination distribution, on the power set P( 1 ; n ) (it was first introduced in this manner by Geiringer (1944) to describe the recombination-segregation of gametes in a diploid population). It is natural to assume that the two parents contribute symmetrically to the offspring genotype, that is: Assumption H1: for each subset L of 1 ; n , r L = rL whereL denotes the complementary set of loci, 1 ; n \ L.
With this notation, the probability that, before mutation, the offspring of a pair of individ- Let us express some classical examples of recombination distributions in this notation: Examples 2.1.
Finally we superpose mutation.
Mutation: we assume that mutations occur independently and at the same rate at each locus: µ (N ) 1 will denote the probability that an allele 1 at a given locus of the offspring changes into allele 0 and µ (N ) 0 the probability of the reverse mutation. The probability that the mutation process changes a type k into a type ℓ is:

Expression for the transition probabilities
It is now elementary to write down an expression for the transition probabilities of our model. In the notation above, if z = {z(k), k ∈ A} describes the proportion of individuals of each type in the population at a given time, then the probability that, in the next time step, the number of individuals of type j increases by one and the number of individuals of type i = j decreases by one is

Assumptions on assortative mating, recombination and mutation
In order to obtain a diffusion approximation for a large population, we assume that mutation and assortment parameters are both O(N −1 ), so we set Just as in the two-locus case studied by Ethier & Nagylaki (1989), we can expect diffusion approximations to exist under two quite different assumptions on recombination, corresponding to tight and loose linkage. Here we focus on loose linkage. More precisely, we assume that the recombination distribution does not depend on the size of the population and that recombination can occur between any pair of loci: Assumption H3: For every I ∈ P( 1 ; n ), r I does not depend on N and for any distinct integers h, k ∈ 1 ; n , there exists a subset I ∈ P( 1 ; n ) such that h ∈ I, k ∈ I and r I > 0.
This assumption is satisfied for the last two examples of recombination distribution presented in Example 2.1, but not in the absolute linkage case. In infinite population size multilocus models with random mating, and without selection, this condition is known to ensure that in time the genotype frequencies will converge to linkage equilibrium, where they are products of their respective marginal allelic frequencies (see Geiringer 1944 andNagylaki 1993 for a study of the evolution of multilocus linkage disequilibria under weak selection).
In order that the generator of the limiting diffusion has a tractable form, we shall make two further assumptions on the family of assortment coefficients s = {s i,j , (i, j) ∈ A 2 }: Assumption H4: for every (i, j) ∈ A 2 , 1. s i,j = s j,i 2. the value of s i,j depends only on the set of loci k at which i k = 0 and j k = 1 and on the set of loci ℓ at which i ℓ = 1 and j ℓ = 0.
These conditions mean that the probability of mating between two individuals at a fixed time depends only on the difference between their types. In particular, two individuals of the same type will have a probability of mating that does not depend on their common type: s i,i = s j,j for every i, j ∈ A. In the one-locus case, this assumption means that the model distinguishes only two classes of pairs of individuals since s 0,1 = s 1,0 and s 0,0 = s 1,1 .
To describe positive or negative assortative mating we have to choose how to quantify similarities between two types. Let us present two criteria that provide assortment parameters for which assumption H4 is fulfilled: 1. Hamming Criterion. One simple measure to quantify similarities between two types is the number of loci with distinct alleles: s i,j will be defined as nonnegative reals depending only on the Hamming distance between i and j denoted by d h (i, j) := n ℓ=1 |i ℓ − j ℓ |. A positive assortative mating will be described by a sequence of n + 1 nonnegative reals s 0 ≥ s 1 ≥ . . . ≥ s n by setting s i,j = s d h (i,j) for every i, j ∈ A. This criterion will be called Hamming criterion.
2. Additive Criterion. If we assume that the assortment is based on a phenotypic trait which is determined by the n genes whose effects are similar and additive, then a convenient measure of the difference between individuals of type i and j is d a (i, j) := | n ℓ=1 (i ℓ − j ℓ )|. A positive assortative mating will be described by a sequence of n + 1 nonnegative reals s 0 ≥ s 1 ≥ . . . ≥ s n by setting s i,j = s da(i,j) for every i, j ∈ A. This criterion will be called additive criterion.
The assortative mating in the species formation model of Higgs & Derrida (1992) is a special case of the Hamming criterion. The additive criterion is widely used in models in which assortative mating is determined by an additive genetic trait. is qualitatively similar to that observed in the simulations of Higgs & Derrida (1992) for the Hamming criterion, namely continuous creation of reproductively isolated subgroups.
With the Hamming and additive criteria, every locus is assumed to have an identical positive or negative influence on the assortment. As we have defined a general family of assortment parameters, it is possible to consider more complex situations. For instance, we can take into account that some loci have a greater influence on the mating choice than others by dividing the set of loci into two disjoint subgroups G 1 and G 2 ; we introduce two sets of assortment parameters s (1) and s (2) that satisfy assumption H4 for the subgroups of loci G 1 and G 2 respectively. If we assume that the effects of the two subgroups are additive we set s i,j = s for every i, j ∈ A. This defines a set of assortment parameters that satisfies assumption H4.
More generally, any set of assortment parameters defined as a function of s (1) and s (2) satisfies assumption H4.

The one-locus diffusion approximation
Before studying the multilocus case, for later comparison, in this section we record some properties of the one-locus model.

The generator of the one-locus diffusion
In the case of one locus (n = 1), under assumption H2, the frequency of 0-alleles satisfies: Therefore the distribution of the frequency of 0-alleles at time [N 2 t] is approximately governed, when N is large, by a diffusion with generator: More precisely, if Z If we assume that s satisfies assumption H4, that is s 0,0 = s 1,1 and s 0,1 = s 1,0 , and if we denote their common values by s 0 and s 1 respectively, then the drift has a simpler form and we obtain Remark 3.1. This diffusion can also be obtained as an approximation of a diploid model with random mating, mutation and weak selection in favour of homozygosity 6 (when s 0 − s 1 > 0) or in favour of heterozygosity (when s 0 − s 1 < 0) (see, for example, Ethier & Kurtz 1986, Chapter 10).

Properties of the one-locus diffusion
Stationary measure. If µ 0 and µ 1 are strictly positive, this diffusion has a reversible stationary measure. Its density with respect to Lebesgue measure on [0, 1] is given by Wright's formula: where the constant C µ,s is chosen so that (ii) if µ 1 ≥ 1/2 then 0 is an entrance boundary (started from a point in ]0, 1[ the diffusion will not reach 0 in finite time, but the process started from 0 is well-defined); (iii) if 0 < µ 1 < 1/2 then 0 is a regular boundary (starting from a point z 0 ∈]0, 1[ the diffusion has a positive probability of reaching 0 before any point b ∈]z 0 , 1] in a finite time and the diffusion started from 0 is well-defined); with the obvious symmetric definitions at 1.

Convergence to a diffusion in the multilocus case
In the case of several loci, under assumptions H2 and H3, a Taylor expansion shows that = z] is of order 1 N 2 only inside the set of product distributions on {0, 1} n . This set is often called the Wright manifold or the linkage-equilibrium manifold and denoted by W n (a population is said to be in linkage equilibrium if the geno- . . , k n ) denotes the frequency of individuals having the allele x at the j-th locus).
Outside this manifold, the drift pushes the process towards the Wright manifold at an exponential speed. Therefore to extend the diffusion approximation to the n-locus case, we introduce a change of variables composed of the n 0-allelic frequencies and of 2 n − n − 1 processes that measure the deviation from the linkage equilibrium.
For a nonempty subset L of {1, . . . , n}, between the loci in L at time t. (This is just one of many ways to measure the linkage disequilibrium, see for example Bürger (2000), Chapter V.4.2, for other measures.) The vector of 0-allelic frequencies at time t is X We shall show that if t N tends to +∞ faster than N then Y (N ) [t N ] converges to 0 while if time is sped up by N 2 then X (N ) converges to a diffusion as N tends to +∞.
Before giving a precise statement of the convergence result for the two processes X (N ) and Y (N ) (Theorem 4.1), let us introduce some notation in which to express the parameters of the limiting diffusion.

Mean assortment parameters
For a subset L of loci, consider the set of pairs of genotypes that differ at each locus ℓ ∈ L and are equal at each locus ℓ / ∈ L: Let m L (s) denote the mean value of the assortment parameters for all pairs in this set F L : Examples 4.1.
In each of these expressions the four coefficients are equal by assumption H4.
In this expression the first (resp. last) two coefficients are equal by H4.

Convergence to a diffusion
The following theorem provides convergence results for the two processes X (N ) and Y (N ) as the population size N tends to +∞. The proof, based on Theorem 3.3 of Ethier & Nagylaki (1980), is postponed to §7.
and its closure is the generator of a strongly continuous semigroup of contractions.
) to a diffusion process X with generator G n,s where the polynomial function P i,s (x) has the following expression: Remark 4.1.
1. The recombination distribution (r I ) I⊂ 1;n does not appear in the expression for the limiting diffusion. Nevertheless, the proof of Theorem 4.1 will show that it has an influence on the speed of convergence of the linkage disequilibrium to 0.
2. The limiting diffusion depends on the assortment parameters only via the quantities m A (s) for every A ⊂ 1 ; n . A set of assortment parameters for which m A∪{i} (s) − m A (s) < 0 for every i ∈ 1 ; n and A ⊂ 1 ; n \ {i} favours homozygous mating with respect to the genotype at the i-th locus. It is therefore no surprise that by increasing the value of m A∪{i} (s) − m A (s) for a fixed subset A, we increase the value of the i-th coordinate of the drift at a point x for which x i < 1/2 and decrease it at a point x for which x i > 1/2.

Another expression for the polynomial term
k ∈ 1 ; n \ {i} yields the following expression: The details of the proof are provided in §7.2.
The coefficients α i,L (s) can be compactly expressed using difference operators. Let us introduce some notation: for a function f defined on the subsets of a finite set E and for an Since δ i • δ j = δ j • δ i for every i, j ∈ E, we can, more generally, introduce a difference operator proof by induction on |B| provides the following formula for δ B : Let m(s) denote the function A → m A (s) defined on the subsets of 1 ; n . In this notation, for If, for each subset A of loci, the coefficient m A (s) depends only on the number of loci in A, then it follows from expression (4.3) that P i,s (x) is a symmetric polynomial function, the coefficients of which do not depend on i. This is the case for instance with the Hamming and additive criteria (see Example 4.1 for the corresponding expressions for m A (s)). Let us give the expanded form of P i,s (x) for the Hamming criterion: . As in the general case, the coefficientα k (s) has a compact expression in terms of difference operators. Let δ (1) denote the forward difference operators: δ (1) [s](i) = s i+1 − s i for every i ∈ 0 ; n − 1 . The forward difference operators of higher orders are defined iteratively: Ethier & Nagylaki (1989) established convergence results for a general multiallelic twolocus Wright-Fisher model of a panmictic, monoecious, diploid population of N individuals (identified with 2N haploids) undergoing mutation and selection. In their model, a gamete is described by a pair i = (i 1 , i 2 ) ∈ 1 ; r 1 × 1 ; r 2 where r 1 is the number of alleles in the first locus and r 2 is the number of alleles in the second locus. The parameters of their model are:

Comparison with the two-locus Wright-Fisher diffusion
1. the viability of a pair of gametes (i, j) denoted by w N,i,j = 1− σ N,i,j with the assumption σ N,i,j = σ N,j,i and σ N,i,i = 0 for every i, j ∈ 1 ; r 1 × 1 ; r 2 (after viability selection the proportion of a pair of gametes (i, j) is assumed to be P denotes the frequency of gametes k in the population ∀k ∈ 1 ; r 1 × 1 ; r 2 ); 2. the recombination fraction c N ; j,k that the j-th allele in the i-th locus mutates to the k-th allele.
The population at the generation t + 1 is obtained by choosing 2N gametes uniformly at random with replacement from the pool of gametes of the generation t after the steps of viability selection, recombination and mutation.
They studied the diffusion approximation under several assumptions on selection and recombination coefficients. In the case of weak selection (2N σ N,i,j converges to a real number denoted by σ i,j for every i, j) and loose linkage (c N converges to a finite limit and N c N tends to +∞) they obtained a limiting diffusion for the allelic frequencies (p 1 , . . . , p r 1 −1 , q 1 , . . . , q r 2 −1 ) of the alleles 1, . . . , r 1 − 1 in the first locus and the alleles 1, . . . , r 2 − 1 in the second locus. In the case of two alleles at each locus (r 1 = r 2 = 2), the generator of the limiting diffusion is Accordingly, the generator L coincides with G 2,s if we assume (a) that the mutation rates ν and set σ i,j = − 1 2 s i−1,j−1 , for every i, j ∈ {1, 2} 2 (with the notation 1 = (1, . . . , 1)). This comparison suggests that the effect of assortative mating on the genotype evolution of a large population in our model is similar to the effect of weak viability selection in a diploid Wright-Fisher model with mutation.

Description of the limiting diffusion
This section collects some properties that can be deduced from the form of the generator, G n,s , of the limiting diffusion.

The set of generators arising from the model
where {α A , A ⊂ 1 ; n , A = ∅} is a family of real numbers, can be interpreted as the generator of the diffusion approximation of an n-locus Moran model as defined in §2.
Proof. We may, for instance, take the following set of assortment parameters {s i,j , i, j ∈ A}: • s i,i = 0 for every i ∈ A.
Let us check that this family satisfies 2 |L|−1 δ L [m(s)](∅) = α L for every nonempty subset L of We invert the double sum and use the formula In particular, the n-locus Moran model with assortative mating based on the Hamming criterion allows us to obtain, through diffusion approximation, any generator on C 2 ([0, 1] n ) of the form: To see this, given any sequence α 0 , . . . , α n−1 of n reals, we have to find n + 1 real numbers s 0 , . . . , s n such that α ℓ = 2 ℓ δ (ℓ+1) [s](0). These are given by the inversion formula (A.3) in the Appendix, from which we see that we may set s 0 = 0 and s k = k ℓ=1 2 1−ℓ k ℓ α ℓ−1 for k ∈ 1 ; n .

The generator for two groups of loci
Let us consider a partition of the set of loci into two subgroups, G 1 = 1 ; k and G 2 = k + 1 ; n , say. We introduce two sets of assortment parameters s (1) and s (2) depending on subgroups of loci from G 1 and from G 2 respectively and satisfying assumption H4.
If we assume that the assortment parameter between two individuals of type i and j is for every i, j ∈ A, then m L (s) = m L∩G 1 (s (1) ) + m L∩G 2 (s (2) ) for every subset L of 1 ; n . This implies that the first k coordinates of diffusion limit evolve independently of the last n − k coordinates and that the generator of the diffusion limit is: Therefore, with these choices we can reduce our study to subgroups of loci having the same influence on assortment.

Conditions for independent coordinates
For some patterns of assortment, the allelic frequencies at each locus in a large population evolve approximately as independent diffusions: Proposition 5.1. Assume that the assortment parameters s = {s i,j , i, j ∈ A} satisfy the assumption H4. 2. If condition (H5) holds, the i-th coordinate behaves as the one-locus diffusion with assortment coefficients s 0 = s 1,1 and s 1 = s u i ,1 where u i = (0 {i} , 1 1;n \{i} ) denotes the genotype which differs from the genotype 1 only on the locus i; its generator is 3. In particular, Proof. First note that G n,s is the generator of n independent diffusions if and only if the polynomial term P i,s (x) is a constant function for every i ∈ 1 ; n .
1. According to the formula (4.2), the polynomial term P i,s (x) is a constant function for every i ∈ 1 ; n whenever condition H5 holds. Conversely, assume that the polynomial term P i,s (x) is a constant function for every i ∈ 1 ; n . By formulae (4.3) and (4.6), δ L [m(s)](∅) = 0 for every subset L of 1 ; n having at least two elements. We derive, from the inversion formula (A.2) stated in the Appendix, that for every subset A ∈ P( 1 ; n ), Therefore, condition H5 is satisfied.
After some computation, we obtain for i ∈ L, It follows from (5.1) that for every c ∈ IR, the system defined by

Behaviour at the boundaries
In this section the trajectories of the coordinates of the limiting diffusion are compared with those of one-dimensional diffusions in order to investigate whether an allele can be (instantaneously) fixed at one of the loci.
The following proposition shows that, just as for the one-locus case, the boundary behaviour of the solution to (5.2) depends only on the values of the mutation rates µ 0 and µ 1 .
(ii) If µ 1 = 0 and µ 0 > 0 then each coordinate of (x t ) t reaches the point 0 in a finite time almost surely.
Similar statements to (ii), (iii) and (iv) hold for the point 1 on exchanging the rôles of µ 1 and µ 0 .
Proof. Let i ∈ 1 ; n . On [0, 1] n the polynomial function P i,s is bounded above by and is bounded below by for every u ∈ [0, 1]. For every i ∈ 1 ; n , pathwise uniqueness holds for the following two stochastic differential equations: Let ξ + t (i) and ξ − t (i) be the solution starting from x 0 (i) of the stochastic differential equations (5.3) and (5.4) respectively. As the i-th coordinate of the drift is bounded above by b + i and is bounded below by b − i , the comparison theorem of Ikeda & Watanabe (1977) ensures that the following inequalities hold with probability one: The nature of the points 0 and 1 as described by the Feller classification is the same for (ξ − t (i)) t and (ξ + t (i)) t and depends only on µ 1 and µ 0 . To describe their behaviours near 0, let τ ±,i z (a, b) denote the first time the process (ξ ± t (i)) t , starting from z, exits (a, b) for 0 ≤ a < z < b ≤ 1.
1. If µ 1 = µ 0 = 0 then 0 and 1 are absorbing points; (ξ ± t (i)) t reaches 0 or 1 in a finite time with probability one and IP lim 2. If µ 1 = 0 and µ 0 > 0 then 0 is the only absorbing point and (ξ ± t (i)) t reaches 0 in a finite time with probability one.
These properties are sufficient to prove the boundary behaviour claimed for (x t ) t . Since i)) t reaches 0 in a finite time then so must (x t (i)) t . Similarly, if 0 is attainable for (ξ + t (i)) t then 0 is also attainable for (x t (i)) t . In the same way, It remains to prove that (x t ) t exits from ]0, 1[ in a finite time with probability one if µ 1 = µ 0 = 0. Let ǫ > 0 be small enough that x 0 ∈ [ǫ, 1 − ǫ] n . The diffusion x t exits from the compact [ǫ, 1 − ǫ] n in a finite time with probability one. Let x ǫ be a point on the boundary of [ǫ, 1 − ǫ] n . There exists i ∈ 1 ; n such that x ǫ (i) ∈ {ǫ, 1 − ǫ}. For z ∈]0, 1[, . By the comparison theorem applied to the solutions of the stochastic differential equations (5.2), (5.3) and (5.4) starting from x ǫ , the probability that the solution of (5.2) starting at x ǫ reaches the boundary of [0, 1] n in a finite time is By the strong Markov property, the probability that (x t ) reaches the boundary in a finite time is greater ǫ)), i ∈ 1 ; n } for every ǫ > 0. Therefore (x t ) t reaches the boundary in a finite time with probability one.
6 The stationary measure of the limiting diffusion 6.1 Existence of a stationary distribution and an expression for its density As in the one-locus case, when the mutation rates are strictly positive, the diffusion has a reversible stationary distribution: Proposition 6.1. Assume that the hypothesis H4 holds and that the mutation rates µ 0 and µ 1 are strictly positive. Sets i,j = s i,j − s 1,1 for every pair of types i, j ∈ A. The diffusion with generator G n,s has a unique reversible stationary distribution which has the following density with respect to the Lebesgue measure on [0, 1] n : • C n,µ,s is chosen so that [0,1] n g n,µ,s (x 1 , . . . , x n )dx 1 · · · dx n = 1.
Remark 6.1. An expansion of the polynomial function H n,s yields: Proof of Proposition 6.1. Let G n,0 denote the generator of the limiting diffusion in the random mating case (s i,j = 0 for every i, j ∈ A). The diffusion associated with this generator is ergodic and has a reversible stationary distribution m µ,0 which is the product of Beta distributions: m µ,0 := (Beta(2µ 0 , 2µ 1 )) ⊗n . In the general case, the generator G n,s can be decomposed as Therefore, as explained in Ethier & Nagylaki (1989), we can apply a result of Fukushima & Stroock (1986) to deduce that the diffusion associated with G n,s has a unique reversible stationary distribution m µ,s given by where C is chosen so that m µ,s is a probability distribution.

Description of the density of the stationary measure
We analyse the density of the stationary measure under two supplementary assumptions: Assumption H6: The two mutation rates µ 0 and µ 1 are assumed to be equal to a strictly positive real number µ.
Assumption H7: For every L ∈ P( 1 ; n ), m L (s) depends only on |L|. We write m(ℓ) for the common value of m L (s) for L ∈ P( 1 ; n ) such that |L| = ℓ.
Assumption H7 holds if the assortment parameters satisfy the Hamming criterion or the additive criterion.
Under the hypotheses H1, H2, H3, H4, H6 and H7, the density of the invariant measure can be written as g n,µ,s (x) = C exp(h n,µ,s (x)) with The study of the invariant measure in the one-locus case already provides a precise image of the graph of g n,µ,s when the n coordinates of the diffusion are independent, that is when the assortment coefficients are chosen so that for every ℓ ∈ {0, . . . , n − 1}, m(ℓ + 1) − m(ℓ) = m(1) − m(0).
Proposition 6.2 gives conditions on the assortment parameters under which (1/2, . . . , 1/2) is the only critical point of the density, as in the random mating case. Proposition 6.3 deals with situations far from the random mating case (the proofs are postponed to §6.4).

If
Remark 6.2. The statement of Proposition 6.2 can be easily extended to a family of assortment parameters for which H7 does not hold: V n must be replaced by The following proposition describes the properties of the critical points of the density in two cases, (1) µ > 1/2 and a condition on the assortment parameters which strongly favours mating between individuals carrying similar types: [m](0) ≤ 0 and δ (1) [m](n − 2) < 0, and (2) 0 < µ < 1/2 and a condition on the assortment parameters which strongly favours mating between individuals with dissimilar types: To simplify the statement, the description is limited to the hypercube [0, 1/2] n . The description on the whole space [0, 1] n can be deduced from this since g n,µ,s (x) is invariant if we replace any coordinate x i with 1 − x i . Proposition 6.3. Assume that conditions H1, H2, H3, H4, H6 and H7 hold. Set 1. Case µ > 1/2. Assume furthermore that: (a) If V n > 0 then (1/2, . . . , 1/2) is a global maximum and is the only critical point of the density g n,µ,s .
iii. The other critical points of g n,µ,s in [0, 1/2] n are saddle points: for every ℓ ∈ 1 ; n − 1 , g n,µ,s has n ℓ saddle points of index n − ℓ in [0, 1/2] n . The saddle points of index n − ℓ have ℓ coordinates equal to 1/2 and the other coordinates have the same value denoted by ξ ℓ . iv. The relative positions of the coordinates of the critical points in [0, 1/2] n satisfy 0 < ξ n−1 < · · · < ξ 0 < 1/2. v. The value of g n,µ,s is the same at any saddle point of index n − ℓ and decreases as ℓ increases.
iii. The other critical points of g n,µ,s in [0, 1/2] n are saddle points: for every ℓ ∈ 1 ; n − 1 , g n,µ,s has n ℓ saddle points of index ℓ in [0, 1/2] n . The saddle points of index ℓ have ℓ coordinates equal to 1/2 and the other coordinates have the same value denoted by ξ ℓ .
iv. The relative positions of the coordinates of the critical points in [0, 1/2] n satisfy 0 < ξ n−1 < · · · < ξ 0 < 1/2. v. The value of g n,µ,s is the same at any saddle point of index n − ℓ and increases as ℓ increases.
3. If x i is the proportion of the population with allele 0 at the i-th locus, 2x i (1 − x i ) is the probability that two individuals sampled at random from the population carry different alleles at the ith locus. The density function of the reversible measure takes its maximum value at a point x such that for each i ∈ 1 ; n , x i (1 − x i ) = λ 0 .
Example 6.2. Let us consider a quadratic sequence of parameters s ℓ = s 0 − (bℓ + cℓ 2 ) ∀ℓ ∈ 0 ; n and let us define the assortment with this sequence by means of the Hamming criterion. If c > 0, b+c ≥ 0 and µ > 1/2 then g n,µ,s has 3 n critical points if and only if b+nc > 8µ−4.
In this case, λ 0 = n −1/2 2µ−1 4c + O(n −1 ). If h n,k denotes the value of the function h n,µ,s at a critical point of index n−k then h n,0 −h n,n ∼ n→+∞ c 8 n 2 and h n,0 −h n,1 ∼ n→+∞ n 1/2 1/2 c(2µ − 1) (see Appendix A.2 for more details). To illustrate the evolution of the 0-allelic frequency when µ > 1/2 and the assortative mating strongly favours pairing between similar types, simulations were run in a population of size N = 10 3 with the two-locus model (Fig. 5) and with the three-locus model (Fig. 6). For these simulations, every individual initially carries the allele 0 at every locus, recombination occurs independently at each locus and the assortative mating is defined by the Hamming criterion.

Graphs of the density and simulations of trajectories in the two and three locus cases
The trajectory is plotted at intervals of size N between the iterations N 2 and 33N 2 . To help to visualize the evolution, the colour of the plot changes every 1 2 N 2 iterations. The form of the density of the stationary measure here is highly reminiscent of that of the fitness landscapes studied in the adaptive evolution literature in modelling additive traits under frequency dependent intraspecific competition, see e.g. Schneider (2007) and references therein. In the deterministic setting the existence of multiple 'long term equilibria' renders the behaviour of the system very sensitive to assumptions about the initial conditions. In our setting, the presence of genetic drift is sufficient for the population to (eventually) explore the neighbourhoods of all the maxima, irrespective of its starting point. The time spent by the population in the neighbourhood of a maximum depends on the assortment parameters ( Fig. 6a and 6b).  6.4 Proofs of Propositions 6.2 and 6.3 Proof of Proposition 6.2 Let us introduce some notation in order to shorten the expres- and h(x) =h(ρ(x)) for x = (x 1 , . . . , x n ) ∈]0, 1[ n . With this notation, g n,µ,s (x) = C n,µ,s exp(h(x)).
1. For every x ∈]0, 1[ n and i ∈ 1 ; n , First, the point u n = (1/2, . . . , 1/2) is a critical point of h n,µ,s and the Hessian matrix at this point is the diagonal matrix −2∆I n where This proves the first two assertions of the proposition. In this case, ∆ = 0 and the stationary density has only one critical point at (1/2, . . . , 1/2) which is a maximum. If {s i,j , (i, j) ∈ A 2 } is a family of assortment parameters such that ∂ ih (x) is nonnegative for every x ∈]0, 1/4] n and the density g n,s,µ has a unique critical point at (1/2, . . . , 1/2) which is a maximum, then the same is true for any family of assortment parameters for every ℓ ∈ 0 ; n − 1 .
Proof of Proposition 6.3 We retain the notation introduced in the proof of Proposition 6.2.
In the proof we shall use (several times) the following identity for elementary symmetric polynomial functions: Lemma 6.1. Let n be an integer greater than 1 and let k ∈ 0 ; n − 2 . For every x ∈ IR n , set . . , x n ) for i ∈ 1 ; n and Then, We shall also use the following alternative expression for symmetric polynomial functions that are similar to the polynomial term in h: Lemma 6.2. Let n ∈ IN * and let a 0 , . . . , a n be real numbers. Then for every x ∈ IR n , (1 − 2x j ).
1. Let us assume that x = (x 1 , . . . , x n ) is a critical point of g n,µ,s different from u n . Let ℓ denote the number of coordinates equal to 1/2 (ℓ ∈ 0 ; n − 1 ). Every coordinate x i different from 1/2 has to satisfy: ∂ i h n,µ,s (ρ(x)) = 0, that is In particular, it follows from Lemma 6.1 that if x i and x j are two coordinates of the critical point x not equal to 1/2 then By Lemma 6.2, where Q ℓ denotes a polynomial function which is positive on Therefore, (E ′ ℓ ) coincides with (E ℓ ) of Remark 6.3. The derivative of φ ℓ is equal to: For every pair of disjoint subsets I and J of 1 ; n , let us introduce the following point: We have shown that if V n and ν have opposite signs, then every point u I,J is a critical point and any critical point is one of these points u I,J .
So that we may use our conclusions above, from now on, we assume that the hypotheses stated in point 1 of the proposition are satisfied. However, the computations that follow do not depend on these hypotheses, and so our proof is easily modified to the setting of point 2.
2. Let us study the Hessian matrix of h n,µ,s at a critical point u I,J such that |I| ≤ n − 1.
For that, set ℓ = |I|, ℓ + = |J| and ℓ − = n − ℓ − ℓ + and let us introduce the following notations: The Hessian matrix of h n,µ,s at u I,J is permutation-similar to the following block matrix: • B ℓ,k denotes the following k-by-k matrix : • C ℓ denotes the ℓ + -by-ℓ − matrix all the elements of which are equal to −c ℓ .
By assumption on µ, b ℓ < 0. To complete the proof of assertions (i) and (ii) of 1-(b), we shall prove that a ℓ < 0 and that b ℓ < c ℓ < 0. That will imply that the submatrix B ℓ,ℓ + C ℓ C ℓ B ℓ,ℓ − is negative definite (for more details, see Lemma A.2) hence that the Hessian matrix of h n,µ,s at a point u I,J has |I| positive eigenvalues and n − |I| negative eigenvalues.

Proof of convergence to the diffusion
In this section, we prove convergence to the diffusion approximation in the n-locus case (Theorem 4.1). We also establish the two simple expressions for the drift presented in §4.
Then the operator is closable in C([0, 1] n ) and its closure is the generator of a strongly continuous semigroup of contractions.
To prove the convergence result, we use the following theorem, due to Ethier and Nagylaki, on diffusion approximations for Markov chains with two time scales.
Theorem 7.2 (Ethier & Nagylaki 1980, Theorem 3.3). For N ∈ IN * , let {Z N k , k ∈ IN} be a homogeneous Markov chain in a metric space E N with Feller transition function. Let F 1 and F 2 be compact convex subsets of IR n and IR m respectively, having non-empty interiors. Assume further that 0 ∈ • F 2 . Let Φ N : E N → F 1 and Ψ N : E N → F 2 be continuous functions. Define Let (ǫ N ) N and (δ N ) N be two positive sequences such that δ N → 0 and ǫ N /δ N → 0.
Assume that there exist continuous functions a : F 1 × IR m → IR n ⊗ IR n , b : F 1 × IR m → IR n and c : F 1 × IR m → IR m such that for i, j ∈ 1 ; n and ℓ ∈ 1 ; m the following properties (a)-(e) hold as N → +∞ uniformly in z ∈ E N where x = Φ N (z) and y = Ψ N (z): Assume further that (f ) c is of class C 2 , c(x, 0) = 0 for all x ∈ IR m and the solution of the differential equation d dt u(t, x, y) = c(x, u(t, x, y)), u(0, x, y) = y.
(g) The closure of the following operator generates a strongly continuous semigroup on C(F 1 ) corresponding to a diffusion process Then the following conclusions in which the symbol ⇒ denotes convergence in distribution, hold: Remark 7.1. We have only stated the part of Ethier and Nagylaki's theorem that we need.
The full statement also gives a convergence result when the sequence (δ N ) N converges to a positive real number.
To apply this theorem, we consider the two sequences ǫ N = N −2 and δ N = N −1 , we set where u i = ℓ,ℓ i =0 z(ℓ) for i ∈ 1 ; n and u I = i∈I u i − ℓ,ℓ |I ≡0 z(ℓ) for each I ⊂ 1 ; n having at least two elements.
First (in §7.1), we shall check that X 1 ) satisfy the conditions (a)-(f) of Ethier and Nagylaki's theorem with the following expressions for the functions a i,j (x, 0) and b i (x, 0): and, for two subsets I and J of 1 ; n , s I,J denotes the assortment parameter s i,j for the types i = (0 I , 1Ī ) and j = (0 J , 1J ).
In §7.2 we shall show that P i,s has the following two equivalent expressions:

Verification of the conditions (a)-(f) of Ethier and Nagylaki's theorem
As the proportion of individuals of a given type i can only change by ±1/N in one step: • If r ∈ IN * and i ∈ A, then • if r ≥ 3 and i (1) ,. . . ,i (r) ∈ A so that at least three of them are distinct, then using assumption H2 yields the following formula: Lemma 7.1. For every i ∈ A, Proof. By assumption H2, for two different types i, j ∈ A To prove Lemma 7.1, it suffices to use these expansions in and to simplify.
Let u ∈ 1 ; n . To establish an expression for the drift of X (N ) (u), we must compute i∈A,iu=0 B i (z) and i∈A,iu=0 B (1) i (z). Direct computations yield: Lemma 7.2. For every u ∈ 1 ; n and z ∈ E N , i∈A, iu=0 Proof. For ǫ ∈ {0, 1} and i ∈ A, let σ (ǫ) u (i) denote the type i modified by setting the allele ǫ at the locus u. We shall use the following formula several times: withr(u) = I⊂ 1;n \{u} r I = 1 2 by assumption H1. First, formula (7.8) with ǫ = 0 provides i∈A, iu=0 (1,j) i (z) denote the j-th line of the expression of B (1) i (z) for j ∈ {1, 2, 3}. As i∈A, iu=0,ix=a q((j, k); σ ǫ x (i)) does not depend on the value of a if u = x: Applying (7.8) again, we obtain: Due to the symmetry of the parameters: s i,j = s j,i for i, j ∈ A, we have: Finally, computations using (7.8) yet again yield: To obtain condition (a), it remains to express G u (z) in the new coordinates. The following lemma describes the inverse of the change of coordinates (Φ N , Ψ N ): Lemma 7.3. For z ∈ E N and L ⊂ 1 ; n , set x(L) = i, i |L ≡0 z(i) with the convention x(∅) = 1 and y(L) = ℓ∈L x(ℓ) − x(L) if |L| ≥ 2. Then for every J ⊂ 1 ; n , (−1) |I|−|J| y(I). (7.9) Proof. First, by induction on n − |J|, we show that (7.10) Since z(0) = x( 1 ; n ), the equality (7.10) holds for J = 1 ; n .
Let m ∈ 1 ; n . Assume that the formula (7.10) holds for every subset J of 1 ; n such that |J| ≥ m. Let K be a subset of 1 ; n with m − 1 elements.
We apply the formula (7.10) to every term in the sum and we invert the double sum we have obtained: The sum between parentheses is equal to Thus the formula (7.10) is also satisfied for the subset K which completes the induction.
To complete the proof, we replace x(I) in (7.10) with i∈I x(i) − y(I) for every subset I having at least two elements and use the following equality: (1 − x(i)).
With this notation, for every J ⊂ Λ u , where R J (y) and R J∪{u} (y) denote polynomial functions that vanish at y ≡ 0. Therefore, where R u (x, y) is a polynomial function in the variables x(1), . . . , x(n) and y(I) for I ⊂ 1 ; n such that |I| ≥ 2, that vanishes in the equilibrium manifold: R u (x, 0) = 0.
The expression for G u (z) can be simplified by using the two assumptions H4 on the assortment parameters, that is s J,H = s H,J for every J, H ⊂ 1 ; n and s J∪{u},H∪{u} = s J,H for every u ∈ 1 ; n and J, H ⊂ Λ u : In summary, we have established the following expansion of the drift of X (N ) : Lemma 7.4. Assume that hypotheses H1, H2, H3 and H4 hold. For every i ∈ 1 ; n , and R i (x, y) is a polynomial function in the variables x(1), . . . , x(n) and y(I) for I ∈ P( 1 ; n ) with at least two elements such that R i (x, 0) = 0.

Condition (b).
Computations similar to those used to obtain (7.6) lead to the following expansion of the second moments of X showing that condition (b) holds: Proof. Let i, j ∈ 1 ; n and z ∈ E N . By definition of X (N ) , Using formulae (7.3) and (7.4) and assumption H2, we obtain i,j = x({i, j}) and it follows from assumption H1 (r I = rĪ for every I ⊂ 1 ; n ) that Therefore, for every i, j ∈ 1 ; n , I⊂ 1;n \{i,j} r I = 1 2 .
Condition (d). Let I be a subset of 1 ; n with at least two elements. To compute the drift of Y (N ) (I), we use the following lemma and formulae (7.3), (7.4) and (7.5) describing the Lemma 7.6. Let J be a finite set. Consider two families of reals {a j , j ∈ J} and {b j , j ∈ J}.
The following identity holds: Computations yield: uniformly on z ∈ E N . As we have shown that j∈A, j i =0 B (0) j (z) = 0 for every i ∈ 1 ; n (equation (7.6)), Direct computations provide the following expression of the sum on the right-hand side of (7.14) using the variables x(L) = j∈A, j |L ≡0 x(j) for L ∈ P( 1 ; n ): Condition (f ). The following lemma shows that the condition (f) holds under the assumption H3: Lemma 7.8. For two distinct loci k, ℓ, let r k,ℓ denote the probability that the offspring does not inherit the genes at the loci k and ℓ from the same parent, r k,ℓ = I⊂ 1;n , k∈I and ℓ ∈I (r I + rĪ ), and set r(n) = min(r k,h k, h ∈ 1 ; n and h = k).
If r(n) > 0 then the following system of differential equations (S n,I ) dv n,I dt (t, x, y) = c n,I (x, v n,I (t, x, y)) v n,I (0, x, y) = y(I) ∀I ⊂ 1 ; n s. t. |I| ≥ 2 has a unique solution v n = {v n,I , I ⊂ 1 ; n and |I| ≥ 2} which is of the form: Remark 7.2. For every subset I ⊂ 1 ; n with two elements say k and ℓ, dv n,I dt (t, x, y) = −r k,ℓ v n,I (t, x, y).
Therefore if r(n) = 0 then there exists a subset I of 1 ; n with two elements such that v n,I (t, x, y) = y(I). Thus the assumption r(n) > 0 is a necessary condition for the solution of (S n,I ) to converge to 0 as t tends to +∞ for any initial values.
Proof. Let n ≥ 2 and let I ⊂ 1 ; n be such that |I| ≥ 2. As c n,I (x, y) depends only on the coordinates x(ℓ) for ℓ ∈ I and y(L) for L ⊂ I such that |L| ≥ 2, we shall prove by induction on the number of elements of I that for any J ⊂ I, (S n,J ) has a unique solution of the form v n,J (t, x, y) = exp(−r(n)t)f n,J (t, x, y), where f n,J is a continuous and bounded function on IR ×[0, 1] n × [−1, 1] 2 n −n−1 such that the value of f n,J (t, x, y) depends on x and y only through the values of the coordinates x(j) for j ∈ J and y(L) for L ⊂ J such that |L| ≥ 2.
• If I has two elements say k and ℓ, then (S n,I ) is the following differential equation: dv n,I dt (t, x, y) = −r k,ℓ v n,I (t, x, y) v n,I (0, x, y) = y(I) It has a unique solution v n,I (t, x, y) = y(I)e −r(2)t f n,I (t, x, y) where f n,I (t, x, y) = e −(r k,ℓ −r(2))t y(I).
• Let 2 ≤ m < n. Assume that the inductive hypothesis holds for any subsets J with m elements. Let I be a subset of 1 ; n with m + 1 elements. Then r L e −tr(n) f n,I∩L (t, x, y)f n,I∩L (t, x, y) Asr I is the probability that the offspring does not inherit all the genes at loci i ∈ I from the same parent,r I ≥ r(n). Therefore the differential equation (S n,I ) has a unique solution: v n,I (t, x, y) = y(I)e −r I t + e −r I t t 0 g(s, x, y)e (r I −r(n))s ds.
By our assumptions on the functions f n,J for J I, g is a bounded continuous function on IR + ×[0, 1] n × [−1, 1] 2 n −n−1 such that the value of g(t, x, y) depends on x and y only through the coordinates x(i) for i ∈ I and y(L) for L ⊂ I such that |L| ≥ 2. Therefore, the function f n,I (t, x, y) = e r(n)t v n,I (t, x, y) has the asserted properties.

Expressions for the drift
We have shown that the i-th coordinate of the drift of the limiting diffusion is and, for two subsets I and J of 1 ; n , s I,J denotes the assortment parameter s i,j for the types i = (0 I , 1Ī ) and j = (0 J , 1J ). The following lemma states that P i,s (x) is actually a polynomial function in the variables x(i)(1 − x(i)) for i ∈ 1 ; n \ {u}: Proof. Let P Λ (β) denote the polynomial function on the right-hand side. The proof is by induction on |Λ|. First, P ∅ (β)(x) = β ∅,∅ = C ∅ (β).
Let n ∈ IN. Assume that the equality (7.17) holds for every subset Λ of IN with at most n elements and every family of reals β satisfying the assumptions of the lemma.
Let Λ be a subset of IN with n+1 elements, let j be an element of Λ and let η = {η I,J , I, J ⊂ Λ} be a family of reals such that η I,J = η I\J,J\I for every I, J ⊂ Λ. We split P Λ (η) into a sum over the subsets of Λ containing j and a sum over the subsets of Λ \ {j} to obtain the following expression: x(j) 2 η K∪{j},L∪{j} + (1 − x(j)) 2 η K,L + x(j)(1 − x(j))(η K∪{j},L + η K,L∪{j} ) .
The following factorised form of the polynomial function P i,s can be derived from a general identity stated in Lemma A.1:

A.1 Combinatorial formulae for difference operators
This section collects some combinatorial formulae used to study the limiting diffusion. Let E be a finite set and t be a real. For a function f defined on P(E), we set This shows the equality between the expanded form (4.3) and factorised form (4.1) of the polynomial term P i,s (x) appearing in the drift of the limiting diffusion.

A.3 Property of a symmetric matrix
The following lemma is used to determine the nature of the critical points of the density of the invariant measure (Proposition 6.3).
Lemma A.2. For a real a and two integers k and n so that n ≥ 1 and 0 ≤ k ≤ n, let M n,k (a) denote the following symmetric matrix: • B k 1 ,k 2 denotes the k 1 -by-k 2 matrix all the elements of which are equal to −a.
If 0 ≤ a < 1 then M n,k (a) is positive definite.