{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import os\n",
    "import math\n",
    "import importlib\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# from Unckless et al, \"Modeling the Manipulation of Natural Populations by the Mutagenic Chain Reaction\"\n",
    "# doi: 10.1534/genetics.115.177592\n",
    "conversion_rate = 1 # 100% converted per generation\n",
    "# s = fitness = 0 (neutral), so Se simplifies to just\n",
    "Se = conversion_rate\n",
    "Ne = 500000 # !!! changed: to match Simulate.ipynb's population size\n",
    "q_critical = 1/(2*Ne*Se)\n",
    "# how many people for gene drive to not die out?\n",
    "q_critical * Ne"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc(q0):\n",
    "    generations_to_half = (1/(2*Se))*math.log(Ne/q0)\n",
    "    # how many years is that\n",
    "    # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2602656/\n",
    "    # https://en.wikipedia.org/w/index.php?title=Generation_time&oldid=859057208\n",
    "    generation_time = 26\n",
    "    generation_years = generations_to_half*generation_time\n",
    "    print(\"q0=\"+str(int(q0*100)) + \"%\", \"generations to half:\", generations_to_half, \"or\", generation_years, \"years\") "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "q0=1% generations to half: 8.86376678169621 or 230.45793632410147 years\n",
      "q0=5% generations to half: 8.05904782547916 or 209.53524346245817 years\n",
      "q0=10% generations to half: 7.7124742351991875 or 200.52433011517888 years\n",
      "q0=15% generations to half: 7.509741681145105 or 195.25328370977272 years\n",
      "q0=20% generations to half: 7.365900644919215 or 191.51341676789957 years\n"
     ]
    }
   ],
   "source": [
    "for ratio in [0.01, 0.05, 0.10, 0.15, 0.20]:\n",
    "    calc(ratio)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}