A M Etheridge
email: etheridg@stats.ox.ac.uk.

S Lemaire
email: sophie.lemaire@math.u-psud.fr

Diffusion approximation of a multilocus model with assortative mating

Keywords: AMS 2000 Subject Classification: 60J20, 92D25, 60J70 Moran model, population genetics, multilocus models, assortative mating, diffusion approximation

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). [START_REF] Ethier | Diffusion approximations of the two-locus Wright-Fisher model[END_REF] 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) 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 [START_REF] Baake | Single-crossover dynamics: finite versus infinite populations[END_REF]. A review of several simulation studies can be found in the introduction of [START_REF] Devaux | Incipient allochronic speciation due to non-selective assortative mating by flowering time, mutation and genetic drift[END_REF]. Among these, the 'species formation model', introduced by [START_REF] Higgs | Genetic distance and species formation in evolving populations[END_REF] 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 1s 0 = s 2s 1 = . . . = s ns 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 1s 0 ); that is they solve the following stochastic differential equation:

dx t = x t (1 -x t )dW t + µ 1 (1 -x t ) -µ 0 x t + (s 1 -s 0 )(1/2 -x t )x t (1 -x t ) dt,
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:

dx t (i) = x t (i)(1 -x t (i))dW t (i)
+ µ 1 (1x t (i))µ 0 x t (i) + (1/2x t (i))x t (i)(1x t (i))P i,s (x t ) dt, (1.1)

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)(1x(j)), j ∈ {1, . . . , n} \ {i} whose coefficients depend only on the parameters s 1s 0 , . . . , s ns n-1 . More precisely, P (x (i) )

is an increasing function of each parameter s 1s 0 , . . . , s ns 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 1s 0 , . . . , s ns 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 ns n-1 ≤ . . . ≤ s 1s 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

m L∪{i} (s) -m L (s) for L ⊂ {1, . . . , n} replace s 1 -s 0 ,. . . , s n -s n-1 .
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 [START_REF] Ethier | Diffusion approximations of the two-locus Wright-Fisher model[END_REF].

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) so that the composition of the population is described by the set Z

k := (k 1 , ..., k n ) ∈ {0, 1} n . Let A = {0, 1}
(N ) t = {Z (N ) t (k), k ∈ A}.
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 (N ) i,j Z (N ) t (j) k∈A 1 + s (N ) i,k Z (N ) t (k) ,
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 [START_REF] Geiringer | On the probability theory of linkage in Mendelian heredity[END_REF] 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 where L denotes the complementary set of loci, 1 ; n \ L.

With this notation, the probability that, before mutation, the offspring of a pair of individ-

uals of types (i, j) is of type k is q((i, j); k) = L⊂ 1;n r L 1I {k=(i |L ,j | L)} .
Let us express some classical examples of recombination distributions in this notation:

Examples 2.1.

1. r ∅ = r 1;n = 1 2 (no recombination, also called absolute linkage) 2. r I = 2 -n for each I ∈ P( 1 ; n ) (free recombination) 3. r 1;x = r x+1;n = r 2(n-1) for 1 ≤ x ≤ n -1 and r ∅ = r 1;n = 1/2(1 -r)
where r denotes an element of ]0, 1] (at most one exchange between the sequence of loci which occurs with equal probability at each position).

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:

µ (N ) (k, ℓ) := n i=1 (µ (N ) k i ) |ℓ i -k i | (1 -µ (N ) k i ) 1-|ℓ i -k i | .

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

f N (z, i, j) := k,ℓ∈A z(i)z(k)w (N ) (z, i, k)q((i, k); ℓ)µ (N ) (ℓ, j)
where

w (N ) (z, i, k) = 1 + s (N ) i,k ℓ∈A (1 + s (N ) i,ℓ )z(ℓ)
.

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 Assumption H2: µ

(N ) ǫ = µǫ N for ǫ ∈ {0, 1} and s (N ) i,j = s i,j
N for i, j ∈ A. Just as in the two-locus case studied by [START_REF] Ethier | Diffusion approximations of the two-locus Wright-Fisher model[END_REF], 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 [START_REF] Geiringer | On the probability theory of linkage in Mendelian heredity[END_REF][START_REF] Nagylaki | The evolution of multilocus systems under the weak selection[END_REF] 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 .

In the two-locus case, this assumption leads to a model with five assortment parameters:

• one parameter, s 00,00 = s 11,11 = s 10,10 = s 01,01 , for pairs of individuals having the same genotype,

• one parameter s 00,10 = s 10,00 = s 11,01 = s 01,11 for pairs of individuals whose genotypes only differ on the first locus,

• one parameter, s 00,01 = s 01,00 = s 11,10 = s 10,11 , for pairs of individuals whose genotypes only differ on the second locus,

• two parameters s 01,10 = s 10,01 and s 00,11 = s 11,00 for pairs of individuals whose genotypes differ on the two loci.

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 [START_REF] Higgs | Genetic distance and species formation in evolving populations[END_REF] 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. For example, [START_REF] Devaux | Incipient allochronic speciation due to non-selective assortative mating by flowering time, mutation and genetic drift[END_REF] use it to investigate speciation in flowering plants due to assortative mating determined by flowering time. Flowers can only be pollinated by other flowers that are open at the same time. Modelling flowering time as an additive trait, they observe an effect that is qualitatively similar to that observed in the simulations of [START_REF] Higgs | Genetic distance and species formation in evolving populations[END_REF] 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

(1) i |G 1 ,j |G 1 + s (2) i |G 2 ,j |G 2
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:

IE z [Z (N ) 1 (0) -z] = 1 N 2 (1 -z)µ 1 -zµ 0 + 1 2 z(1 -z)((s 1,0 -s 1,1 )(1 -z) -(s 0,1 -s 0,0 )z) + O(1/N 3 ) IE z [(Z (N ) 1 (0) -z) 2 ] = 1 N 2 z(1 -z) + O(1/N 3 ) IE z [(Z (N ) 1 (0) -z) 4 ] = O(1/N 4 ) uniformly in z.
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:

G 1,s = 1 2 x(1 -x) d 2 dx 2 + (1 -x)µ 1 -xµ 0 + 1/2x(1 -x)((s 1,0 -s 1,1 )(1 -x) -(s 0,1 -s 0,0 )x) d dx . (3.1) More precisely, if Z (N ) 0 converges in distribution in [0, 1] as N tends to +∞, then (Z (N )
[N 2 t] ) t≥0 converges in distribution in the Skorohod space of càdlàg functions D [0,1] ([0, +∞)) to a diffusion with generator G 1,s (see, for example, Ethier & Kurtz 1986, Chapter 10).

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

G 1,s = 1 2 x(1 -x) d 2 dx 2 + (1 -x)µ 1 -xµ 0 + (s 1 -s 0 )(1/2 -x)x(1 -x) d dx .
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 homozygosity6 (when s 0s 1 > 0) or in favour of heterozygosity (when s 0s 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:

g µ,s (x) = C µ,s x 2µ 1 -1 (1 -x) 2µ 0 -1 exp -1/2((s 1,0 -s 1,1 )(1 -x) 2 + (s 0,1 -s 0,0 )x 2 )
where the constant C µ,s is chosen so that 1 0 g µ,s (x)dx = 1. This is plotted, for various parameter values, in Fig. 1 under the assumptions µ 1 = µ 0 = µ, s 0,0 = s 1,1 = s 0 and s 0,1 = s 1,0 = s 1 . (i) if µ 1 = 0 then 0 is an absorbing state and the diffusion exits from ]0, 1[ in a finite time almost surely;

(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

the drift IE[Z (N ) t+1 (i) -Z (N ) t (i) | Z (N ) t
= 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-

type distribution z is in W n , that is if z(i) = z 1 (i 1 ) • • • z n (i n ) ∀i = (i 1 , . . . , i n ) ∈ A where z j (x) = k∈{0,1} n-1 z(k 1 , . . . , k j-1 , x, k j+1 , .
. . , 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 nn -1 processes that measure the deviation from the linkage equilibrium.

For a nonempty subset L of {1, . . . , n},

• let X (N ) t (L) = j∈A, j |L ≡0 Z (N ) t
(j) denote the proportion of individuals having the allele 0 on all loci in L at time t;

• let Y (N ) t (L) = i∈L X (N ) t ({i}) -X (N ) t
(L) for |L| ≥ 2 describe the linkage disequilibrium between the loci in L at time t. (This is just one of many ways to measure the linkage disequilibrium, see for example [START_REF] Bürger | The mathematical theory of selection, recombination, and mutation[END_REF], Chapter V.4.2, for other measures.)

The vector of 0-allelic frequencies at time t is X

(N ) t := X (N ) t ({1}), ..., X (N ) t ({n}) . The process Y (N ) defined by Y (N ) t := Y (N ) t (L), L ⊂ 1 ; n such that |L| ≥ 2 for t ≥ 0 vanishes on the Wright manifold.
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:

F L = (i, j) ∈ A 2 : i u = 1 -j u ∀u ∈ L and i u = j u ∀u ∈ L .
Let m L (s) denote the mean value of the assortment parameters for all pairs in this set F L :

m L (s) = 2 -n (i,j)∈F L s i,j .
Examples 4.1.

1. In the two-locus case, ).

m ∅ (s) =
In each of these expressions the four coefficients are equal by assumption H4.

m {1,2} (s) = 1 4 (s 00,11 + s 11,00 + s 01,10 + s 10,01 ).

In this expression the first (resp. last) two coefficients are equal by H4.

2. With the Hamming criterion, m L (s) = s |L| for every L ⊂ 1 ; n , where |L| denotes the number of loci in L.

With the additive criterion

, m ∅ (s) = s 0 , m {ℓ} (s) = s 1 ∀ℓ ∈ 1 ; n and more generally m L (s) = 2 -|L| |L| k=0 |L| k s |2k-|L|| for every L ⊂ 1 ; n .

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 [START_REF] Ethier | Diffusion approximations of Markov chains with two time scales and applications to population genetics[END_REF], is postponed to §7.

Theorem 4.1. Assume that hypotheses H1, H2, H3 and H4 hold.

(a) For i ∈ 1 ;

n let P i,s (x) denote a polynomial function in the n -1 variables x k (1 -x k ) for k ∈ 1 ; n \ {i}.
Then, the operator

G n,s = 1 2 n i=1 x i (1 -x i ) ∂ 2 ∂ x i ∂ x i + n i=1 (1 -x i )µ 1 -x i µ 0 + (1/2 -x i )x i (1 -x i )P i,s (x) ∂ ∂ x i (4.1) with domain D(G n,s ) = C 2 ([0, 1] n ) is closable in C([0, 1] n
) and its closure is the generator of a strongly continuous semigroup of contractions.

(b) If X (N ) 0 converges in distribution in [0, 1] n , then (X (N ) [N 2 t] ) t converges in distribution in the Skorohod space of càdlàg functions D [0,1] n ([0, ∞))
to a diffusion process X with generator G n,s where the polynomial function P i,s (x) has the following expression:

P i,s (x) = A∈P( 1;n \{i}) m A∪{i} (s) -m A (s) k∈A (2x k (1 -x k )) ℓ∈ 1;n \{A∪{i}} (1 -2x ℓ (1 -x ℓ )). (4.2) (c) For every positive sequence (t N ) N that converges to +∞, Y (N ) [N t N ] converges in distribution to 0. 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 P i,s (x) of the drift

An expansion of the polynomial function P i,s (x) in terms of the variables x k (1x k ), k ∈ 1 ; n \ {i} yields the following expression:

P i,s (x) = L∈P( 1;n \{i}) α i,L (s) ℓ∈L x ℓ (1 -x ℓ ) (4.3) with α i,L (s) = 2 |L| A⊂L (-1) |L|-|A| (m A∪{i} (s) -m A (s)). (4.4)
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 element i of E, we denote by δ i [f ] the function on P(E) defined by

δ i [f ](A) = f (A ∪ {i}) -f (A), ∀A ∈ P(E).
Since δ i • δ j = δ j • δ i for every i, j ∈ E, we can, more generally, introduce a difference operator δ B for each subset B ∈ P(E) by setting δ ∅ = Id, and

δ B = δ b 1 • • • • • δ br if B = {b 1 , . . . , b r }. A
proof by induction on |B| provides the following formula for δ B :

δ B [f ](A) = J⊂B (-1) |B|-|J| f (A ∪ J) ∀A ⊂ E. (4.5)
Let m(s) denote the function A → m A (s) defined on the subsets of 1 ; n . In this notation, for

every A ⊂ 1 ; n \ {i}, m A∪{i} (s) -m A (s) = δ i [m(s)](A) and α i,A (s) = 2 |A| δ A∪{i} [m(s)](∅). (4.6)
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:

P i,s (x) = n-1 ℓ=0 αℓ L⊂ 1;n \{i}, |L|=ℓ ℓ∈L x ℓ (1 -x ℓ ). (4.7) where αk (s) = 2 k k ℓ=0 (-1) ℓ k ℓ (s k-ℓ+1 -s k-ℓ ).
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+1s i for every i ∈ 0 ; n -1 . The forward difference operators of higher orders are defined iteratively:

δ (k+1) [s] = δ (k) •δ (1) [s] for k ∈ IN * . With this notation, αk (s) = 2 k δ (k+1) [s](0) for k ∈ 0 ; n -1 .

Comparison with the two-locus Wright-Fisher diffusion

Ethier & [START_REF] Ethier | Diffusion approximations of the two-locus Wright-Fisher model[END_REF] 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:

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 * i,j = w N,i,j

P i P j k,ℓ w N,k,ℓ P k P ℓ if P k
denotes the frequency of gametes k in the population ∀k ∈ 1 ; r 1 × 1 ; r 2 );

2. the recombination fraction c N ;

3. the probability (2N ) -1 ν (i) 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

L = 1 2 p 1 (1 -p 1 )∂ 2 p 1 ,p 1 + 1 2 q 1 (1 -q 1 )∂ 2 q 1 ,q 1 + b 1 (p 1 , q 1 )∂ p 1 + b 2 (p 1 , q 1 )∂ q 1 with b 1 (p 1 , q 1 ) =ν (1) 2,1 (1 -p 1 ) -ν (1) 1,2 p 1 -p 1 (1 -p 1 )(1 -2p 1 ) (σ 12,21 + σ 11,22 )q 1 (1 -q 1 ) + σ 11,21 q 2 1 + σ 12,22 (1 -q 1 ) 2 -2p 1 (1 -p 1 )q 1 (1 -q 1 ) σ 11,12 p 1 -σ 21,22 (1 -p 1 ) . b 2 (p 1 , q 1 ) =ν (2) 2,1 (1 -q 1 ) -ν (2) 1,2 q 1 -q 1 (1 -q 1 )(1 -2q 1 ) (σ 12,21 + σ 11,22 )p 1 (1 -p 1 ) + σ 11,12 p 2 1 + σ 21,22 (1 -p 1 ) 2 -2q 1 (1 -q 1 )p 1 (1 -p 1 ) σ 11,21 q 1 -σ 12,22 (1 -q 1 ) .
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

Lemma 5.1. Any generator on C 2 ([0, 1] n ) of the form 1 2 n i=1 x i (1 -x i ) ∂ 2 ∂ x i ∂ x i + n i=1 (1 -x i )µ 1 -x i µ 0 + (1/2 -x i )x i (1 -x i ) L∈P( 1;n \{i}) α L∪{i} k∈L x k (1 -x k ) ∂ ∂ x i ,
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.

• s i,j = B⊂L, |B|≥1 2 -|B|+1 α B for every (i, j) ∈ F L and for every nonempty subset L of 1 ; n .

Let us check that this family satisfies 2

|L|-1 δ L [m(s)](∅) = α L for every nonempty subset L of 1 ; n . First, m L (s) = B⊂L, |B|≥1 2 -|B|+1 α B . For every i ∈ 1 ; n and L ⊂ 1 ; n \ {i} δ L∪{i} [m(s)](∅) = A⊂L (-1) |L|-|A| (m A∪{i} (s) -m A (s)) = A⊂L (-1) |L|-|A| B⊂A 2 -|B| α B∪{i} .
We invert the double sum and use the formula

A⊂L, s. t. B⊂A (-1) |L|-|A| = 1I {L=B} to obtain: δ L∪{i} [m(s)](∅) = B⊂L 2 -|B| α B∪{i} 1I {L=B} = 2 -|L| α L∪{i} .
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:

1 2 n i=1 x i (1 -x i ) ∂ 2 ∂ x i ∂ x i + n i=1 (1 -x i )µ 1 -x i µ 0 + (1/2 -x i )x i (1 -x i ) n-1 ℓ=0 α ℓ L⊂ 1;n \{i} s.t.|L|=ℓ k∈L x k (1 -x k ) ∂ ∂ x i .
To see this, given any sequence α 0 , . . . , α n-1 of n reals, we have to find n + 1 real numbers ℓ+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 0 , . . . , s n such that α ℓ = 2 ℓ δ (
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

s i,j = s (1) i |G 1 ,j |G 1 + s (2) i |G 2 ,j |G 2 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 nk coordinates and that the generator of the diffusion limit is:

G n,s = G k,s 1 ⊗ G n-k,s 2 .
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.

1. The n coordinates of the diffusion associated with the generator G n,s are independent if and only if the following condition holds:

(H5) for every i ∈ 1 ; n , m L∪{i} (s)m L (s) does not depend on the choice of the subset

L of 1 ; n \ {i}.
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

1 2 x(1 -x) d 2 dx 2 + (1 -x)µ 1 -xµ 0 + (s u i ,1 -s 1,1 )(1/2 -x)x(1 -x) d dx .
3. In particular, (a) with the Hamming criterion, G n,s is the generator of n independent one-dimensional diffusions if and only if the value of s ℓ+1s ℓ does not depend on ℓ;

(b) with the additive criterion, G n,s is the generator of n independent one-dimensional diffusions if and only if there exists a constant c such that s ℓ+1s ℓ = c(2ℓ + 1) for every ℓ ∈ 0 ; n -1 .

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 ),

δ i [m(s)](A) = B⊂A δ B∪{i} [m(s)](∅) = δ i [m(s)](∅).
Therefore, condition H5 is satisfied.

2. With the Hamming criterion, m A (s) = s |A| and condition H5 is equivalent to

s k+1 -s k = s 1 -s 0 for every k ∈ 1 ; n -1 .
3. With the additive criterion, for a subset L with ℓ elements m L (s) = 2 -ℓ ℓ j=0 ℓ j s |2j-ℓ| .

After some computation, we obtain for i ∈ L,

m L∪{i} (s) -m L (s) = 2 -ℓ ℓ j=0 ℓ j (s |2j-ℓ+1| -s |2j-ℓ| ) =            2 -ℓ ℓ+1 2 j=1 ℓ ℓ+1 2 -j δ (2) [s](2j -2) if ℓ is odd, 2 -ℓ ℓ 2 j=1 ℓ ℓ 2 -j δ (2) [s](2j -1) + ℓ ℓ 2 δ (1) [s](0)
if ℓ is even.

(5.1)

It follows from (5.1) that for every c ∈ IR, the system defined by m L∪{i} (s)m L (s) = c for every i ∈ 1 ; n and L ⊂ 1 ; n \ {i} has a unique solution which is δ (1) [s](k) = c(2k + 1) for every k ∈ 0 ; n -1 .

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.

Consider the stochastic differential equations associated with the generator G n,s :

dx t (i) = x t (i)(1 -x t (i))dW t (i) + b i (x t )dt ∀i ∈ 1 ; n , (5.2)
where (W t (1)) t≥0 ,. . . , (W t (n)) t≥0 denote n independent standard Brownian motions, and

b i (x) = µ 1 (1 -x(i)) -µ 0 x(i) + (1/2 -x(i))x(i)(1 -x(i))P i,s (x) for i ∈ 1 ; n .
Theorem 1 of [START_REF] Yamada | On the uniqueness of solutions of stochastic differential equations[END_REF] ensures pathwise uniqueness for the stochastic differential equation (5.2), since the drift is Lipschitz and the diffusion matrix is a diagonal matrix of the form

σ(x) = diag(σ 1 (x(1)), . . . , σ n (x(n))),
where the functions σ i are 1/2-Hölder continuous functions.

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 .

Proposition 5.2. Let (x t ) t≥0 denote a solution of the stochastic differential equation (5.2)

starting from a point x 0 ∈]0, 1[ n .
(i) If µ 1 = µ 0 = 0 then the diffusion process (x t ) t exits from ]0, 1[ n in a finite time almost surely.

(ii) If µ 1 = 0 and µ 0 > 0 then each coordinate of (x t ) t reaches the point 0 in a finite time almost surely.

(iii) If 0 < µ 1 < 1/2 then 0 is attainable for each coordinate of the diffusion process:

IP[∃t > 0, x t (i) = 0] > 0 ∀i ∈ 1 ; n .
(iv) If µ 1 ≥ 1/2 then 0 is inaccessible for each coordinate of the diffusion process:

IP[∃t > 0, x t (i) = 0] = 0 and IP[ lim t→+∞ x t (i) = 0] = 0 for every i ∈ 1 ; n .
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

M + i = A⊂ 1;n \{i} 2 -|A| max m A∪{i} (s) -m A (s), 0
and is bounded below by

M - i = - A⊂ 1;n \{i} 2 -|A| max -(m A∪{i} (s) -m A (s)), 0 . Let b + i and b - i denote the functions defined on [0, 1] by b + i (u) = µ 1 (1 -u) -µ 0 u + (1/2 -u)u(1 -u)(M + i 1I {u<1/2} +M - i 1I {u>1/2} ), b - i (u) = µ 1 (1 -u) -µ 0 u + (1/2 -u)u(1 -u)(M + i 1I {u>1/2} +M - i 1I {u<1/2} ),
for every u ∈ [0, 1]. For every i ∈ 1 ; n , pathwise uniqueness holds for the following two stochastic differential equations:

du t = u t (1 -u t )dW t (i) + b + i (u t )dt (5.3)
and

du t = u t (1 -u t )dW t (i) + b - i (u t )dt.
(5.4) 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 [START_REF] Ikeda | A comparison theorem for solutions of stochastic differential equations and its applications[END_REF] ensures that the following inequalities hold with probability one:

ξ - t (i) ≤ x t (i) ≤ ξ + t (i) ∀t ≥ 0, ∀i ∈ 1 ; n .
(5.5)

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 t→τ ±,i z (0,1)

ξ ± t (i) = 0 = 1 z exp - x 1/2 2b ± i (u) u(1-u) du dx 1 0 exp - x 1/2 2b ± i (u)
u(1-u) du dx .

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.

3. If 0 < µ 1 < 1/2 and µ 0 > 0 then 0 is attainable: for every 0 < z < b < 1, IP τ ±,i z (0, b) < +∞ and lim

t→τ ±,i z (0,b) ξ ± t (i) = 0 > 0.
4. If µ 1 ≥ 1/2 and µ 0 > 0 then 0 is inaccessible: for every 0 < z < 1,

IP z [∃t > 0, ξ ± t (i) = 0] = 0 and IP z lim t→+∞ ξ ± t (i) = 0 = 0.
Similar properties hold for the behaviour near the point 1.

These properties are sufficient to prove the boundary behaviour claimed for (x t ) t . Since

x t (i) ≤ ξ + t (i) for every t ≥ 0, if (ξ + t (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, since x t (i) ≥ ξ - t (i) for every t ≥ 0, if 0 is inaccessible for (ξ - t (i)) t then 0 is also inaccessible for (x t (i)) t .

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[,

set φ ± i (z) := IP z [lim t→τ ±,i (0,1) ξ ± t (i) = 0]
. 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

greater than φ + i (ǫ) if x ǫ i = ǫ and is greater than 1 -φ - i (1 -ǫ) if x ǫ i = 1 -ǫ.
By the strong Markov property, the probability that (x t ) reaches the boundary in a finite time is greater than min{min(φ + i (ǫ), 1φ - i (1ǫ)), 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

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. Set si,j = s i,js 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 :

g n,µ,s (x) = C n,µ,s n i=1 x 2µ 1 -1 i (1 -x i ) 2µ 0 -1 exp(H n,s (x))
where

• H n,s (x) = 1 2 L⊂ 1;n , |L|≥1 m L (s) ℓ∈L (2x ℓ (1 -x ℓ )) k∈ 1;n \L (1 -2x k (1 -x k )); • 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:

H n,s (x) = L⊂ 1;n , |L|≥1 2 |L|-1 δ L [m(s)](∅) ℓ∈L x ℓ (1 -x ℓ ).
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

G n,s = G n,0 + 1 2 n i=1 x i (1 -x i )∂ i h(x)∂ i where h(x) = L⊂ 1;n , |L|≥1 2 |L|-1 δ L [m(s)](∅) ℓ∈L x ℓ (1 -x ℓ ).
Therefore, as explained in [START_REF] Ethier | Diffusion approximations of the two-locus Wright-Fisher model[END_REF], we can apply a result of [START_REF] Fukushima | Reversibility of solutions to martingale problems[END_REF] to deduce that the diffusion associated with G n,s has a unique reversible stationary distribution m µ,s given by

m µ,s (dx) = C exp(h(x))m µ,0 (dx),
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

h n,µ,s (x) = (2µ -1) n i=1 ln(ρ(x i )) + n-1 ℓ=0 α ℓ L⊂ 1;n , |L|=ℓ+1 k∈L ρ(x k ), where ρ(x i ) = x i (1 -x i ) and α ℓ = 2 ℓ δ (ℓ+1) [m](0).
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).

There are then at least four different types of graph depending on the respective contributions to allelic diversity of mutations (µ > 1/2 or 0 < µ < 1/2) and assortment parameters (m(1)-m(0) smaller than |8µ -4| or not) as shown in Fig. 1. 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).

Proposition 6.2. We assume that the hypotheses H1, H2, H3, H4, H6 and H7 hold. Set

V n = 2µ -1 + 2 -(n+1) n-1 k=0 n-1 k δ (1) [m](k). 1. If V n > 0, then (1/2, . . . , 1/2) is a local maximum of g n,µ,s . 2. If V n < 0, then (1/2, . . . , 1/2) is a local minimum of g n,µ,s . 3. If µ > 1/2 and if δ (1) [m](ℓ) ≥ -(8µ -4) ∀ℓ ∈ 0 ; n -1 , then (1/2, . . . , 1/2) is a global
maximum and is the only critical point of g n,µ,s .

4. If 0 < µ < 1/2 and if δ (1) [m](ℓ) ≤ -(8µ -4) ∀ℓ ∈ 0 ; n -1 , then (1/2, . . . , 1/2) is a
global minimum and is the only critical point of g n,µ,s .

Example 6.1. Let us consider the additive criterion with the assortment sequence s ℓ = bℓ for ℓ ∈ 0 ; n . Then δ (1) [m](ℓ) = 2 -ℓ ℓ ℓ/2 b 1I {ℓ is even} . As 2 -ℓ ℓ ℓ/2 is a strictly decreasing sequence smaller than 1, b < 0 implies V n > 2µ -1 + 1 8 b. Thus, it follows from Proposition 6.2 that if µ > 1/2 and b ≥ -8(2µ -1), the point (1/2, . . . , 1/2) is a local maximum of g n,s,µ . Let us note that if we consider the same sequence s ℓ = bℓ but with the Hamming criterion, then for µ > 1/2 and b < -4(2µ -1), (1/2, . . . , 1/2) is a local minimum of g n,s,µ .

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

V n,i = 2µ -1 + 2 -(n+1) B⊂ 1;n \{i} δ i [m(s)](B)
for every i ∈ 1 ; n and the conditions on δ[m](ℓ) in assertions 3 and 4 are replaced by a condition on δ i [m(s)](A) for every i ∈ 1 ; n and A ∈ P( 1 ; n \ {i}).

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: 1) [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: 1) [m](0) ≥ 0 and δ (1) [m](n -2) > 0.

δ (1) [m](n -1) ≤ δ (1) [m](n -2) ≤ . . . ≤ δ (
δ (1) [m](n -1) ≥ δ (1) [m](n -2) ≥ . . . ≥ δ (
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 1x i . Proposition 6.3. Assume that conditions H1, H2, H3, H4, H6 and H7 hold. Set

V n = 2µ -1 + 2 -(n+1) n-1 k=0 n -1 k δ (1) [m](k).
1. Case µ > 1/2. Assume furthermore that: ii. In [0, 1/2] n , g n,µ,s takes its maximum value at a unique point of the form (ξ 0 , . . . , ξ 0 ).

δ (1) [m](n -1) ≤ δ (1) [m](n -2) ≤ . . . ≤ δ (1) [m](0) ≤ 0 and δ (1) [m](n -2) < 0. (a) If V n > 0 then (1/2, . . . , 1/2) is
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.

2. Case 0 < µ < 1/2. Assume furthermore that: ii. In [0, 1/2] n , g n,µ,s takes its minimum value at a unique point of the form (ξ 0 , . . . , ξ 0 ).

δ (1) [m](n -1) ≥ δ (1) [m](n -2) ≥ . . . ≥ δ (1) [m](0) ≥ 0 and δ (1) [m](n -2) > 0. (a) If V n < 0 then (1/2, . . . , 1/2)
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.

Remark 6.3.

1. ξ 0 = 1/2 -1/2 √ 1 -4λ 0 where λ 0 is the unique solution in ]0, 1/4[ of the equation: 2µ -1 + x n-1 i=0 δ (1) [m](i) n -1 i (2x) i (1 -2x) n-1-i = 0 (E 0 )
More generally, for every ℓ ∈ 0 ; n -1 , ξ ℓ = 1/2 -1/2 √ 1 -4λ ℓ where λ ℓ is the unique solution in ]0, 1/4[ of the equation:

2µ -1 + x n-1 i=0 B n-1,ℓ,i (2x)δ (1) [m](i) = 0 (E ℓ ) and B n,ℓ,i (x) = 2 -ℓ min(i,ℓ) j=max(0,i-n+ℓ) ℓ j n -ℓ i -j x i-j (1 -x) n-ℓ-(i-j) .
Let us note that (B n,ℓ,i (x)) i=0,...,n are positive on ]0, 1[ and their sum is equal to 1.

2.

The assumption that δ (1) [m](i) is a decreasing function of i cannot be removed since one can find examples of assortment parameters satisfying δ (1) [m](i) < 0 for every i ∈ 0 ; n -1 and such that:

(a) µ > 1/2, V n > 0, but (1/2, . . . , 1/2) is not the only local maximum, (b) V n < 0 and g n,µ,s has more than 2 n local maxima.

3. If x i is the proportion of the population with allele 0 at the i-th locus, 2x i (1x 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).

Graphs of the density and simulations of trajectories in the two and three locus cases

Figures 2 to 4 show graphs of the density of the reversible stationary measure in the twolocus case for µ = 0.6 and for several values of s 1s 0 and s 2s 1 , the assortative mating being defined by the Hamming distance. Figures 2 and3 illustrate the two situations considered in Proposition 6.3 when µ > 1/2. When s 1s 0 = 0, the density may have a continuum of critical points as in Fig. 4; this corresponds to a case in which the assumption δ (1) [m](n -2) < 0 of Proposition 6.3 is not satisfied.

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.

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. [START_REF] Schneider | Long-term evolution of polygenic traits under frequency-dependent intraspecific competition[END_REF] 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 and6b).

λ i = ξ i (1 -ξ i )
provides the position of the critical points of index 3i (see Proposition 6.3) and h i is the value of the log-density h n,µ,s at a critical point of index 3i.

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 expressions. We set

ν = 2µ -1, ρ(u) = u(1 -u) for u ∈ [0, 1], ρ(x) = (ρ(x 1 ), . . . , ρ(x n )), h(x) = ν n i=1 log(x i ) + 1 2 n ℓ=1 (m(ℓ) -m(0)) L⊂ 1;n , |L|=ℓ j∈L (2x j ) k∈ 1;n \L (1 -2x k ) 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 , ∂ i h(x) = (1 -2x i )∂ i h(ρ(x))
where

∂ i h(x) = ν x i + n-1 ℓ=0 δ (1) [m](ℓ) L⊂ 1;n \{i},|L|=ℓ j∈L (2x j ) k∈ 1;n \(L∪{i}) (1 -2x k ).
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

∆ = 4ν + 2 -(n-1) n-1 i=0 n -1 i δ (1) [m](i) = 4V n .
This proves the first two assertions of the proposition.

2. The last two assertions follow from the fact that ∂ i h(x) and ∆ are increasing functions of δ (1) [m](ℓ) for every ℓ. Let us prove assertion 3 to illustrate the method. First, if δ (1) [m(s)](ℓ) = -(8µ -4) for every ℓ ∈ 0 ; n -1 then the n coordinates of the diffusion are independent. 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 ∂ i h(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 {ŝ i,j , (i,

j) ∈ A 2 } such that δ (1) [m(ŝ)](ℓ) ≥ δ (1) [m(s)](ℓ)
for every ℓ ∈ 0 ; n -1 .

Proof of Proposition 6.3 We retain the notation introduced in the proof of Proposition 6.2. k+1) [m](0) and denote by e n,k the elementary symmetric polynomial function in n variables of degree k:

For k ∈ 1 ; n , we set α k = 2 k δ (
e n,0 (x) = 1 and e n,k (x) =

L⊂ 1;n , |L|=k ℓ∈L x ℓ for k ∈ 1 ; n .
For instance, e n,1 (x) = x 1 + . . . + x n , e n,2 (x) = 1≤i<j≤n x i x j .

With this notation

h(x) = ν n i=1 ln(x i ) + n-1 ℓ=0 α ℓ e n,ℓ+1 (x).
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(i) = (x 1 , . . . , x i-1 , x i+1 , . . . , x n ) for i ∈ 1 ; n and
x(i,j) = x(j,i) = (x 1 , . . . , x i-1 , x i+1 , . . . , x j-1 , x j+1 , . . . , x n ) for i, j ∈ 1 ; n such that i < j.

Then,

x i e n-1,k (x (i) ) -x j e n-1,k (x (j) ) = (x i -x j )e n-2,k (x (i,j) ). (6.1)
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 , n k=0

2 k δ (k) [a](0)e n,k (x) = n i=0 a i I⊂ 1;n , |I|=i i∈I 2x i j ∈I
(1 -2x j ).

In particular, for every y ∈ IR and ℓ ∈ 0 ; n ,

n k=0 2 k δ (k) [a](0)e n,k ((1/4) ⊗ℓ , y ⊗(n-ℓ) ) = n i=0 a i B n,ℓ,i (2y)
where

B n,ℓ,i (y) = 2 -ℓ min(i,ℓ) j=max(0,i-n+ℓ) ℓ j n -ℓ i -j y i-j (1 -y) n-ℓ-(i-j) .
Proof. See Corollary A.2.

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 ν + ρ(x i ) n-1 k=0 α k e n-1,k ( ρ(x) (i) ) = 0.
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

ρ(x i ) = ρ(x j ) or n-2 k=0 α k e n-2,k ( ρ(x) (i,j) ) = 0. By Lemma 6.2, n-2 k=0 α k e n-2,k (x) = n-2 ℓ=0 δ (1) [a](ℓ)Q ℓ (x),
where Q ℓ denotes a polynomial function which is positive on x ∈]0, 1/4[ n-2 for every ℓ ∈ 0 ; n -2 . Thus this sum cannot vanish in ]0, 1/4[ n-2 under the assumption that all coefficients δ (1) [m](i) have the same sign and that for at

least one i ≤ n -2, δ (1) [m](i)
is non-zero. Therefore, such a critical point exists only if there exists a solution in the interval ]0, 1/4[ of

ν + y n-1 k=0 α k e n-1,k ( 1 4 ) ⊗ℓ , y ⊗(n-ℓ-1) = 0. (E ′ ℓ )
In order to study the solutions of (E ′ ℓ ), let φ ℓ (y) denote the left-hand side of (E ′ ℓ ):

φ ℓ (y) = ν + y n-1 k=0 α k e n-1,k (( 1 4
) ⊗ℓ , y ⊗(n-ℓ-1) ) (6.2) By Lemma 6.2, 1) [m](i). (6.3) Therefore, (E ′ ℓ ) coincides with (E ℓ ) of Remark 6.3. The derivative of φ ℓ is equal to:

φ ℓ (y) = ν + y n-1 i=0 B n-1,ℓ,i (2y)δ (
φ ′ ℓ (y) = n-1 i=0 B n-1,ℓ,i (2y)δ (1) [m](i) + 2y(n -1 -ℓ) n-2 i=0 B n-2,ℓ,i (2y)(δ (1) [m](i + 1) -δ (1) [m](i)). If δ (1) [m](n-1) ≤ • • • ≤ δ (1) [m](0) ≤ 0 (respectively δ (1) [m](n-1) ≥ • • • ≥ δ (1) [m](0) ≥ 0),
φ ℓ is a decreasing function on the interval [0, 1/2] (resp. an increasing function on the interval [0, 1/2]). The value of φ ℓ at 0 is ν and the value at 1/4 is V n . Therefore, under the assumptions of 1 or 2 of the proposition, for every ℓ ∈ {0, . . . , n -1} (E ′ ℓ ) has no solution in ]0, 1/4[ if V n and ν have the same sign and has exactly one solution in ]0, 1/4[ denoted by λ ℓ if V n and ν have opposite signs. This proves assertions 1.(a) and 2.(a).

For every pair of disjoint subsets I and J of 1 ; n , let us introduce the following point:

u I,J = (x 1 , . . . , x n ) with x i =          1/2 if i ∈ I, 1/2 + 1/2 1 -4λ |I| if i ∈ J, 1/2 -1/2 1 -4λ |I| if i ∈ 1 ; n \ (I ∪ J).
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:

a ℓ = ∂ 1 h(( 1 4 ) ⊗ℓ , (λ ℓ ) ⊗(n-ℓ) ), b ℓ = -(1 -4λ ℓ ) ν λ 2 ℓ , c ℓ = (1 -4λ ℓ )∂ 2 n,n-1 h(( 1 4
) ⊗ℓ , (λ ℓ ) ⊗(n-ℓ) ).

The Hessian matrix of h n,µ,s at u I,J is permutation-similar to the following block matrix:

H I,J =     A ℓ 0 0 0 B ℓ,ℓ + C ℓ 0 C ℓ B ℓ,ℓ -    
where

• A ℓ denotes the scalar matrix -2a ℓ I ℓ with a ℓ = ∂ 1 h(( 1 4 ) ⊗ℓ , (λ ℓ ) ⊗(n-ℓ) ),
• B ℓ,k denotes the following k-by-k matrix :

B ℓ,k =        b ℓ c ℓ • • • c ℓ c ℓ . . . . . . . . . . . . . . . . . . c ℓ c ℓ • • • c ℓ b ℓ        ,
• 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.

First, let us study the sign of

a ℓ = 4ν + n-1 i=0 α i e n-1,i (( 1 4 ) ⊗(ℓ-1) , λ ⊗(n-ℓ) ℓ
). As φ ℓ (λ ℓ ) = 0, an application of Lemma 6.1 yields:

a ℓ = (1 -4λ ℓ ) n-2 i=0 α i e n-2,i (( 1 4 ) ⊗(ℓ-1) , λ ⊗(n-1-ℓ) ℓ
).

The right-hand side can be rewritten using Lemma 6.2:

a ℓ = (1 -4λ ℓ ) n-2 i=0 δ (1) [m](i)B n-2,ℓ-1,i (2λ ℓ ).
The conditions on δ (1) [m](i) imply that a ℓ is negative.

Let us now study the coefficients bℓ = (1 -4λ ℓ ) -1 b ℓ and cℓ = (1 -4λ ℓ ) -1 c ℓ . As in the study of a ℓ we use that φ ℓ (λ ℓ ) = 0 and Lemma 6.2 to write bℓ and cℓ in terms of the coefficients δ (1) [m](s)(i):

bℓ = 1 λ ℓ n-1 i=0 δ (1) [m](i)B n-1,ℓ,i (2λ ℓ ), cℓ = 2 n-2 i=0 (δ (1) [m](i + 1) -δ (1) [m](i))B n-2,ℓ,i (2λ ℓ ).
As δ (1) [m(s)](i) is assumed to be a decreasing sequence, cℓ < 0. After some computations, we obtain:

λ ℓ (c ℓ -bℓ ) = - n-2 i=0 δ (1) [m](i)B n-2,ℓ,i (2λ ℓ ).
The conditions on δ (1) [m](i) imply that cℓ > bℓ .

3. Let us prove that 0 < λ n-1 < • • • < λ 0 < 1/4 , which gives the relative positions of the coordinates of the critical points.

Let ℓ ∈ 0 ; n -2 . If we return to the expression (6.2) of φ ℓ , use Lemma 6.1 and then Lemma 6.2, we obtain: 2) [m](i)B n-2,ℓ,i (2y).

φ ℓ+1 (y) -φ ℓ (y) = y(1/4 -y) n-2 i=0 α i+1 e n-2,i ((1/4) ⊗ℓ , y ⊗(n-2-ℓ) ) = 2y(1/4 -y) n-2 i=0 δ (
By assumption, δ (2) [m](i) ≤ 0 for every i ∈ 0 ; n -2 hence φ ℓ+1 (y) ≤ φ ℓ (y) for every y ∈ [0, 1/4]. As the functions φ ℓ are decreasing on [0, 1/4], we deduce that λ ℓ+1 ≤ λ ℓ for every ℓ ∈ 0 ; n -2 . As the two critical points u 1;ℓ ,∅ and u 1;ℓ+1 ,∅ have not the same properties, they cannot coincide and thus λ ℓ+1 < λ ℓ for every ℓ ∈ 0 ; n -2 .

4. Proof of assertion 1.(b).v: let h ℓ denote the value of h n,µ,s at a saddle point of index n-ℓ: -ℓ) ). To prove that h ℓ > h ℓ+1 for every ℓ ∈ 0 ; n -2 , we shall use the properties of the gradient dynamical system dx(t) dt = -∇ h(x) with h = -h n,µ,s . Fix a positive value M large enough so that U M = h-1 ([-M, M ]) contains all critical points of h (such an M exists since h(x) tends to infinity as x tends to the boundary of [0, 1] n ).

p ℓ = ((1/2) ⊗ℓ , (ξ ℓ ) ⊗(n
The function h decreases along trajectories and a trajectory of a point x ∈ M converges to a critical point of h as t tends towards +∞, since h has only isolated critical points.

For k ∈ {0, . . . , n -1}, let U (k) M denote the subset:

U (k) M = {x ∈ U M , x 1 = • • • = x k = 1/2 and x i < 1/2 ∀i > k}.

Every subset U (k)

M contains exactly one critical point, the saddle point p k . As ∂ i h(x) = 0 at points x such that x i = 1/2, the subset U (k) M is positively invariant by the gradient flow. Therefore, to prove that h k > h k+1 , it is enough to show that there exists 0 < y 0 < 1/2 such that for y ∈]y 0 , 1/2[, h((1/2) ⊗k , y, ξ ⊗n-k-1 k+1 ) < h(p k+1 ).

As h((1/2) ⊗k , y, ξ ⊗n-k-1 k+1 ) = -hn,µ,s ((1/4) ⊗k , y(1y), λ ⊗n-k-1 k+1 ), it is enough to show that ∂ k+1 hn,µ,s ((1/4) ⊗(k+1) , λ ⊗(n-k-1) k+1

) < 0. Using that λ k+1 is solution of the equation (E k+1 ), we obtain

∂ k+1 hn,µ,s ((1/4) ⊗(k+1) , λ ⊗(n-k-1) k+1 ) = (1 -4λ k+1 ) n-2 i=0 δ (1) [m](i)B n-2,k,i (2y) < 0.

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.

First, the properties of the generator G n,s stated in assertion (a) of Theorem 4.1 can be obtained by applying the following theorem established by Cerrai and Clément: Theorem 7.1 [START_REF] Cerrai | Well-posedness of the martingale problem for some degenerate diffusion processes occuring in dynamics of populations[END_REF]. Let S + (IR n ) be the space of symmetric, non-negative definite, n × n matrices. Let A : [0, 1] n → S + (IR n ) and b : [0, 1] n → IR n be mappings of class C 2 . For i ∈ {1, . . . , n} and ǫ ∈ {0, 1}, let ν i ǫ denote the unit inward normal vector of the hypercube

C i ǫ = {x ∈ [0, 1] n , x i = ǫ}.
Let us assume the following two conditions:

• for every i ∈ {1, . . . , n}, ǫ ∈ {0, 1} and x ∈ C i ǫ , A(x)ν i ǫ (x) = 0 and b(x), ν i ǫ (x) ≥ 0;

• for every i, j ∈ {1, . . . , n}, A i,j (x) depends only on x i and x j .

Then the operator

L = 1 2 n i=1 n j=1 A i,j (x) ∂ 2 ∂ x i ∂ x j + n i=1 b i (x) ∂ 2 ∂ x i 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 X N k = Φ N (Z N k ) and Y N k = Ψ N (Z N k ) for each k ∈ IN.
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
(a) ǫ -1 N IE z [X N 1 (i) -x(i)] = b i (x, y) + o(1), (b) ǫ -1 N IE z (X N 1 (i) -x(i))(X N 1 (j) -x(j)) = a i,j (x, y) + o(1), (c) ǫ -1 N IE z [(X N 1 (i) -x(i)) 4 ] = o(1), (d) δ -1 N IE z [Y N 1 (ℓ) -y(ℓ)] = c ℓ (x, y) + o(1), (e) δ -1 N IE z [(Y N 1 (ℓ) -y(ℓ)) 2 ] = o(1).
Assume further that (g) The closure of the following operator

L = 1 2 n i,j=1 a i,j (x, 0) ∂ 2 ∂ x i ∂ x j + n i=1 b i (x, 0) ∂ ∂ x i , D(L) = C 2 (F 1 ),
generates a strongly continuous semigroup on C(F 1 ) corresponding to a diffusion process

X in F 1 .
Then the following conclusions in which the symbol ⇒ denotes convergence in distribution, hold:

(i) If X N 0 ⇒ X(0) then {X N [t/ǫ N ] , t ≥ 0} ⇒ X(•) in D F 1 ([0, +∞[) (where D F 1 ([0, +∞[
) is the space of càdlàg paths ω : [0, ∞) → F 1 with the Skorohod topology), (ii) For every positive sequence (t N ) N that converges to +∞, Y N [t N /δ N ] ⇒ 0.

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

E N = {z ∈ (N -1 IN) A , i∈A z(i) = 1}
, and we define by (Φ N , Ψ N ) a change of coordinates such that Ψ -1 N ({0}) is the linkage equilibrium manifold:

Φ N : E N → [0, 1] n and Ψ N : E N → [-1, 1] 2 n -n-1 z → (u 1 , . . . , u n ) z → (u I , I ⊂ 1 ; n s. t. |I| ≥ 2)
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

(N ) 1 = Φ N (Z (N ) 1 ) and Y (N ) 1 = Ψ N (Z (N )
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):

a i,j (x, 0) = x(i)(1 -x(i)) 1I {i=j} , (7.1) b i (x, 0) = (1 -x(i))µ 1 -x(i)µ 0 + (1/2 -x(i))x(i)(1 -x(i))P i,s (x), (7.2)
where

P i,s (x) = J⊂ 1;n \{i} H⊂ 1;n \{i} (s J∪{i},H -s J,H ) j∈J x(j) h∈H x(h) j∈ 1;n , j ∈J ∪{i} (1 -x(j)) h∈ 1;n , h ∈H∪{i} (1 -x(h)),
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 , 1 J ).

In §7.2 we shall show that P i,s has the following two equivalent expressions:

P i,s (x) = A⊂ 1;n \{i} 2 |A| δ A∪{i} [m(s)](∅) ℓ∈A x(ℓ)(1 -x(ℓ)) = A⊂ 1;n \{i} δ i [m(s)](A) k∈A 2x(k)(1 -x(k)) ℓ ∈A∪{i} 1 -2x(ℓ)(1 -x(ℓ)) .

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 IE z (Z (N ) 1 (i) -z(i)) r = N -r j∈A\{i} f N (z, j, i) + (-1) r f N (z, i, j) (7.3) • if r, u ∈ IN * , i, j ∈ A so that i = j, then IE z (Z (N ) 1 (i) -z(i)) r (Z (N ) 1 (j) -z(j)) u = N -(r+u) (-1) r f N (z, i, j) + (-1) u f N (z, j, i) (7.4)
• if r ≥ 3 and i (1) ,. . . ,i (r) ∈ A so that at least three of them are distinct, then

IE z r u=1 Z (N ) 1 (i (u) ) -z(i (u) ) = 0. (7.5) Lemma 7.2. For every u ∈ 1 ; n and z ∈ E N , i∈A, iu=0 B (0) i (z) = 0, (7.6) i∈A, iu=0 B (1) i (z) = (1 -x(u))µ 1 -x(u)µ 0 + 1 2 G u (z), (7.7)
where

x(u) = i∈A, iu=0 z(i) and G u (z) = j∈A h∈A z(j)z(h)s j,h (1 I {ju=0} -x(u)).
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:

i∈A, iu=0 q((j, k); σ (ǫ) u (i)) = 1I {ju=ǫ} +r(u)(1I {ku=ǫ} -1I {ju=ǫ} ) (7.8) with r(u) = I⊂ 1;n \{u} r I = 1 2 by assumption H1. First, formula (7.8) with ǫ = 0 provides i∈A, iu=0 B (0) i (z) = j∈A, ju=0 z(j) + r(u) j∈A k∈A (1I {ku=ǫ} -1I {ju=ǫ} ) - i∈A, iu=0 z(i) = 0. Let B (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:

i∈A, iu=0 B (1,1) i (z) = k∈A j∈A z(j)z(k) i∈A, iu=0 q((j, k); σ (1) u (i))µ 1 -q((j, k); σ (0) u (i))µ 0 .
Applying (7.8) again, we obtain:

i∈A, iu=0 B (1,1) i (z) = (1 -x(u))µ 1 -x(u)µ 0 .
Due to the symmetry of the parameters: s i,j = s j,i for i, j ∈ A, we have:

i∈A, iu=0 B (1,2) i (z) = 0.
Finally, computations using (7.8) yet again yield:

i∈A, iu=0 B (1,3) i (z) = 1 2 G u (z).
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 , z(0 J , 1 J ) = i∈J x(i) i∈ J(1 -x(i)) - I⊂ 1;n s. t. J⊂I, |I|≥2
(-1) |I|-|J| y(I).

(7.9)

Proof. First, by induction on n -|J|, we show that

z(0 J , 1 J ) = I⊂ 1;n s. t. J⊂I (-1) |I|-|J| x(I). (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.

z(0 K , 1 K ) = x(K) - L⊂ 1;n s. t. K L z(0 L , 1L)
We apply the formula (7.10) to every term in the sum and we invert the double sum we have obtained:

z(0 K , 1 K ) = x(K) - H⊂ 1;n s. t. K H x(H) L⊂ 1;n s. t. K L⊂H (-1) |H|-|L| .
The sum between parentheses is equal to

|H|-|K| v=1 (-1) |H|-|K|-v |H| -|K| v = -(-1) |H|-|K| .
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:

I⊂ 1;n , J⊂I (-1) |I|-|J| i∈I x(i) = j∈J x(j) L⊂ 1;n \J (-1) |L| ℓ∈L x(ℓ) = j∈J x(j) i∈ 1;n \J (1 -x(i)).
To shorten the notation, set Lemma 7.7. For a subset I of 1 ; n having at least two elements,

N IE z [Y (N ) 1 (I) -y(I)] = c n,I (x, y) + O(N -1 ) (7.16) where c n,I (x, y) = - L⊂ 1;n , L∩I =∅, L∩I =∅ r L y(I) -1 I {|I|≥4} L⊂ 1;n , |I∩L|≥2,|I∩ L|≥2 r L y(I ∩ L)y(I ∩ L) + 1 I {|I|≥3} L⊂ 1;n , |I∩L|≥2,|I∩ L|≥1 (r L + rL)y(L ∩ I) ℓ∈I∩ L x(ℓ).
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, 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).

By assumption r(k, ℓ) ≥ r(2) > 0, hence the result holds.

• 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. r L e -tr(n) f n,I∩L (t, x, y)f n,I∩ L(t, x, y)

+ 1I {|I|≥3} L⊂ 1;n s. t. |I∩L|≥2,|I∩ L|≥1 (r L + rL)f n,L∩I (t, x, y) ℓ∈I∩ L x(ℓ).
As rI is the probability that the offspring does not inherit all the genes at loci i ∈ I from the same parent, rI ≥ 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.

Conditions (c) and (e). Condition (c) is easy to verify using formulae (7.3), (7.4), (7.5) describing the moments of Z (N ) 1 z. This leads to:

N 2 IE z [(X (N ) 1 (i) -x(i)) 4 ] = O(N -2 ) ∀i ∈ 1 ; n , uniformly on z ∈ E N .
Similarly, using Lemma 7.6, we obtain

N IE z [(Y

Expressions for the drift

We have shown that the i-th coordinate of the drift of the limiting diffusion is (1x(i))µ 1x(i)µ 0 + (1/2x(i))x(i)(1x(i))P i,s (x) where (1x(h)), 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 , 1 J ). The following lemma states that P i,s (x) is actually a polynomial function in the variables x(i)(1x(i)) for i ∈ 1 ; n \ {u}: 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:

P Λ (η)(x) = K⊂Λ\{j} L⊂Λ\{j} k∈K x(k) ℓ∈L x(ℓ) k∈Λ\K (1 -x(k)) h∈Λ\L (1 -x(h))×
x(j) 2 η K∪{j},L∪{j} + (1x(j)) 2 η K,L + x(j)(1x(j))(η K∪{j},L + η K,L∪{j} ) . This expression can be simplified by using that η K∪{j},L∪{j} = η K,L : P Λ (η)(x) = P Λ\{j} (η (0) )(x) + x(j)(1x(j)) P Λ\{j} (η (1) )(x) + P Λ\{j} (η (2) )(x) -2P Λ\{j} (η (0) )(x) , where η (0) , η (1) and η (2) are the following three families of reals indexed by the pairs of subsets of Λ \ {j}: By Lemma 7.9, the expanded form of P i,s as a polynomial function of the n -1 variables x(j)(1x(j)), j = i is:

P i,s (

A Appendix

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.

By taking x i = -1/t for every i ∈ U in Lemma A.1, we can deduce the inverse of the operator S t . This gives a useful formula for inverting a relation between two sequences indexed by the subsets of a finite set.

Corollary A.1. The inverse of the operator S t is S -t , that is (1 + tx j )

f (A) =
where e n,j denotes the elementary polynomial of degree j in n variables:

e n,j (x) =

J⊂ 1;n s. t. |J|=j i∈J

x i .

2. The operator s -t is the inverse of the operator s t :

f (k) = k ℓ=0 k ℓ (-t) k-ℓ s t (f )(ℓ) for every k ∈ 1 ; n .
This corollary provides identities for the forward finite difference operators of any orders since s -1 (f )(k) = δ (k) [f ](0) for every k ∈ 0 ; n . In particular, this leads to the following formula used in the proof of Proposition 5.1: ℓ) [f ](0) = f (k) for every k ∈ 1 ; n (A.3) and Lemma 6.2 used in the proof of Proposition 6.3.

k ℓ=0 k ℓ δ (

A.2 Example 6.2

Under the hypotheses of the assertion 1-(b) of Proposition 6.3, the logarithm of the stationary density h n,s,µ takes its maximum value in [0, 1/2] n at a unique point (ξ 0 , . . . , ξ 0 ) such that λ 0 = ξ 0 (1ξ 0 ) is the unique solution in ]0, 1/4[ of the equation k+1) [m](0) n -1 k y k+1 = 0.

E ′ 0 : 2µ -1 + n-1 k=0 2 k δ (
In [0, 1/2] n the saddle points of index n -1 has one coordinate equal to 1/2 and (n -1) coordinates equal to ξ 1 where λ 1 = ξ 1 (1ξ 1 ) is the unique solution in ]0, 1/4[ of the equation k+1) [m](0) y k+1 = 0.

E ′ 1 : 2µ -1 + n-2 k=0 n -2 k 2 k-1 δ (k+2) [m](0) + 2 k δ (

Figure 1 :

 1 Figure1: Representation of the invariant density g µ,s for the one-locus diffusion when the two mutations rates µ 1 = µ 0 = µ, s 0,0 = s 1,1 = s 0 and s 0,1 = s 1,0 = s 1 . In the figure on the left, µ > 1/2 and matings between individuals of the same allelic type are favoured. The density is bimodal if and only if s 0s 1 > 8µ -4. In the figure on the right, 0 < µ < 1/2 and matings between individuals of different allelic types are favoured. The density tends to +∞ at the boundaries and has a global minimum at 1/2 if and only if s 1s 0 ≤ 4 -8µ.

  do not depend on the locus i and set ν that the coefficients of selection satisfy σ 11,21 = σ 12,22 and σ 11,12 = σ 21,22 (second condition of assumption H4)

  a global maximum and is the only critical point of the density g n,µ,s . (b) If V n < 0 then i. g n,µ,s has a local minimum at (1/2, . . . , 1/2).

  is a global minimum and is the only critical point of the density g n,µ,s .(b) If V n > 0 then i. g n,µ,s has a local maximum at (1/2, . . . , 1/2).

Figure 2 :Figure 3 :Figure 4 :Figure 5 :

 2345 Figure2: Graph of g 2,µ,s when µ = 0.6, s 1s 0 = -0.4 and s 2s 1 = -0.6 so that the point (1/2, 1/2) is the only critical point of the density g 2,s,µ .

Figure 6 :

 6 Figure 6: Simulations of the evolution of the 0-allelic frequency with the three-locus model for two different sets of assortment parameters. The assortative mating favours more strongly pairing between similar types in Fig. 6b. The size of the population is N = 10 3 and the mutation rate is µ = 1. A black dot marks the position of each global maximum of the stationary density, a cross the position of each saddle point of index 2 and a diamond the position of each saddle point of index 1. Some numerical characteristics of the stationary density are presented to the right of each figure: for i ∈ {1, 2, 3}, the value of λ i = ξ i (1ξ i )

  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):

  (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. exists for all (t, x, y) ∈ [0, +∞[×F 1 × F 2 and satisfies lim t→+∞ sup (x,y)∈F 1 ×F 2 |u(t, x, y)| = 0.

  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:v n,I (t, x, y) = exp(-r(n)t)f n,I (t, x, y),where f n,I is a continuous and bounded function on IR ×[0, 1] n × [-1, 1] 2 n -n-1 so that the value of f n,I (t, x, y) depends on x and y only via the coordinates x(i) for i ∈ I and y(J) for J ⊂ I such that |J| ≥ 2.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).

  y(I)) 2 ] = O(N -1 ) ∀I ⊂ 1 ; n , s.t. |I| ≥ 2, uniformly on z ∈ E N .

P

  i,s (x) = J⊂ 1;n \{i} H⊂ 1;n \{i} (s J∪{i},Hs J,H )

  Lemma 7.9. Let Λ be a finite subset of IN. Consider a family of reals β = {β I,J , I, J ⊂ Λ} such that β I,J = β I\J,J\I for every I, J ⊂ Λ. Then, \A .Proof. Let P Λ (β) denote the polynomial function on the right-hand side. The proof is by induction on |Λ|. First, P ∅ (β)(x) = β ∅,∅ = C ∅ (β).

  = η A,B∪{j} for every A, B ⊂ Λ \ {j}.The inductive hypothesis applies to Λ \ {j} and the three families of reals η (0) , η (1) and η (2) :P Λ (η)(x) = L⊂Λ\{j} C L (η) ℓ∈L x(ℓ)(1x(ℓ)) + A∪{j},T \A + η A,(T ∪{j})\A -2η A,T \A ).The double sum of the terms η A∪{j},T \A + η A,(T ∪{j})\A is equal to:T ⊂L, j∈T (-2) |L|-|T | A⊂T η A,T \A .Therefore, CL = C L (η) and P Λ (η)(x) = L⊂Λ C L (η) ℓ∈L x(ℓ)(1x(ℓ)) which completes the proof by induction.

  A∪{i},T \As A,T \A ).The coefficient α i,L (s) can be rewritten in terms of the mean values of the assortment parameters m T (s) for T ⊂ L:α i,L (s) = 2 |L| T ⊂L (-1) |L|-|T | (m T ∪{i} (s)m T (s)) = 2 |L| T ⊂L (-1) |L|-|T | δ i [m(s)](T ).Indeed, it follows from the assumption H4 that for every i ∈ 1 ; n and T ⊂ 1 ; n \ {i},m T (s) = 2 -|T | A⊂T s A,T \A and m T ∪{i} (s) = 2 -|T | A⊂T s A∪{i},T \A .Using formula (4.5), we obtainα i,L (s) = 2 |L| δ L∪{i} [m(s)](∅).The following factorised form of the polynomial function P i,s can be derived from a general identity stated in Lemma A.1:

S

  t (f )(A) = B⊂A t |A|-|B| f (B) for every A ∈ P(E)(with the usual convention a 0 = 1 for every a ∈ IR). Most of the combinatorial formulae used in the paper can be deduced from this general identity: Lemma A.1. Let U be a subset of E and let {x u , u ∈ U } be a family of reals. One way to derive this equality is to interchange the sum on the right-hand side of the equation with the sum that appears in the definition of S t (f )(A), to use the new summation index C = A \ B and to recognize the following expansion of the product of the terms 1 + tx i :As S -1 (f )(A) is nothing other than δ A [f ](∅) by (4.5), if we apply Lemma A.1 with t = -1, f (A) = δ i [m(s)](A)and the family of reals {2x(j)(1x(j)), j ∈ 1 ; n \ {i}}, we obtain the following equalityA⊂ 1;n \{i} 2 |A| δ A∪{i} [m(s)](∅) ℓ∈A x(ℓ)(1x(ℓ)) = A⊂ 1;n \{i} δ i [m(s)](A) k∈A 2x(k)(1x(k)) ℓ ∈A∪{i} 1 -2x(ℓ)(1x(ℓ)) .

B⊂A(

  -t) |A|-|B| S t (f )(B) for every A ⊂ E.From Corollary A.1 we can deduce the following identity for the finite difference operator:f (A) = B⊂A δ B [f ](∅) for every A ∈ P(E). (A.2)By considering the operator S t for a function f which is constant on subsets having the same number of elements, we can rewrite the previous relations to obtain useful formulae relating two sequences indexed by the integers 0, 1, . . . , n.Corollary A.2. Let t be a real number. Let n ∈ IN * . For a function f defined on 0 ; n , let s t (f ) be the function defined by:

  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

  1 4 (s 00,00 + s 01,01 + s 10,10 + s 11,11 ), 00,01 + s 01,00 + s 11,10 + s 10,11

	m {1} (s) =	1 4	(s 00,10 + s 10,00 + s 01,11 + s 11,01 ),
	m {2} (s) =	1 4	(s

A diploid individual is homozygous at a gene locus when its cells contain two identical alleles at the locus.

AME supported in part by EPSRC grant no. EP/G052026/1

Condition (a). To show that condition (a) of Theorem 7.2 holds, we first examine the drift of Z (N ) . A Taylor expansion of the transition probabilities of the Markov chain (Z (N ) t

) t∈IN using assumption H2 yields the following formula: Lemma 7.1. For every i ∈ A,

q((j, k); (1i u , i 1;n \{u} ))µ 1-iuq((j, k); i)

Proof. By assumption H2, for two different types i, j ∈ A

where

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 (z) and i∈A,iu=0 B

(1) i (z). Direct computations yield:

• s I,J = s i,j for i = (0 I , 1Ī ) and j = (0 J , 1 J ).

With this notation, for every J ⊂ Λ u ,

• z(0

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 ,

uniformly on z ∈ E N where

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 (N ) 1

x, showing that condition (b) holds:

Using formulae (7.3) and (7.4) and assumption H2, we obtain

With the convention x({i, j}) = x(i) if i = j, we have T

(2) i,j = x({i, j}) and it follows from assumption H1 (r I = rĪ for every

Therefore, for every i, j ∈ 1 ; n ,

and

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 moments of Z (N ) 1

z.

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:

). (7.13) 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

To obtain an expression for IE z [Y (N ) 1

(I)y(I)] in the new coordinates, it remains to replace each term x(L) for |L| ≥ 2 with ℓ∈L x(ℓ)y(L) in (7.15). This leads to the following lemma and shows that condition (d) holds.

If we denote by h n,i the value of h n,s,µ at a critical point of index ni then

If we define the assortment by means of the Hamming criterion with the quadratic sequence of parameters: s k = s 0 -(bk + ck 2 ) ∀k ∈ 0 ; n with c > 0 and b + c > 0, then 2) [m](0) = -2c and δ (r) [m](0) = 0 ∀r ≥ 3.

In this case, λ 0 and λ 1 are solutions of quadratic functions: 2µ-1-(b+c)λ 0 -4c(n-1)λ 2 0 = 0 and 2µ-1-(b+2c)λ 1 -4c(n-2)λ 2 1 = 0. After some computations, we obtain:

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.

Proof. Let Q n,k,a denote the quadratic form with matrix M n,k (a) in the canonical basis. For every x ∈ IR n , Q n,k,a (x) = n i=1 x 2 i + 2a 1≤i<j≤n ǫ i ǫ j x i x j , where ǫ 1 = . . . = ǫ k = 1 and ǫ k+1 = . . . = ǫ n = -1. This lemma can be established by induction on n by using the following decomposition of Q n,k,a (x):

1≤i<j≤n-1

where b = a 1+a ∈ [0, 1[.