Comparison island model formula and mean_coalescence_time()

[krissiliv]

Hi,

I was trying to compare the formula for the mean coalescence time for the island model (2*d*N + (d-1)*(d-1)/(2*m*d)) with the mean_coalescence_time() from msprime. I wanted to plot their difference for different migration rates. The formula for the island model was derived under the condition that m is small, which would mean that it has to perform better for small m. But my code gives me exactly the opposite picture:

import msprime
import matplotlib.pyplot as plt
import numpy as np
d=5
N=10000
N_T=d*N
L=set((i,j) for i in np.arange(d) for j in np.arange(d)) 
def get_mean_coaltimes(m):
    demographyIM = msprime.Demography.island_model([N]*d, migration_rate=m/(d-1))
    debugIM = demographyIM.debug()
    R=np.zeros((d,d))
    for (i,j) in L:
        R[i,j] = debugIM.mean_coalescence_time({i: 1,j: 1})
        R[i,i] = debugIM.mean_coalescence_time({i: 2})
    return [(2*d*N + (d-1)*(d-1)/(2*m*d))-int(np.mean(R))]  
M = np.linspace(0.001, 0.9, 11)
T = np.zeros((len(M), 2))
for j, m in enumerate(M):
    T[j, :] = get_mean_coaltimes(m)
plt.plot(M, T[:,0], label="mean coalescence time IM Code-Formula")
plt.legend()
plt.xlabel("migration rate")
plt.ylabel("coalecence time");

The outcome is the following plot:

I would appreciate any hints on what I did wrong
many thanks in advance

[petrelharp]

Hm, interesting. To understand the assumptions, where's that formula from?

I think there is a small bug in your code - to get the distribution of mean coalescence times for a randomly chosen pair, you should be doing

L=[(i,j) for i in np.arange(d) for j in np.arange(d)]

(i.e., without the set( )). As it is I think you are double-counting comparisons within the same location.

[krissiliv]

Thank you very much for your fast reply!
The formula can be found for example here:
https://igem.temple.edu/labs/nei/downloads/publications/1993%20Publications/1993-nei-takahata.pdf on the 2nd page
or here https://www.cambridge.org/core/services/aop-cambridge-core/content/view/FCC418CBC6F021B741C83FDE6A0E7558/S0016672300029827a.pdf/inbreeding_coefficients_and_coalescence_times.pdf on page 3
If I am not mistaken, the formula assumes diploidy, a sample of size 2, and it gives the expected coalescence time in generations.
I tried to re-run the code with the new L and it gives me this:

[petrelharp]

Do you agree that you had a bug in your code, and if so, did you fix it?

[krissiliv]

No, unfortunately the plot still looks the same and I could not find the problem. I would appreciate any help very much. Many thanks in advance!