The effect of the dispersal kernel on isolation-by-distance in a continuous population

Simulation

In our individual-based simulation, a population exists on a 100 × 100 rectangular lattice (cf. Epperson, 1995; Epperson & Li, 1997; Epperson, 2007; Hardy & Vekemans, 1999). Individuals are uniformly spaced with a single individual per cell. Each individual contains one haploid locus. The initial population of 10,000 individuals each carry a unique allele. Generations are discrete, and individuals reproduce by generating 15 clonal offspring that experience mutations according to the infinite alleles model at rate μ. All starting and mutant alleles are selectively neutral.

The offspring disperse from the parent cell following a given dispersal distribution. The landscape boundaries are absorbing, and when offspring disperse off of the lattice they are lost. Offspring that land in the same cell will compete to become a parent in the next generation. Because all alleles are selectively neutral, a single successful offspring is uniformly selected for each cell. To avoid storing all the offspring in memory until dispersal is completed, we use a reservoir sampling method to immediately accept or reject offspring when they land on a cell (Vitter, 1985). This method allows us to keep track of two randomly chosen offspring per cell. The first offspring becomes a parent in the next generation and the second individual is recorded to measure the probability of identity-by-descent for offspring competing for the same cell. While it is possible for a cell to remain empty after dispersal, we determined that when each parent generates 15 offspring the number of empty cells per generation is negligible so we assume a constant homogeneous population density.

Modeling Dispersal

The simulation is spatially-explicit with space represented on a rectangular lattice. Due to the discrete nature of the lattice, the dispersal kernels will be discretized approximations of continuous distributions (Chesson & Lee, 2005; Chipperfield et al., 2011). The dispersal kernel function, f(r, θ; σ), takes a parameter σ and returns continuous polar coordinates. The σ parameter is the square root of one-half the second moment of dispersal distance. The polar coordinates include the angle, θ ∈ [0, 2π], which is uniformly distributed to ensure isotropic dispersal and distance, r > 0, which is drawn from a continuous distribution.

Once the angle and distance are drawn, the final position is determined by converting the polar coordinates into rectangular coordinates and adding them to the parent’s position. The new coordinates are then rounded to determine the integer coordinates of the destination cell. This dispersal scheme is similar to the centroid-to-area approximation of continuous dispersal kernels described by Chipperfield et al. (2011), which showed minimal deviation from expectations especially when cell length is less than the expected distance.

We looked at seven different dispersal distance kernels (Fig. S1): Rayleigh, exponential, half-normal, triangular, gamma, Pareto, and Lomax. We chose these distributions because they provide a range of shapes for short, intermediate, and long distance dispersal.

The Rayleigh is a distribution of Euclidean distances that result from bivariate normal displacement along the x and y axis. The Rayleigh distribution follows the assumptions of Wright (1946)’s two-dimensional isolation-by-distance model.

The exponential distribution is more leptokurtic with higher probability of dispersal at short and long distances and less at intermediate distances. The exponential tail is the boundary that separates truly heavy-tailed distributions with potentially infinite higher moments from distributions with all moments finite. The distinction is important because leptokurtic, heavy-tailed dispersal kernels are typically a better fit to observed dispersal in nature (Clark, 1998).

The half-normal distribution is equivalent to a normal distribution that has been folded over the y-axis. In this case, Euclidean distance is simply the absolute value of normally distributed random variables. The half-normal is a monotonically decreasing distribution with a convex shoulder near zero. This distribution has a higher probability of dispersal at intermediate distances compared to the exponential.

The triangular distribution is typically defined using three points: a lower limit, a, an upper limit, b, and a mode, c. Here we use a special case of the triangular distribution where a = 0 and b = c = 2σ. We chose this special case of the triangular distribution because in our dispersal function it will return polar coordinates that are uniformly sampled from within a circle with area 4πσ² which is the same as the neighborhood area (see proof in Appendix A). The triangular distribution is also the only one of our distributions that has a finite range, r ∈ [0, 2σ].

Unlike the previous single parameter distributions, the final three distributions have an additional α shape parameter. The gamma distribution is equivalent to the exponential distribution when α = 1, and as α increases the distribution becomes more symmetrical with a higher probability for intermediate distances and a lower probability for short distances.

The Lomax and Pareto distributions are both heavy-tailed power-law distributions. The n-th moments are finite only when α > n. The support for the Pareto distribution, r ∈ [x_min, + ∞), begins at a parameter x_min > 0. The Lomax distribution is a Pareto distribution that has been shifted so that the support begins at zero. We chose values of α between 2 and 3 so that the second moment of the distribution would be finite but higher moments are infinite.

The dispersal function is executed over 100-billion times per simulation, and thus it is important to make the implementation as efficient as possible. With this aim, we used an xorshift algorithm for uniform pseudo-random number generation and the ziggurat rejection sampling algorithm when applicable (Marsaglia & Tsang, 2000b; Marsaglia, 2003). We used two different versions of the ziggurat algorithm to draw distances from the exponential and half-normal distribution. For the gamma distribution we used a rejection sampling method that uses the ziggurat algorithm to draw normal variates (Marsaglia & Tsang, 2000a). Random variables from the Pareto distribution are generated by x_mine^X∕a where X is an exponentially distributed random variable that we draw using the ziggurat algorithm. The Lomax distribution is sampled similarly: x_mine^X∕a − x_min.

In addition to generating random distances, the dispersal function requires costly conversions from polar to Cartesian coordinates. We were able to avoid this conversion for the Rayleigh and triangular distributions. We simulated the Rayleigh distribution by drawing vertical and horizontal offsets from independent normal distributions using the ziggurat algorithm. For the triangular distribution, we developed a discrete sampling algorithm using the alias method that allows us to draw the vertical and horizontal offsets simultaneously in constant time (Vose, 1991). See Appendix D for a description of the algorithm, and Appendix E for an analysis of its superior efficiency compared to the other distributions we used.

Analysis

A simulation was run for each of the seven types of dispersal distributions under 4 levels of dispersal (σ = 1, 1.5, 2 and 4) with a mutation rate of μ = 10⁻⁴. Each simulation was run for a burn-in period of 10,000 generations to allow the population to reach a mutation-drift equilibrium. After the burn-in, data was collected from populations that were 1,000 generations apart to decrease the correlation between them for a total of 2,000 replicate populations per simulation. In each population, a straight transect of 50 individuals was sampled from the center of the landscape to avoid measuring edge effects.

From the transect, all possible pairs of individuals were placed into distance classes based on the geographical distance between the pair. The number of pairs that shared an identical allele was determined and recorded as a proportion of the total number of pairs in the distance class. The probabilities for each distance class were then averaged over all sampled populations. Under this sampling scheme, the number of pairs per distance class decreases as distance increases so in distance class 50 there is only one pair sampled per population.

The parameters for each dispersal distribution were calculated so that E[X²] = 2σ²; the calculations are reflected in the probability distribution functions in Fig. S1. Due to the discrete nature of the lattice, some parameters values were adjusted slightly until the simulations produced an average, observed, squared distance between parent and offspring, s², that was within 5% of the expected value, σ². Three of the distributions require a second α parameter. For the gamma distribution we used α = 1, 2, 4, 8. For the Lomax and Pareto distributions we used α = 2.4, 2.6, 2.8, 3.0 all of which are infinite in the 3rd and higher moments.

Under isolation-by-distance, individuals geographically near one another will tend to be genetically similar, and this similarity will decrease as the distance between pairs of individuals increases. Therefore, isolation-by-distance is described by constructing correlograms of genetic similarity between individuals versus the distance between them. Genetic similarity can be measured using identity-by-descent, identity-by-state, relatedness, conditional kinship, or F-coefficients and can be based on coalescent times, an ancestral population, or the current population (Hardy & Vekemans, 1999; Hardy, 2003; Malécot, 1969; Rousset, 1997; Rousset, 2002; Wang, 2014). For two-dimensional populations, genetic similarity is often plotted against the log-distance separating pairs because theory predicts that this relationship is approximately linear (Barton et al., 2013; Hardy & Vekemans, 1999; Rousset, 2000).

We recorded the probability of identity-by-descent for pairs of individuals in each distance class. Under the infinite alleles model, pairs of individuals were considered identical-by-descent if they shared the same allele. The probability of identity-by-descent in each distance class depends on the mutation rate; the probability will be greater when there are fewer alleles. For more consistent results that are nearly independent of mutation rate, the probability of identity is often calculated as a ratio that measures genetic similarity (or differentiation) relative to a particular reference group. We calculated the kinship coefficient which measures the correlation of genetic similarity between pairs of individuals a certain distance apart relative to the genetic similarity in the whole sample. (1)${\displaystyle F_r=\frac{p_{ij}-\bar{p}}{1-\bar{p}}\approx \frac{E\left[T\right]-E_{ij}\left[T\right]}{E\left[T\right]}.}$

Here p_ij is the probability of identity-by-descent between haploid individuals i and j at distance r and $\bar{p}$ is the probability of identity-by-descent between random haploid individuals in the current sample (Hardy & Vekemans, 1999). The kinship coefficient is related to differences in the expected coalescent times, T, between a specific pair of individuals and a random pair in the population (Barton et al., 2013). Kinship coefficients were calculated for each transect and then averaged across transects for each distance class. Since this statistic is highly dependent on the sampling scheme, we sampled the same transect in all simulations.

We also calculated the a_r parameter of Rousset (2000): (2)${\displaystyle a_r=\frac{p_0-p_{ij}}{1-p_0}}$ which measures genetic differentiation over distance relative to the probability of identity-by-descent within a location. The a_r parameter is independent of sampling scheme, but it does depend on the level of local identity-by-descent, p₀, in the population such that a_r approaches infinity as p₀ approaches one (Vekemans & Hardy, 2004). Typically, p₀ is estimated from the amount of autozygosity in the population; however, we estimated p₀ as the probability that an individual shared an allele with one of the offspring that it competed with for the cell, which is suitable for haploid organisms and better fits its definition (Vekemans & Hardy, 2004).

For each simulation, we calculated the average number of unique alleles in a 50-individual transect ($\bar{k}$) and the average squared distance between parents and offspring (2s²). Using $\bar{k}$, we estimated the population-level diversity, ${\hat{\theta }}_k$ (Ewens, 2004 eq. 9.32) and estimated effective haploid population size as ${\hat{N}}_e={\hat{\theta }}_k∕2\mu$ and effective density as ${\hat{D}}_e={\hat{N}}_e∕A$, where A = 10, 000.

Finally, we estimated neighborhood size using two different methods. First we used our estimated demographic parameters to calculate neighborhood size as the product ${\hat{N}}_b=4\pi s^2{\hat{D}}_e$. We calculated an estimate from samples from each population and calculated an average over all populations. We then estimated neighborhood size using the regressions of both F_r and a_r on the log of distance. The slope of the a_r regression is an estimate of 1∕2πσ²D_e and the slope of F_r regression is an estimate of −(1 − F₀)∕2πσ²D_e (Barton et al., 2013; Hardy & Vekemans, 1999; Rousset, 2000). We performed the regression for distance classes between 5 and 35. We estimated the slope from each population sample then pooled the data from all the samples to get a combined slope estimate.

Behavior of dispersal distributions

The more leptokurtic distributions (exponential, gamma-1 and Lomax) with a high probability peak near zero have a much higher probability of not dispersing from the original cell, especially when σ is low. The Pareto distribution, which has a fat tail but has been shifted so it does not have a peak at zero, has a very low probability of not dispersing. Under the gamma distribution as the α parameter increases, the probability of remaining at the origin decreases; when α = 8 the probability is nearly zero for all values of σ. Figure S2 shows the empirical cumulative distributions generated from 10,000 simulated dispersal events from each distribution. The probability of not dispersing from the original cell is indicated by the height of the left-most horizontal line for each distribution.

The average squared parent–offspring dispersal distance, s², observed for each distribution was very similar with a relative error of less than 5% from the expected σ² value (Table 1); however, the distribution of these values over sampled generations varied (Fig. S3A). Expectedly, the thin tailed or no-tail (triangular) dispersal distributions have the smallest variance because their properties are easier to represent with a small number of samples. The Lomax distribution has the highest variance with the median falling slightly below the expected value.

Figure S3B shows the distribution of the average cubed parent–offspring dispersal distances, s³, for each transect. The theoretical third moment of the Lomax and Pareto distributions is infinite, while it is not possible to simulate this on a finite landscape, we do observe values of s³ that are several orders of magnitude larger than distributions with finite third moments. The distribution of s² and s³ for the Lomax and Pareto distributions both have a large positive skew.

Allelic diversity

The distribution of the number of unique alleles is similar for most of the dispersal kernels with the median falling near the expected value under the infinite alleles model (Fig. S4). The expected number of alleles under the infinite alleles model is equal to ${\sum }_{i=0}^{n-1}\theta ∕\left(\theta +i\right)=7.03$ where n = 50 is the number of individuals in the sampled transect. The Lomax distributions have a higher median number of alleles at lower values of σ but this gets closer to the expected value when σ > 2. The average diversity is also slightly elevated for the exponential and gamma-1 simulations.

Differences in effective population size between simulations can be measured by comparing the number of unique alleles observed in the transects. Different dispersal kernels produce similar levels of diversity, except for the Lomax distributions which have a higher θ_k and consequently a larger effective population size (Table 1).

Spatial autocorrelation and isolation-by-distance

To describe the patterns of isolation-by-distance, we first measured the average probability of identity-by-descent for each sampled population as a function of distance. All dispersal kernels produced very similar patterns of isolation-by-distance especially for larger distance classes (Fig. 1). The probability of identity-by-descent is higher at small distance when σ is small but the relationship flattens out when σ = 4. Differences between the different dispersal distributions become apparent when the distance between individuals is small. The more leptokurtic dispersal distributions have a steeper incline as distance decreases and they have a higher probability of autozygosity at distance class zero. The plots for the triangular distribution nearly perfectly overlap the plots for the Rayleigh distribution in all cases.

Because the probability of identity-by-descent is sensitive to differences in the number of alleles present in the sample, we also calculated the pairwise-kinship coefficient over the log of distance (Fig. 2). The kinship coefficient shows how much more or less similar pairs of individuals in a given distance class are compared to the sample as a whole. The kinship coefficient is nearly independent of differences in allele number and there is much better overlap of the plots for the different dispersal distributions. When the kinship coefficient is plotted against the log of distance there is a negative linear relationship over a certain range of distances (Hardy & Vekemans, 1999). The slope of this linear range is also similar across distributions for each value of σ.

Finally, we plot Rousset (2000)’s a_r parameter against the log of distance. There is a positive linear relationship between a_r and the log of distance (Fig. 3). The slope of a_r is fairly similar among the dispersal distributions for a given value of σ. However, there is less overlap in the plots for the different dispersal distributions because the overall amount of differentiation varies.

Estimated neighborhood size

The ${\hat{N}}_{b\left({\theta }_k\right)}$ estimates are shown in Table 1 and Fig. 4A. Table 1 shows the average estimate over all population samples. The colored dots in Fig. 4A show this same average relative to the expected values and the bars represent the middle 50% of the individual sample estimates. As mentioned previously, the populations with Lomax dispersal tend to have a greater number of unique alleles and this translates to higher ${\hat{\theta }}_k$, higher effective population size, and ultimately higher effective density. The estimates for s² were highly variable but skewed towards lower values. As a result, the estimates of ${\hat{N}}_{b\left({\theta }_k\right)}$ for the Lomax distribution appear to be higher on average but the estimates are skewed. Otherwise, the estimates for the other dispersal distributions are similar and close to the expected values.

Table 1 shows the ${\hat{N}}_{b\left(a_r\right)}$ estimates calculated as the twice the inverse of the regression of a_r and the log of distance for the pooled sample data. Estimates using the slope of the F_r statistics were identical so they are not shown. The colored dots in Fig. 4B show the slope estimate of the combined data relative to the expected slope and the bars represent the middle 50% of the slopes from individual populations. All of the dispersal distributions have similar slopes. When σ = 4, the actual spread of the slope values is smaller than the the spread of the slopes for the other values of σ (not shown), but in Fig. 4 the values are relative so the middle 50% is wider.