Decay of linkage disequilibrium

In recombining populations, genetic linkage decays with time. This script initializes a population with high linkage disequilibrium (LD) and tracks how LD. The example can be found in decay_of_LD.py.

First, we load the FFPopSim module, along with the number crunching and plotting tools:

import numpy as np
import matplotlib.pyplot as plt
import FFPopSim as h

Next, we set up the population:

# specify parameters
N = 500000                          # Population size
L = 4                               # number of loci
mu = 0.0                            # no new mutations
r = 0.01                            # recombination rate

### set up
pop = h.haploid_lowd(L)             # produce an instance of haploid_lowd with L loci
pop.carrying_capacity = N           # set the steady-state population size
pop.set_recombination_rates(r)      # assign the recombination rates
pop.set_mutation_rates(mu)          # assign the mutation rate

# initialize the population with N/2 individuals with genotypes 0, that is ----
# and N/2 with the opposite genotype 2**L -1, that is ++++
pop.set_genotypes([0, 2**L-1],[N/2, N/2])

Third, we let the population evolve and we track linkage disequilibrium via the get_LD function:

max_gen = 50
#get LD for locus pairs (0,1), (0,2) and (0,3)
LD_trajectories = [[pop.generation,pop.get_LD(0,1), pop.get_LD(0,2), pop.get_LD(0,3)]]
for ii in range(max_gen):
    pop.evolve(5)               #evolve 5 generations
    #get a new set of LD values
    LD_trajectories.append([pop.generation, pop.get_LD(0,1), pop.get_LD(0,2), pop.get_LD(0,3)])
LD_trajectories=np.array(LD_trajectories)      #cast to an array for plotting

Fourth, we plot the resulting linkage disequilibrium curves:

cols = ['r', 'b', 'g', 'm', 'c']
for ii in range(LD_trajectories.shape[1]-1):
    #plot the LD from simulations and compare it to the exponential decay expected from theory
    plt.plot(LD_trajectories[:,0], LD_trajectories[:,ii+1], color=cols[ii], label=r'$D_{1'+str(ii+1)+'}$')
    plt.plot(LD_trajectories[:,0], 0.25*np.exp(-LD_trajectories[:,0]* r * (ii+1)), ls='--', color=cols[ii])

plt.legend()
plt.xlabel('Time [generations]')
plt.ylabel('LD $D_{ij}$')

plt.ion()
plt.show()

As expected, LD decays faster if loci are further apart. The typical plot we obtain is the following and shows complete concordance of theory and simulations: