import sys
import os
import math
import importlib
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# from Unckless et al, "Modeling the Manipulation of Natural Populations by the Mutagenic Chain Reaction"
# doi: 10.1534/genetics.115.177592
conversion_rate = 1 # 100% converted per generation
# s = fitness = 0 (neutral), so Se simplifies to just
Se = conversion_rate
Ne = 2726151 # number of people alive in the Familinx dataset in 1809. It's our population size.
q_critical = 1/(2*Ne*Se)
# how many people for gene drive to not die out?
q_critical * Ne
0.5
def calc(q0):
generations_to_half = (1/(2*Se))*math.log(Ne/q0)
# how many years is that
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2602656/
# https://en.wikipedia.org/w/index.php?title=Generation_time&oldid=859057208
generation_time = 26
generation_years = generations_to_half*generation_time
print("q0="+str(int(q0*100)) + "%", "generations to half:", generations_to_half, "or", generation_years, "years")
for ratio in [0.01, 0.05, 0.10, 0.15, 0.20]:
calc(ratio)
q0=1% generations to half: 9.711785734104552 or 252.50642908671836 years q0=5% generations to half: 8.907066777887502 or 231.58373622507506 years q0=10% generations to half: 8.56049318760753 or 222.5728228777958 years q0=15% generations to half: 8.357760633553447 or 217.30177647238963 years q0=20% generations to half: 8.213919597327557 or 213.5619095305165 years