Bahram Houchmandzadeh

Neutral Aggregation in Finite Length Genotype space

The advent of modern genome sequencing techniques allows for a more stringent test of the neutrality hypothesis of Darwinian evolution, where all individuals have the same fitness. Using the individual based model of Wright and Fisher, we compute the amplitude of neutral aggregation in the genome space, i.e., the probability of finding two individuals at genetic (hamming) distance k as a function of genome size L, population size N and mutation probability per base ν. In well mixed populations, we show that for N ν < 1/L, neutral aggregation is the dominant force and most individuals are found at short genetic distances from each other. For N ν > 1 on the contrary, individuals are randomly dispersed in genome space. The results are extended to geographically dispersed population, where the controlling parameter is shown to be a combination of mutation and migration probability. The theory we develop can be used to test the neutrality hypothesis in various ecological and evolutionary systems.

I. INTRODUCTION.

Aggregation of individuals is a common observation in evolutionary and ecological systems. By aggregation, we mean the observation that some areas of space contain large numbers of individuals while other parts contain relatively few. To be more precise, the variance of population distribution is much larger than its mean. Consider for example ecological communities where observations are made in the real space. Since the seminal work of Taylor et al., [START_REF] Taylor | The density dependence of spatial behaviour and the rarity of randomness[END_REF] who surveyed around 4000 samples from 100 species across different kingdoms, it has been established that all species tend to spatially aggregate.

There are many causes for aggregation and determining these factors is the fundamental subject of ecological and evolutionary theories and of our understanding of the natural world. One cause that is often disregarded is encoded in the very nature of life: individuals appear by birth close to their parents, but can die anywhere. Therefore, each birth event enriches the spatial pair correlation function at short distances, while a death event depletes all distances. Of course, individuals (or their seed for plants) move randomly in space and the diffusion phenomena tends to counteract and dilute the effect of correlation created at short distance by birth. It can be shown however [START_REF] Young | Reproductive pair correlations and the clustering of organisms[END_REF][START_REF] Houchmandzadeh | Clustering of diffusing organisms[END_REF][START_REF] Houchmandzadeh | Clustering in neutral ecology[END_REF][START_REF] Houchmandzadeh | Theory of neutral clustering for growing populations[END_REF]that diffusion is not enough to efficiently dilute the correlation creation at short distances for spatial dimensions d ≤ 3 and large system size extension. This phenomena is called neutral clustering (see [START_REF] Korolev | Genetic demixing and evolution in linear stepping stone models[END_REF] for a review); it has been demonstrated experimentally [START_REF] Houchmandzadeh | Neutral clustering in a simple experimental ecological community[END_REF] and the main concept has been applied to other systems such as neutrons in nuclear reactors [START_REF] Dumonteil | Particle clustering in Monte Carlo criticality simulations[END_REF][START_REF] Houchmandzadeh | Neutron fluctuations: The importance of being delayed[END_REF] and evolution of bimolecular networks [START_REF] Stein | Logistic map analysis of biomolecular network evolution[END_REF].

The same arguments can be applied to genome space. Consider individuals characterized by the sequence of their genome, of length L. Each duplication event can give rise to a new individual due to mutations at one position in the sequence. The genetic distance k between two individuals being defined as the number of differences between their sequence (the hamming distance), we see that birth events again tend to enrich correlations at short distances, under the assumptions that all mutations are neutral, i.e., don't affect the fitness of the individual. There is however a marked difference between diffusion in real space due to random movements and diffusion in genotype space due to mutations. The former happens in a low dimensional space (d ≤ 3) and large extension; the latter occurs in high dimensional space (d = L) but small extension (the number of values such as A,T,C,G that a position along the sequence can take). Therefore, it is not clear a priori which effect (neutral aggregation or diffusion) is dominant. The purpose of this article is to weight these factors precisely in genotype space.

In the evolutionary field, neutral models were first proposed by Kimura [START_REF] Kimura | The Neutral Theory of Molecular Evolution[END_REF]as the main driving force in evolution. In the ecology field, Hubbell [START_REF] Stephen | The unified neutral theory of Biodiversity and Biogeography[END_REF] used a special version of the neutral model (UNTB), called the infinite allele model [START_REF] Kimura | The number of alleles that can be maintained in a finite population[END_REF], to explain the pattern of biodiversity in nature in terms of neutral mutations . Both theories are passionately debated in the literature (see [START_REF] Azaele | Statistical mechanics of ecological systems: Neutral theory and beyond[END_REF] for a review of UNTB, [START_REF] Linquist | Applying ecological models to communities of genetic elements: the case of neutral theory[END_REF] for a critical review of both theories and their interconnections, or [START_REF] Rosindell | The unified neutral theory of biodiversity and biogeography at age ten[END_REF] defending the importance of neutral theories in ecology).

The advent of modern gene sequencing tools has opened the possibility for more stringent tests of the neutral hypothesis in ecological communities, combining the standard measurement of abundances of species with the histogram of genetic distances between species. Jeraldo et al. [START_REF] Jeraldo | Quantification of the relative roles of niche and neutral processes in structuring gastrointestinal microbiomes[END_REF] for example measured the genetic distance between OTUs of six gastrointestinal microbiomes of different mammals and found that the histogram of distances is sharply peaked toward short distances. They thus concluded that neutral processes play a negligible role in these communities and selection is the dominant force. This argument of Jeraldo et al. is however problematic, as they equate neutral with "drawn at random", ignoring the importance of neutral forces discussed above. This point was raised by D'Andrea and Ostling [START_REF] Rafael | Can clusterting in genotype space reveal niches ?[END_REF] who developed a simple, interaction free model to demonstrate that neutral causes for this model will lead to genetic distance histograms that are peaked at short distances.

We will show here, by exactly computing the distance pair correlation function (normalized histogram) under the neutral hypothesis, that the picture is more nuanced and depends crucially on the mutation number, i.e. the product of mutation probability and the population size. We show that in well mixed populations, for small mutation number, neutral aggregation is the dominant force and the distribution of distances is nearly geometric, showing a sharp peak at the origin. In this limit, most individuals are very close to each other in genotype space. For large mutation numbers however, the distribution of distances is more binomial-like and shows a peak at ∼ L/2, as was assumed by Jeraldo et al. (figure 5).

A further complication is that real populations are geographically distributed and spatial migration also plays an important role. Histogram measurements have thus to take into account this factor in order to weight the effects of neutral factors.

The aim of this article is to compute precisely the genetic pair correlation function under the neutral hypothesis, using an individual based model. This quantity is of prime importance in population genetics [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF] and is nowadays investigated mainly by coalescent theory [START_REF] Hein | Gene Genealogies, Variation and Evolution: A primer in Coalescent Theory[END_REF]. However, application of coalescent theory to finite sequences and to geographically dispersed populations is rather difficult. In contrast, the method we develop below uses only straightforward mathematical tools and the results are derived by simple means. The computations presented below are confirmed by individual based numerical simulations.

This article is organized as follow: in the next section, we define the model, the pair correlation function u k and summarize the main findings of the article. The next three sections are devoted to the mathematical derivations of these results: in section III, we use the well known "infinite site" model to illustrate how the equations governing u k can be obtained directly from the model and be solved. Section IV generalizes this approach to finite sequences when back mutations are explicitly taken into account and the u k are found by their probability generating function. This section constitutes the heart of this article. Section V generalizes this approach to take into account the geographic extension of populations and the influence of migrations. The final section is devoted to placing this work in perspective and to the concluding remarks.

These computations should constitute a useful tool in assessing the importance of neutral phenomena in ecological communities in general, where results analogous to that obtained by Jeraldo et al. can be evaluated.

II. MODEL DEFINITION AND MAIN RESULTS.

In this section, we define the model and summarize the main results of the manuscripts. The mathematical derivations of these results are given in the following The habitat contains N haploid individuals with various alleles of a gene of length L. At generation Gt, each individual produces q progeny, some of whom may contain a mutation with respect to their parent (marked by a star). N individuals are selected among the N q progeny to constitute the Gt+1 generation.

sections.

The two fundamental individual based models of population genetics are the non-overlapping generations model of Wright-Fisher (WF) [START_REF] Fisher | The genetical theory of natural selection, a complete variorum edition[END_REF] and the continuous time model of Moran [START_REF] Moran | The Statistical processes of evolutionary theory[END_REF]. These models are equivalent and obey the same diffusion equation in the large population limit [START_REF] Houchmandzadeh | Alternative to the diffusion equation in population genetics[END_REF]. Their use is dictated by the respective ease by which relevant quantities are obtained. For the present computations, we use the WF model which allows for the straightforward computation of the pair correlation function and its extension to the spatial case.

Consider a habitat containing N haploid individuals at each generation (figure 1). The individuals are characterized by the sequence of their genome of length L. Without loss of generality, we will assume that each base can take only two values, 0 and 1: When there are 2 n possible values for each base (such as A, T, C, G where n = 2 and 2 n = 4), we can map the problem to a binary system with sequence length nL.

At each generation, each individual produces q > 1 progeny (the model can be trivially generalized to a random number of progeny). The progeny are then sampled at random to constitute the N individuals of the next generation (figure 1). We assume the system to be neutral : all individuals have equal fitness, and each progeny, regardless of its genome, has the same probability of getting to the next generation.

Each progeny can differ from its parent because of mutations. In the following, we will assume the mutation probability per base ν to be small and neglect the probability of having two mutations on the same progeny. This limitation can be relaxed if needed. We further assume ν to be the same for all bases of the genome sequence, regardless of their position in the sequence and the value of neighboring sites. This mutation model, known as K80, was first proposed by Kimura [START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF] and various refinements of this model are widely used in the field of molecular evolution, for example to deduce phylogenetic trees from sequence data [START_REF] Yang | Computational molecular evolution[END_REF].

The genetic distance between two individuals is defined as the number of bases in which they differ (hamming distance). For example, the two individuals shown in figure 2 have a hamming distance of k = 4.

The quantity we compute in the following sections is u k , i.e., the probability that, for large times, two individuals drawn at random are at distance k from each other. The computations are done in the limit of small mutation probability per sequence λ = Lν 1. For a gene of length L = 1000 base pair and mutation probability per base of ν = 10 -10 , this is a reasonable hypothesis. On the other hand, for viruses with high mutation probabilities such as ν = 10 -4 , this approximation is more problematic and the model has to be extended to take into account higher order perturbations.

If the distribution u k is peaked at k = 0 and is a fast decreasing function of k, the system is clustered: there is one dominant sequence in the population, and most other sequences are at short distances from the dominant one. On the other hand, if sequences were distributed totally randomly, we would expect u k to have a binomial distribution

u k = 2 -L L k (1)
and display a peak at k = L/2. We find that the actual distribution u k depends only on two parameters: the sequence length L and the mutation number per base

Ω = 2νN (2)
The mutation number combines in a single number the contribution of mutations and the population size. We find that for small mutation numbers Ω < 1/L, the population is clustered, u k is peaked at k = 0 and decreases nearly exponentially with k. On the other hand, for Ω > 1, the distribution of distances u k tends toward a binomial one, showing peaks at k > 1 (Figure 5). To measure experimentally u k necessitates a large quantity of data. More robust quantities are the moments of the distribution such as the mean distance k =

ku k and the variance

V = k 2 u k -k 2 .
The moments of the distribution contain the same information as the distribution itself. The results for the moments are particularly simple. For example, the mean is given by

k L ≈ Ω 1 + 2Ω (3)
and all the other moments are deduced from the mean by simple recurrence ( equation 33). The probability distribution u k can be determined from these moments (equation 36). We see again from expression (3) that for Ω < 1/L, k < 1 and that most sequences are located at very short distances from the dominant one. On the other hand, for Ω 1, k ≈ L/2, a result which is expected from the binomial distribution of u k (relation 1) and which denotes a random distribution of sequences.

We observe that the distribution depends crucially on the size of the population N . For large populations, the hypothesis of well mixed populations ceases to be valid. We must explicitly take into account physical space and migrations in order to determine the effective size of the population.

For this purpose, we resort to the classical stepping stone model (see [START_REF] Korolev | Genetic demixing and evolution in linear stepping stone models[END_REF] for a review), dividing the space into demes (or patches) of N individuals, and denote by m the migration probability, i.e. the probability that an individual in a deme at generation t + 1 has its parent in a neighboring deme at generation t (figure 3). The migration number is defined as M = N m and combines into a single number the size of the deme and the migration probability.

Let y be the number of patches between two demes. We compute u(k, y), the probability that two individuals drawn at random in two demes at physical distance y are at genetic distance k. We restrict the computation to the mean hamming distance

k(y) = L k=0 ku(k, y)
We show that the mean hamming distance in the same patch is similar to expression ( 3)

k(0) L = Ω e 1 + 2Ω e (4)
Where the effective mutation number Ω e is defined as

Ω e = Ω 2 + 2ΩM
The amount of neutral clustering depends crucially on competition between migration and mutations. Moreover, k(y) tends exponentially in y toward L/2 (relation 47).

In the following sections, we detail the mathematical derivations of the above results and provide a rough estimation of the expected clustering for gut bacteria.

III. INFINITE SEQUENCE MODEL.

Consider a well mixed population of N individuals, where each individual is characterized by the sequence of its genome of length L; the population follows a Wright-Fisher competition (figure 1). At each reproduction event, a base (a position along the sequence) can be mutated with probability ν. We define the mutation probability per sequence as λ = Lν and we suppose λ 1 in this article (see below for the limit of validity of this hypothesis).

We wish to compute the probability u k (t + 1) that two individuals drawn at random at generation G t+1 are at (genetic) distance k, i.e. they differ in exactly k bases.

Before doing the full computation, let us consider the classical "infinite site" model [START_REF] Watterson | On the number of segregating sites in genetical models without recombination[END_REF][START_REF] Tajima | Infinite-allele model and infinite-site model in population genetics[END_REF], where the main mathematical concepts are straightforward to introduce. The infinite site model assumes that each mutation gives rise to a new individual that does not already exist in the habitat; this consists in assuming L → ∞ while λ remains finite.

We denote

a = q -1 N q -1 (5)
as the probability that two individuals drawn at random are from the same parent, where q is the number of progeny of each individual. Note that for q 1, a ≈ 1/N .

The probability that two individuals are at distance k = 0 from each other is :

1. The two individuals are from the same parent (probability a) AND no mutation has taken place in them (probability (1 -λ) 2 ) OR 2. The two individuals are from different parents which themselves are genetically identical (probability (1 -a)u 0 (t) ) AND no mutation has taken place in them (probability (1 -λ) 2 ).

In other words,

u 0 (t + 1) = (1 -λ) 2 [a + (1 -a)u 0 (t)] (6)
Following the same line of argument (see appendix A), we can write generally the linear system:

u 0 (t + 1) = (1 -a)Bu 0 (t) + b 0 (7) u 1 (t + 1) = (1 -a) (Au 0 (t) + Bu 1 (t)) + b 1 (8) u k (t + 1) = (1 -a) (Au k-1 (t) + Bu k (t)) (k > 1) (9)
where

A = 2λ(1 -λ), B = (1 -λ) 2 , b 0 = (1 -λ) 2 a and b 1 = 2λ(1 -λ)a.
Note that the probability that two individuals from the same parent are at distance k = 2 is aλ 2 and therefore it has been neglected in relation [START_REF] Houchmandzadeh | Neutron fluctuations: The importance of being delayed[END_REF]. It should also be noted that because we suppose λ 1,

A/B = b 1 /b 0 ∼ 2λ 1.
The coefficient A captures mutation events that increase the genetic distances, while B relates to events that maintain the genetic distance.

We define

Φ = (1 -a)A 1 -(1 -a)B ≈ 2λ a + 2λ = Θ 1 + Θ ∈ [0, 1)
where

Θ = 2λ/a ≈ 2λN (10)
is the per sequence mutation number. The stationary solution u k for large times of relations [START_REF] Houchmandzadeh | Neutral clustering in a simple experimental ecological community[END_REF][START_REF] Dumonteil | Particle clustering in Monte Carlo criticality simulations[END_REF][START_REF] Houchmandzadeh | Neutron fluctuations: The importance of being delayed[END_REF], which is a one term recurrence relation, is:

u k = Φ k (u 0 + a 1 -a ) (k ≥ 1) (11)
≈ Φ k (u 0 + 1 N ) (12)
where

u 0 = b 0 1 -(1 -a)B ≈ a a + 2λ = 1 1 + Θ (13)
Note that for the infinite allele model, the autocorrelation function is always peaked at k = 0 (or at k = 1 for Θ 1). In the low mutation regime (Θ 1), nearly all individuals are identical (u 0 ≈ 1 -Θ), and u k drops sharply as a function of k (Φ 1). On the other hand, in the high mutation regime Θ 1, the distribution is nearly flat for a wide range of k, and many different genomes coexist. The mean genetic distance for this model is :

k = ∞ k=1 ku k ≈ Θ(1 + aΘ)
In the above computation, we have assumed λ 1 and therefore have neglected terms in λ 2 . We can evaluate the accuracy of this approximation by estimating S(t) = k u k (t) and its deviation from unity. Summing the lines of relations (7-9), we have

S(t + 1) = (1 -a)(A + B)S(t) + (b 0 + b 1 ) = (1 -λ 2 ) [(1 -a)S(t) + a]
and for large times, the stationary value of S is

S ≈ a a + λ 2 ≈ 1 - λ 2 a ≈ 1 -N λ 2
Therefore the approximation we are making is valid when N λ 2 1. For gene sequences of length L ∼ 10 3 and mutation probability per base ν ≈ 10 -9 , the limit of validity of this computation is N 10 12 . The approximation is therefore valid for large populations, especially when we consider spatially distributed populations as in section V.

IV. FINITE SEQUENCE LENGTH.

We now investigate the case of finite sequence length L, which is the main point of this article. Note that if sequences were generated randomly with equal probabilities for 0 and 1 values of the bases, the probabilities u k would follow the binomial distribution

u k = 2 -L L k ( 14
)
where the mean distance is k = L/2 and the variance V = L/4. This is the implicit hypothesis used by Jeraldo et al. to exclude neutral processes. However sequences are not generated randomly, but are inherited from parents with a small probability for mutations. This fact radically changes the expected distribution of distances as we show below.

A. The evolution equation.

In the infinite site model considered in the preceding section, we captured the basic equations (7-9) in terms of events that increase the genetic distance (coefficient A) and events that maintain the genetic distance (coefficient B). For finite sequences, the picture is more complicated. First, mutations can also decrease genetic distances if they appear in bases which are already different between two sequences (back mutations); we would capture these events in a coefficient C. Second, the coefficients are not constant but depend on the number of bases in which two individuals already differ.

Consider two sequences of length L that differ in k bases (figure 2). The probability that an event increases their difference by one unit is that of one mutation ( probability 2λ(1 -λ) ) occurring in one of their identical bases (probability (L -k)/L ) :

A k = 2λ(1 -λ)(1 - k L ) = 2λ(1 - k L ) + O(λ 2 ) (15)
By the same token, the probability of decreasing their distance by one unit is

C k = 2λ(1 -λ) k L = 2λ k L + O(λ 2 ) (16)
The probability that their distance is conserved is

B k = (1 -λ) 2 + 2λ 2 k(L -k) L 2 = 1 -2λ + O(λ 2 ) (17)
and to the first order in λ, A k + B k + C k = 1. Now, following the same line of argument as in the preceding section, the probability u k (t + 1) of finding two individuals at distance k at time t + 1 is

u 0 (t + 1) = (1 -a) (B 0 u 0 (t) + C 1 u 1 (t)) + b 0 (18) u 1 (t + 1) = (1 -a) (A 0 u 0 (t) + B 1 u 1 (t) + C 2 u 2 (t)) + b 1 (19) u k (t + 1) = (1 -a) (A k-1 u k-1 + B k u k + C k+1 u k+1 )(20) u L (t + 1) = (1 -a)(A L-1 u L-1 + B L u L ) (21)
where b 0 and b 1 are defined as before and 1 < k < L in equation [START_REF] Hein | Gene Genealogies, Variation and Evolution: A primer in Coalescent Theory[END_REF]. Note that obviously, u k = 0 for k > L. The above set of linear equations can be written in vectorial notation as

|u(t + 1) = (1 -a)Q |u(t) + |b (22)
where |u(t) = (u 0 (t), u 1 (t), ..., u L (t)) T ,

Q =            B 0 C 1 0 0 • • • 0 A 0 B 1 C 2 0 • • • 0 0 A 1 B 2 C 3 . . . . . . 0 0 . . . . . . . . . 0 . . . . . . . . . . . . C L 0 • • • 0 A L-1 B L            (23)
and |b = (b 0 , b 1 , 0, . . . , 0) T . The stationary probabilities |u for large times are obtained from

(I -(1 -a)Q) |u = |b ( 24
)
where I is the identity matrix. The above relation is a (L + 1) × (L + 1) linear system that can be solved numerically for practical purposes. We are of course interested in its analytical solution.

B. The probability generating function.

The non-homogeneous linear system ( 22) is a two term recurrence relation where the coefficients are not constant. The stationary solution u k for large time however is surprisingly simple if we use the probability generating function (PGF)

φ(z) = L k=0 u k z k (25)
The PGF contains the most complete information on the system; the probabilities and their moments are obtained from the derivatives of φ at either z = 0 or z = 1. Let

k (n) = k(k -1)...(k -n + 1) (26)
be the factorial moment of order n. For example, the usual mean distance is k = k (1) and the variance is

V = k (2) + k -k
2 and so on. The factorial moments are given by

k (n) = d n φ dz n z=1 (27)
The moments, specially the lower ones, are the most robust quantities that we can estimate from real data analysis. Below, we shall also use the normalized factorial moments µ n = k (n) /n! . On the other hand, the probabilities are

u k = 1 k! d k φ dz k z=0 (28)
although experimentally, their estimations necessitates much more data than what is needed for (lower) moment estimation.

It can be shown that the PGF obeys a simple first order differential equation (see appendix B):

(1-z 2 )φ + Lz -(1 + 1 Θ )L φ = -aL(z -1)- L Θ (29)
The term aL ≈ L/N weights the relative importance of the sequence length compared to the population size. The above differential equation can be exactly solved in terms of the hypergeometric function. However, as we are interested only in the probabilities u k and their moments, we don't even need to solve equation ( 29) and we can extract all the moments from simple arguments as discussed below. Before the full discussion, note that N → ∞ implies Θ → ∞ and a → 0. In this high mutation number regime, equation ( 29) becomes

(1 -z 2 )φ + L(z -1)φ = 0
with the obvious solution φ(z) = 2 -L (z + 1) L satisfying the initial condition φ(1) = 1. This is the PGF for the binomial distribution, which is expected if sequences were drawn at random. Therefore, in the very high mutation number regime (Θ L), the distribution of distances indeed becomes binomial.

On the other hand, the infinite sites results can be recovered from equation ( 29) by letting L → ∞. In this limit, the φ term becomes negligible in equation ( 29) and the PGF is simply

φ(z) = 1 -aΘ + aΘz 1 + Θ -Θz ≈ 1 + aΘz 1 + Θ -Θz (30)
where aΘ = 2λ has been neglected compared to one, as we suppose λ 1. Relation ( 30) is the PGF of the probabilities u k computed for infinite sites model in the preceding section (equation 12). This expression was first computed by Watterson[19] for the infinite allele model.

C. Finding the moments and probabilities.

Let us first consider the point z = 1 in equation ( 29). The first term vanishes at this point and we have trivially φ(1) = 1 which just states that the sum of the probabilities is unity:

φ(1) = L k=0 u k = 1
We see here that computing φ(1) does not require a knowledge of φ (1). This is a general feature of the PGF equation ( 29): φ (n) (1) does not depend on the φ (n+1) (1) and we can deduce all the moments from a hierarchical structure. To see this, we differentiate equation ( 29) in respect to z:

(1 -z 2 )φ" + (L -2)z -(1 + 1 Θ )L φ + Lφ = -aL (31
) As before, the higher order term vanishes at z = 1 and therefore the mean distance between individuals is

k = φ (1) = LΘ(a + 1) L + 2Θ ≈ LΘ L + 2Θ (32)
we see here that for mutation numbers Θ < 1, the mean distance k < 1 and most individuals are clustered very close to each other in genetic space. Only for high mutation numbers Θ = O(L), does the mean distance between individuals become appreciable. For example, Θ = L/2 results in k = L/4. For very large mutation numbers Θ L we reach k = L/2, as in the binomial distribution.

Obtaining higher moments follows the same procedure. Applying d/dz to equation [START_REF] James | Backbones of evolutionary history test biodiversity theory for microbes[END_REF] and computing φ (1), we find

k (2) = 2L(L -1)Θ 2 (L + 2Θ)(L + 4Θ)
For the low mutation regime Θ 1 and long sequences L 1, k (2) ≈ 2Θ 2 and the distance distribution is indeed sharply centered around the origin. For the high mutation regime Θ L, we find k (2) = L(L -1)/4, as expected from the binomial distribution.

Successive differentiation allows us to obtain all the moments. For n ≥ 2

k (n) = n(L -n + 1)Θ L + 2Θn k (n-1) (33)
which leads to

µ n = 1 n! k (n) = 2 -n L(L -1)...(L -n + 1) (γ + 1)(γ + 2)...(γ + n) (34) = 2 -n (L) (n) (γ + n) (n) (n ≥ 1) (35)
where γ = L/(2Θ) and (r) (n) is the descending Pochhammer symbol

(r) (n) = r(r -1)...(r -n + 1) = Γ(r + 1) Γ(r + 1 -n)
The expression [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF]35) for k (n) should be corrected by the factor (1+a) for very small populations, as in relation [START_REF] Gardner | The genetical theory of kin selection[END_REF].

Figure 4 shows the results for normalized factorial moments µ n obtained from relation (35) and their excellent agreement with the moments obtained by numerical resolution of the linear set of equations [START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF].

To find the probabilities u k , we can rearrange the PGF function Expanding (z -1) n and identifying the corresponding powers of z in both sums, we find

φ(z) = L n=0 u n z n = L n=0 µ n (z -1) n
u = L n= (-1) n-µ n n (36)
Defining the matrix C such that (C) n = (-1) n- n , the solution in vectorial notation is

|u = C |µ (37)
Figure 5 shows the excellent agreement between the above theoretical results and individual based stochastic simulations of the neutral model (see appendix D).

It can be observed that as Θ increases , u k transforms from a sharply decreasing distribution with its peak at zero toward a binomial distribution with the peak approaching L/2. The curves cease to decrease monotonically as a function of the genetic distance k for Θ ≈ L/2.

For small mutation numbers Θ 1, µ n given by expression (35) can be approximated by

µ n = (1 + a) Θ L n (L) (n) n ≥ 1 (38)
which decreases faster than exponentially. In this case, only the first term contributes to the sum (36) and we have

u 0 = 1 -(1 + a)Θ u k = µ k (k ≥ 1)
It has been checked numerically that this approximation of the probabilities u k is very good for Θ < 0.1.

To summarize this section, we have obtained the exact solution for the moments µ n (relation 35) and probabilities u k (relation 37) of the finite sequence neutral WF model; we have obtained various limiting expressions (large and small Θ, large L) and have shown the accuracy of these expressions by comparing them with numerical simulations. The parameter Ω that controls the behavior of the system is the ratio of mutation number to sequence length

Ω= Θ L = 2N ν (39)
i.e., the per base mutation number. The expressions for moments could have been derived by using this number instead. For example, the relative mean distance is

k L = Ω 1 + 2Ω (40)
For bacteria such as E.Coli present at 10 10 /ml in the human gut for example [27] and overall mutation probability at ≈ 10 -10 per base [START_REF] Wielgoss | Mutation Rate Inferred From Synonymous Substitutions in a Long-Term Evolution Experiment With Escherichia coli[END_REF], the mutation number is Θ ∼ 2L for a population contained in 1ml. In this regime, we expect to find a large distribution of distances with a peak around L/2, as was indeed supposed by Jeraldo et al. [START_REF] Jeraldo | Quantification of the relative roles of niche and neutral processes in structuring gastrointestinal microbiomes[END_REF]. On the other hand, the 1 ml volume choice is arbitrary, and if we had considered a 1µl volume instead, we would be in the low mutation regime where neutral aggregation dominates. In order to determine the effective mutation number Ω e needed to test the neutral hypothesis, we must explicitly take into account the dispersion of individuals in the real space. The next section deals with this problem.

V. SPATIALLY DISTRIBUTED POPULATIONS.

Consider a population distributed into patches, each site containing N individuals. In each site, the same WF stochastic process as discussed above takes place, with the additional condition that a progeny can descend from a parent in a neighboring site with probability m/2d, where d is the dimension of real space. To take space into account, we need to compute u(x, y; k; t), the probability of finding two individuals drawn at random in sites x and y at time t at genetic distance k.

Let us first consider a one dimensional real space d = 1. We assume translational invariance:

u(x, x ; k; t) = u(|x -x |; k; t)
and will compute u(y; k; t) where y is the absolute discrete distance between two patches. The evolution equation [START_REF] Moran | The Statistical processes of evolutionary theory[END_REF] for d = 0 of the previous section was derived by the usual combination of AND,OR, distinguishing the cases where two individuals descend from the same parent or not. For the spatial case, we must also distinguish the cases where the parent can be from the same patch or not. Again, we collect u(y; k; t) (0 ≤ k ≤ L ) into a vector |u(y; t) = u ((y; 0; t), u(y; 1; t), ..., u(y; L; t)) T Following the same line of arguments as the preceding section, the evolution equation ( 22) is generalized for the spatial case :

|u(0; t) = (1 -2m) {(1 -a)Q |u(0; t) + |b } (41) + 2mQ |u(1; t) |u(1; t) = m {(1 -a)Q |u(0; t) + |b } (42) + (1 -2m)Q |u(1; t) + mQ |u(2; t) |u(y; t) = mQ |u(y -1; t) + (1 -2m)Q |u(y; t) (43) + mQ |u(y + 1; t) (y ≥ 2)
As before, the stationary solution of these equations for large times can be obtained by using the spatial PGF function

φ(y; z) = L k=0 u(y; k)z k (44)
By its very definition, φ(y; 1) = 1 ∀y.

An evolution equation can be obtained for the PGF which has the same property as in the preceding section: (∂ n /∂ z n ) φ(y, z)| z=1 depends only on the lower derivatives, and all the factorial moments can be obtained from a hierarchical structure. For example, the mean

k(y) = L k=0 ku(y, k) = ∂φ(y, z) ∂z z=1
is shown to obey the equation

(-r + m∆) k(y) = - Lr 2 (45)
+ a {(1 -2m)δ y,0 + mδ y,1 } k(0)
where ∆ is the discrete Laplacian operator

∆f (y) = f (y -1) -2f (y) + f (y + 1)
and Equation ( 45) is a non-homogeneous two terms recurrence equation, which can easily be solved. The solution depends on the per sequence mutation number Ω = 2N ν = Θ/L and the migration number

r = 2aΘ L -2aΘ ≈ 4 λ L = 4ν
M = N m (46)
and is written

k(y) L = 1 2 - 1 2 -(1 -a) k(0 L s y (y > 0) ( 47
)
where s is the less than unity solution of

s 2 -2(1 + Ω M )s + 1 = 0 (48) and k(0)
L = Ω + M (1 -s) 1 + 2Ω + 2M (1 -s) (49)
in the regime where a 1 and m 1 (see appendix C for the full expression). Defining the effective mutation number as

Ω e = Ω + M (1 -s) (50) = Ω 2 + 2ΩM (51)
we see that the expression for k(0)/L is similar to the case of a well mixed population (equation 40) where Ω has been replaced by Ω e . Figure 6 shows the excellent agreement between the above theoretical expressions and the results from individual based numerical solution of the spatial, finite sequence length WF model. We first note that for large separations y 1 between sites, k(y) ≈ L/2. This result is expected, as far-distant sites evolve independently of each other and we expect to recover the binomial distribution for u(y, k) when y y * = -1/ log(s).

On the other hand, k(0) is given by a balance between mutation number Ω and migration number M , and their relative amplitude with respect to unity. Even when Ω 1, k(0) can be substantially greater than unity if the migration number is not too small. In particular, note that Ω/M = 2ν/m; when ν m, the effective mutation number is

Ω e = Ω 2 + 2ΩM ≈ √ 2ΩM = 2N √ νm
and the typical distance under which the mean distance notably deviates from L/2 is

y * = -1/ log s ≈ m/(4ν)
Higher moments can be found by the same method if needed, although the algebra becomes increasingly cumbersome. The extension to higher dimensions is trivial and necessitates only the redefinition of the discrete Laplacian operator ∆ (see appendix C).

We can now try to make a rough estimation of the order of magnitude of the neutral clustering for the bacterial communities studied by Jeraldo et al. Let us suppose that individuals (or their seed) diffuse in one generation according to a dispersal law which for simplicity we take to be a Gaussian

p s (x) = (2πσ 2 ) -1/2 e -x 2 /2σ 2
The dispersion length σ determines both N and m. We can model the space as divided into patches of size = 3σ. The exact size of the patch is not crucial, we have to choose it in a way that insures that (i) there is negligible migration to the next nearest neighbor and (ii) every progeny has a non-negligible probability of descending from every parent inside the same patch. = 3σ is a good compromise. For patches of this size, the migration probability between patches is m ≈ 0.27 and the number of individual inside each patch is N = c 0 , where c 0 is the bacterial concentration. Note that the migration number N m is not sensitive to the choice of .

We can, as an order of magnitude, choose 1 µm as the bacterial size; in very dense bacterial concentrations, bacterial movements are reduced and we may broadly estimate σ as some 10 times the bacterial size and therefore N = 30. With this choice of N and m, and ν = 10 -10 , Ω e ≈ 3×10 -4 : this number corresponds to low mutation regime ; k(0) /L is very small and individuals inside the same patch are at small genetic distance of each other. The number of patches over which this clustering is observable is y * ≈ 2.5 × 10 4 patches or ≈0.8m in physical dimension! If we suppose that samples were taken from few µL of bacterial intestinal residues, we see that we can expect large neutral clustering in these samples.

The above rough estimation shows that in principle, Jeraldo et al. data cannot exclude the neutral hypothesis for these communities.

VI. DISCUSSION AND CONCLUSION.

A fundamental question in population genetics is the number of segregating sites in the genome of individuals of a given population submitted to neutral mutations. We have derived in this article the probability u k that two individuals, drawn at random in a population of size N , differ at k sites of a given gene of length L. This information is usually derived in the framework of coalescent theory, which is more general that the approach developed here. For example, using coalescent theory, one derives the probability that n individuals drawn at random have k segregating sites, while the approach here is restricted to n = 2, i.e., the pair correlation function. On the other hand, application of coalescent theory to finite sequence lengths is rather difficult and it is even harder to take into account structured populations like geographically dispersed individuals (see [START_REF] Hein | Gene Genealogies, Variation and Evolution: A primer in Coalescent Theory[END_REF] p45 for a detailed discussion of these difficulties or [START_REF] Wakeley | Coalescent theory has many new branches[END_REF] for recent developments of coalescent theory). The approach we have developed here tackles these problems rather easily, making it interesting for analyzing more realistic systems.

The main assumptions of the approach we have developed is that the population is subject to neutral selection and mutation. In particular, we have supposed that all sites have the same mutation probability. This constraint can be relaxed to some extent and allow for rate heterogeneity, a subject of intense development in the field of molecular evolution [START_REF] Yang | Computational molecular evolution[END_REF][START_REF] Simon | Molecularclock methods for estimating evolutionary rates and timescales[END_REF]. Consider for example a sequence of length L = L 1 + L 2 , where sites belonging to the i-th subsequence have mutation number Θ i = L i Ω i . Then the PGF for the whole sequence is the product of the PGF for each sub-sequence

φ(z) = φ 1 (z)φ 2 (z)
and therefore all the moments of the whole sequence can be deduced from the moments of the subsequences. In particular, the mean is simply given by

k = k 1 + k 2 = L 1 Ω 1 1 + Ω 1 + L 2 Ω 2 1 + Ω 2 when Ω i
1, the system behaves as having an effective mutation rate per base equal to the weighted mean of the subsequences.

Experimental measurement of u k constitutes a stringent test of neutral theories in ecology (UNTB) introduced by Hubbell [START_REF] Stephen | The unified neutral theory of Biodiversity and Biogeography[END_REF]. It has been shown by many authors [START_REF] Rosindell | The unified neutral theory of biodiversity and biogeography at age ten[END_REF] that UNTB predicts abundance distributions that are indeed observed in nature. On the other hand, other scientists have argued that abundance distribution is not a very selective criterion and the observed abundance patterns in natural communities can be predicted just as well by other competing theories [START_REF] Linquist | Applying ecological models to communities of genetic elements: the case of neutral theory[END_REF]. The advent of modern gene sequencing tools allows one to develop more selective criteria. O'Dwyer et al [START_REF] James | Backbones of evolutionary history test biodiversity theory for microbes[END_REF] for example have observed that the reconstructed phylogeny of microbes from various habitats does not correspond to that predicted by UNTB. Jeraldo et al [START_REF] Jeraldo | Quantification of the relative roles of niche and neutral processes in structuring gastrointestinal microbiomes[END_REF] have used a variant of u k measurement to argue against UNTB for various microbial habitats. Their main argument that, as u k is a sharply decreasing function of k and is therefore in contradiction with UNTB, is however weak. The weakness of this argument was pointed out by D'Andrea and Ostling [START_REF] Rafael | Can clusterting in genotype space reveal niches ?[END_REF] who showed, by combining theoretical modeling and numerical simulation, that a neutral theory gives rise to a u k function that is sharply peaked at the origin. In this article, we have derived the explicit expression for u k and we show that the shape of the u k function depends critically on the mutation number : it can be either peaked around the origin or far from it, depending on the value of the mutation and the migration number.

For the microbial communities discussed by Jeraldo et al [START_REF] Jeraldo | Quantification of the relative roles of niche and neutral processes in structuring gastrointestinal microbiomes[END_REF], our rough estimation suggests that the UNTB hypothesis cannot be ruled out. Of course, the current state of DNA sequencing does not allow for the measurement of the spatial correlation function u(y; k) for these communities. But even if we assume that the samples used by Jeraldo et al have been mixed over few millimeters, the amplitude of neutral clustering remains important.

In this article, we have used the neutral WF model of competition where all individuals compete against all the others with equal strength. The genetic distance between individuals does not appear explicitly in this competition. Other models can be formulated where the genetic distance can modify the competition. In kin selection theory for example, competition is reduced for individuals closely related to each other [START_REF] Gardner | The genetical theory of kin selection[END_REF]. On the other hand, Biancalani et al [START_REF] Biancalani | Framework for analyzing ecological trait-based models in multidimensional niche spaces[END_REF], have investigated the problem of competitive exclusion in ecology; in their approach, individuals are represented by a sequence which determines their use of available ecological niches. In their model, competition is enhanced for individuals close to each other in terms of their resource consumption, hence giving rise to patterns of abundance in the sequence space.

The approach we have developed here shares many limitations with other analytical approaches: (i) As we do not track the sequences but only pair differences, it is not possible to compute the abundance curves as in UNTB (which, it should be stressed, was derived, only for infinite allele models). (ii) Deriving more refined results beyond the pair correlation function also seems problematic for the same reasons. The simplicity of the dynamics of pair correlation function (equations 18,21) is lost when n-correlation functions are considered. (iii) We have considered only simple mutations. In the field of molecular evolution however, the substitution rate of a base in a sequence depends on the state of the base; for example, A → T and A → C have in general different probabilities. These different rates are contained in a substitution matrix for which many models coexist [START_REF] Yang | Computational molecular evolution[END_REF]. It could be possible in principle to generalize the approach developped here to include general substitution matrices, but the algebraic cost seems at present to be prohibitive. Despite all the limitations enumerated above, we believe that the formalism developped in this article is a step forward in the search for a better understanding of natural populations. In particular, we are convinced that pair correlation function measurement as proposed by Jeraldo et al will be extended to many more natural communities and the explicit expression derived here can be used to quantify the amplitude of neutral versus non-neutral mutations.

both of them descend from parents of the same patch ( probability (1 -m) 2 ≈ 1 -2m ) or one of them descends from a parent in a neighboring patch ( probability 2 × m(1 -m) ≈ 2m ). Distinguishing these two cases leads to equation (41). Equations 42,43) are derived by following the same arguments. These equations can be written as this equation has two positive solutions 1 and 2 for which 1 2 = 1. As the mean distance is bounded by L, the solution > 1 is unphysical and we consider only the root < 1. Note also that r/m = 2Ω/M , where Ω and M are per-base mutation and migration numbers (relations 39,46). Equation (45) at y = 0 provides a linear relation between k(0) and k(1) . Using solution (C4) at y = 1 provides a second linear relation between these two quantities. Solving this 2 × 2 linear system, we find

k(0 L = Ω + M (1 -) 1 -2m + 2Ω + 2M (1 -(1 -a)) (C6)
which reduces to the expression (49) for m 1 and a 1. The coefficient C, found from the same linear system is

C = (1 -a) k(0) - L 2 (C7)
which completes the solution.

Equations for higher moments can be obtained by the same method and as in the case of d = 0, the moment of order n depends only on moments of order n -1. Their solution involves only straightforward, although cumbersome algebra. As an example, consider the second factorial moment

Figure 1 .

 1 Figure 1. Neutral WF scheme for a well mixed population.The habitat contains N haploid individuals with various alleles of a gene of length L. At generation Gt, each individual produces q progeny, some of whom may contain a mutation with respect to their parent (marked by a star). N individuals are selected among the N q progeny to constitute the Gt+1 generation.

Figure 2 .

 2 Figure 2. Two alleles of a gene of length L, differing in k of their bases (grayed area). A mutation in the grayed area will decrease the difference, while a mutation in the white area will increase it.

Figure 3 .

 3 Figure 3. Scheme of the spatial WF model for finite sequence length. Space is divided into patches of N individuals. In each patch, the WF stochastic process described in fig. 1 takes place, with the additional condition that a progeny can come, with probability m/2 from one of the neighboring patches. u(x, y; k; t) is the probability that two individuals drawn at random in patches x and y at time t are at genetic distance k.

Figure 4 .

 4 Figure 4. (Color online) The normalized factorial moments µn as a function of n for different values of Θ. L = 16, λ = 2.5 × 10 -3 and a = 1/N = 2λ/Θ. Solid lines represent theoretical values given by relation (35). Symbols are obtained by solving the linear system (24) for stationary probabilities and then computing their factorial moments from equation[START_REF] Tajima | Infinite-allele model and infinite-site model in population genetics[END_REF].

Figure 5 .

 5 Figure 5. (Color online) Individual based stochastic simulation (symbols) of the finite sequence neutral WF model and its comparison with the theoretical predictions (solid lines). Figures show the probabilities u k as a function of k for different values of N (and hence Θ) in (a) linear, (b) logarithmic representations. L = 16, λ = 2.5 × 10 -3 , N ∈ [16, 4096], Θ ∈ [0.08, 20.48]. Simulation values are obtained by generating between 10 6 (for large N ) to 10 8 (for smaller N ) paths for 10N generations (see appendix D).

Figure 6 .

 6 Figure 6. (Color online) Comparison between the theoretical result for k(y) (solid line) and stochastic simulations of the spatial, neutral, finite sequence length of WF model (symbols) for various value of (m, λ). There are Q = 128 one dimensional, circular sites with N = 32 individual per site and L = 16.

2

 2 |u(y) = (1 + m∆)Q |u(y) (C1) + a {(1 -2m)δ y,0 + maδ y,1 } (-Q |u(0) + |b )The PGF isφ(y; z) = L k=0 u(y, k)z k = η|u(y)where η| was defined in the previous appendix and is has been shown thatη|Q|u(y) = 2λ L (1 -z 2 )∂ z φ(y; z) + (1 + 2λ(z -1)) φ(y; z) = L[φ(y, z)] (C2)where the operator L[φ] captures the right hand side of equation (C2). As η| and ∆ commute, applying η| to equation (C1) leads toφ(y; z) = (1 + m∆)L[φ(y, z)] (C3) + a {(1 -2m)δ y,0 + maδ y,1 } (-L[φ(0, z)] + f (z))wheref (z) = η|b = a (1 + 2λ(z -1))Note that f (1) = a and f (1) = 2λa. As L[φ(y, 1)] = φ(y, 1) = 1, it can be checked that equation (C3) trivially verifies φ(y, 1) = 1.In order to compute k(y) , the mean genetic distance between two patches k(y) = L k=0 ku(y, k) = ∂φ(y; z) ∂z z=1 we need to apply the projection operator D = ∂ z | z=1 to equation (C3). Setting r = 4λ/L = 4ν and noting that DL[φ(y, z)] = (1 -r) k(y) + 2λ we find that k(y) must obey the relation k(y) = (1 + m∆) {(1 -r) k(y) + 2λ} + a {(1 -2m)δ y,0 + mδ y,1 } {(1 -r) k(0) -2λ(1 -a)} Grouping the terms in k(y) , dividing by (1 -r) and approximating r/(1 -r) ≈ r, (1 + a) ≈ 1 finally leads to equation (45) of the main text. The equation (45) can be solved by noting that the solution for the bulk equation (y ≥ 1) is k(y) = L

k ( 2 )

 2 (y) = L k=0 k(k -1)u(y; k) = ∂ 2 φ(y, z) ∂z 2which is obtained by applying the operatorD 2 = ∂ 2 z z=1to equation (C3):(-2r + m∆) k 2 (y) = -Lr(1 + m∆) k(y) + a {(1 -2m)δ y,0 + mδ y,1 } × k (2) (0) + Lr k(0)As k(y) is known (relation C4), this is again a solvable two term recurrence relation. The extension to higher dimensions is easily obtained by redefining the discrete Laplace operator ∆. For example, for d = 2,∆ 2 f (x, y) = -2f (x, y) 1 2 {f (x + 1, y) + f (x -1, y) + f (x, y -1) + f (x, y+1)} and the moments k(x, y) are given by k(x, y) = L 2 + C x+y

ACKNOWLEDGMENTS

We thanks M. Vallade, E. Geissler, O. Rivoire and I. Junier for the critical reading of the manuscript and fruitful discussions.

Appendix A: The infinite allele model.

The other terms of the linear system [START_REF] Houchmandzadeh | Neutral clustering in a simple experimental ecological community[END_REF][START_REF] Dumonteil | Particle clustering in Monte Carlo criticality simulations[END_REF][START_REF] Houchmandzadeh | Neutron fluctuations: The importance of being delayed[END_REF] are obtained by the same argument as the u 0 (t + 1) terms. The probability u 1 (t + 1) is obtained by considering that:

1. The two individuals are from the same parent (probability a) AND a mutation has taken place in only one of them (probability 2λ(1 -λ) ).

2. The two individuals are from different k = 0 parents (probability 1 -a) AND one mutation has taken place in one of them (probability 2λ(1 -λ) ).

3. The two individuals are from different k = 1 parents AND no mutation has taken place in any of them.

or, in other words

A1) As we will neglect terms of order λ 2 , two individuals at distance k ≥ 2 will not be from the same parent. Therefore the probability

A2) These relations are summarized in the linear system [START_REF] Houchmandzadeh | Neutral clustering in a simple experimental ecological community[END_REF][START_REF] Dumonteil | Particle clustering in Monte Carlo criticality simulations[END_REF][START_REF] Houchmandzadeh | Neutron fluctuations: The importance of being delayed[END_REF]. We derive here the stationary PGF equation ( 29) of the main text. All notations are identical to the main text (eqs. [START_REF] Moran | The Statistical processes of evolutionary theory[END_REF][START_REF] Houchmandzadeh | Alternative to the diffusion equation in population genetics[END_REF]. Consider the linear form η| = (1, z, z 2 , ..., z L ). By its very definition,

Moreover, straightforward matrix multiplication shows that

and therefore,

The stationary solution for the probabilities |u derived from equation ( 22) is given by

Applying the linear form η| to the above relation, we obtain

as we suppose that a 1 and neglect terms of O(a 2 ). Multiplying both side of the relation (B5) by L/2λ leads to equation (29) of the main text. All numerical simulations are written in C++, and data analysis is performed by the high level language Julia [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF]. For numerical simulation of well mixed populations (0 dimension), individuals are represented by the binary sequence of their genome. Using sequence lengths of powers of 2 (such as 16 and 32) allows us to represent each sequence as an integer, and to use the tools of the C language to manipulate the bits directly. If for example the powers of 2 integers are stored in an array p2, flipping the n-th bit of the integer b is performed by the operation b^p2 [n] where "^" represents the binary XOR.

For WF simulation of a population of size N , two one-dimensional integer arrays of size N are considered (present and future generations). Each element in the second array chooses randomly a parent in the first array and inherits its integer (genome), accompanied possibly by a mutation with probability λ. The process is iterated (by exchanging the role of the two arrays) T times (usually T = 10N ). Hamming distances between all elements of the final array are computed and stored in a new array H. The process is then repeated M times (M ∼ 10 6 -10 8 ) to obtain a statistically significant array H, which represent the probabilities u k .

For one dimensional spatial simulations, the same process is applied to two (present and future generation) two-dimensional integer arrays B of size Q × N , where Q is the spatial extension of the system and B[q][n] represent the genome of individual n at position q. This time, each progeny chooses its parent from a neighboring site with probability m/2, and from the same site with probability (1 -m).