Simulation of the fractional coalescent with two populations

The simulation code simtree.py will be available in our public GitHub. Here, a few examples of the output are shown. The key code is given in the snippet in Figure 1. Figure 2 gives examples for two populations with different $\alpha$, all histograms were drawn from 5000 independent replicates using the same effective population sizes ($\Theta_1=0.01$,$\Theta_2=0.01$) and immigration rates ($M_{2\rightarrow1}=100$, $M_{1\rightarrow2}=100$), but different $\alpha$. Each histogram is compared with the standard Kingman coalescent.

# generates Mittag-Leffler time interval based on mylambda and alpha
# for each evolutionary force: Theta_1, Theta_2, M_21, M_12
# for the future, I assume this also will work for growth and population divergence 
# with the correct lambda 
def randommlftime(mylambda, alpha):
    pia = 3.1415926 * alpha
    r1 = np.random.uniform(0,1)
    r2 = np.random.uniform(0,1)
    denoma = 1.0 / alpha
    denomlambda = 1.0 / mylambda
    return -denomlambda**denoma * 
         (np.sin(pia)/(np.tan(pia*(1.0-r1)))-np.cos(pia))**denoma * 
         np.log(r2)

# creates the time for a migration or coalescent event,
# evaluating the time intervals for each force and picks the smallest
# looping through Y , the alphas vector here needs and entry for every force
# see the function fill_Yalphas()                                                                                                 
def randomtime(Y,alphas,t0):
    smallu = 1e100
    for yi,ai in zip(Y,alphas):
        u = randommlftime(yi,ai)
        if u < smallu:
            smallu = u
    return t0 + smallu

Key python function to draw new event times using the fractional coalescent with different α per population.

Comparison of six different scenarios for two populations with different α: top left: both populations are essentially following the Kingman coalescent; top middle and right: one population deviates from the Kingman coalescent; bottom left: both populations deviate similarly from Kingman coalescent, bottom middle and right: Same scenario as ‘top middle and right’ except that the α are reversed.

usage: simtree.py [-h] [-l LOCI] [-s SITES] [-i INDIVIDUALS] [-t THETA] [-m MIG]
              [-a ALPHA] [-f FILE] [-p]

Simulate a tree

options:
-h, --help            show this help message and exit
-l LOCI, --loci LOCI  number of loci
-s SITES, --sites SITES
                    number of sites
-i INDIVIDUALS, --individuals INDIVIDUALS
                    Number of samples for each  population
-t THETA, --theta THETA
                    thetas for each population
-m MIG, --mig MIG     immigration rate matrix the diagonal values are ignored
                    so for two population specify: 0,100,100,0; for 3
                    populations stepping stone: 0,100, 0,0,0,100,0,0,0, this
                    is from 3->2->1; there need to be n*n values!
-a ALPHA, --alpha ALPHA
                    alpha for each population
-f FILE, --file FILE  treefile to be used with migdata, default is NONE which
                    is a placeholder for sys.stdout
-p, --plot            Plots density histogram of TMRCA

Example: 3 populations with stepping stone migration one-way from 3->2->1 all with alpha=0.9


python simtree.py -i 5,5,5 -l 2 -t 0.01,0.01,0.01 -m 0,1,0,0,0,1,0,0,0  -a 0.9,0.9,0.9 -f test.tre