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 = 500000 # !!! changed: to match Simulate.ipynb's 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: 8.86376678169621 or 230.45793632410147 years q0=5% generations to half: 8.05904782547916 or 209.53524346245817 years q0=10% generations to half: 7.7124742351991875 or 200.52433011517888 years q0=15% generations to half: 7.509741681145105 or 195.25328370977272 years q0=20% generations to half: 7.365900644919215 or 191.51341676789957 years