{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Fluctuation Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "In this notebook, we will simulate a virtual experiment mimicking that of Luria & Delbruck (1943) on bacterial resistance to phage treatment. In this classic experiment, they exposed bacterial cells to phages and measured statistics of the number of resistant cells, in order to test whether these resistant cells emerged in response to the treatment or preexisted before that. In the first scenario, mutations would be induced by phages and happen only after the exposure, whereas in the second scenario, mutations would happen irrespectively of the phages and some cells would have become resistant before being exposed. Their experimental results and statistical analyses confirmed the second scenario, which helped establish that genetic mutations occur randomly, as we now know since they are caused by random alterations of DNA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "Let us recapitulate the method they used to distinguish between those two scenarios. Suppose we have grown a large number $N$ of bacterial cells, and we treat them with phages to find that a small fraction $\\mu$ could survive. Therefore, on average we would find $\\mu N$ resistant cells. But if we repeat the experiment, the actual number of resistant cells would fluctuate from replicate to replicate. The key point is to measure the variance of such fluctuation and compare it to the mean ($\\mu N$). Luria & Delbruck argued that, in the first scenario, the number of resistant cells should have a Poisson distribution and hence the variance should be equal to the mean, whereas in the second scenario, the variance should be much larger than the mean. Their data showed a large variance, thus supporting the second scenario. Because the analyses relied on quantifying the fluctuation of event numbers, their method has become known as the \"fluctuation test\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "We would like to understand why the two scenarios would have different variance-to-mean ratios. For that purpose, we will simulate the growth of bacterial cells and the occurrence of mutant cells. Remember that in a typical experiment, we start from a single cell and let it grow to a large population by successive cell divisions. For simplicity, we will assume that cells divide synchronously and so the population grows geometrically. We will study the scenarios where resistant mutants appear either after or before the phage exposure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Scenario 1: resistant cells emerge after exposure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "In the first scenario, there are no mutants before the cells are treated with phages. Therefore, cells simply grow in numbers and we need to keep track of the number of cells after a given number of generations. When we treat them with phages, each cell has the same probability $\\mu$ of becoming resistant. Here is a Python class that describes such a process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class MutationAfterExposure:\n",
    "    \"\"\"\n",
    "    simulate the growth of bacterial cells in the scenario where resistant cells appear after exposure.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, mu, N0=1):\n",
    "        \"\"\"\n",
    "        initialize experiment, assuming all cells are normal in the beginning.\n",
    "        inputs:\n",
    "        mu: float, probability that a normal cell becomes resistant *after exposure*.\n",
    "        N0: int, initial number of normal cells.\n",
    "        \"\"\"\n",
    "        self.mu = mu\n",
    "        self.N = N0    # number of normal cells\n",
    "        self.M = 0     # number of mutant cells\n",
    "    \n",
    "    def grow(self, gen=1):\n",
    "        \"\"\"\n",
    "        grow cells for a given number of generations.\n",
    "        inputs:\n",
    "        gen: int, number of generations to grow.\n",
    "        \"\"\"\n",
    "        for g in range(gen):\n",
    "            self.N = self.N * 2    # cell division doubles their number in every generation\n",
    "    \n",
    "    def treat(self):\n",
    "        \"\"\"\n",
    "        treat cells, during which each cell has a fixed probability of becoming resistant.\n",
    "        \"\"\"\n",
    "        # ran = np.random.rand(self.N)    # draw a random number for each cell\n",
    "        # self.M = np.sum(ran < mu)    # a cell becomes resistant if its random number is less than mu\n",
    "        \"\"\"\n",
    "        The above method becomes computationally unfeasible when `self.N` is very large.\n",
    "        Instead, we will use an equivalent method knowing that `self.M` obeys the binomial distribution.\n",
    "        \"\"\"\n",
    "        self.M = np.random.binomial(self.N, self.mu)    # draw from a binomial distribution B(N,mu)\n",
    "        self.N = 0    # non-resistant cells die"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "Now we can simulate an experiment where we grow the cells \"overnight\", which can be about 30 generations. This would give us $N = 2^{30} \\sim 10^9$ cells. Let us assume that the mutation rate is $\\mu = 10^{-7}$. Therefore, on average we should have $\\mu N \\sim 10^2$ resistant cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "before treatment: number of normal cells = 1073741824, number of resistant cells = 0\n",
      "after treatment: number of normal cells = 0, number of resistant cells = 99\n"
     ]
    }
   ],
   "source": [
    "mae1 = MutationAfterExposure(1e-7)\n",
    "mae1.grow(gen=30)\n",
    "print(f'before treatment: number of normal cells = {mae1.N}, number of resistant cells = {mae1.M}')\n",
    "mae1.treat()\n",
    "print(f'after treatment: number of normal cells = {mae1.N}, number of resistant cells = {mae1.M}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "We see that there was no resistant cells before the treatment, consistent with the first scenario. After the treatment, no normal cells survived, and the number of resistant cells is close to what we expected. If we repeat the simulation again, we will get a slightly different number due to random fluctuations. To collect statistics on the numbers, let us repeat the \"experiment\" many times and calculate the mean and variance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "G = 30    # number of generations to grow\n",
    "mu = 1e-7    # mutation probability after exposure\n",
    "\n",
    "repeats = 1000    # number of repeat experiments\n",
    "results = []    # collect numbers of resistant cells from experiments\n",
    "\n",
    "for r in range(repeats):\n",
    "    mae1 = MutationAfterExposure(mu)\n",
    "    mae1.grow(gen=G)\n",
    "    mae1.treat()\n",
    "    results.append(mae1.M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean number of resistant cells = 107.748\n",
      "variance = 112.04449599999998\n",
      "variance-to-mean ratio = 1.039875413000705\n"
     ]
    }
   ],
   "source": [
    "mean = np.mean(results)    # mean number of resistant cells in an experiment\n",
    "var = np.var(results)      # variance of number of resistant cells among experiments\n",
    "ratio = var / mean\n",
    "print(f'mean number of resistant cells = {mean}')\n",
    "print(f'variance = {var}')\n",
    "print(f'variance-to-mean ratio = {ratio}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "We see that the variance-to-mean ratio is very close to 1, which signals a Poisson distribution. Mathematically, for $N$ trials with success probability $\\mu$ each, the number of successful trials will obey the binomial distribution (which is what we used in the code above). However, for large $N$ and $\\mu N \\gg 1$, this distribution can be well approximated by a Poisson distribution. Let us plot this distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEGCAYAAABy53LJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUZdrH8e+dRkhCbwIBQpcqQijSFQsoRRBWbICivOJiXeuyKLZdFHtbRQUBCyoqRmmCiIAgHaRDwAChlxACIZByv3/MwMaQMoFMzkxyf65rrpw585wzv1GSe55TnkdUFWOMMcZTAU4HMMYY41+scBhjjMkXKxzGGGPyxQqHMcaYfLHCYYwxJl+CnA5QGCpWrKhRUVFOxzDFxJYtWwBo2LChw0mMuTgrV648rKqVsq4vFoUjKiqKFStWOB3DFBNdu3YFYP78+Y7mMOZiicjO7NbboSpjjDH54tXCISLdRWSLiMSKyJPZvF5CRL50v75URKKyvF5TRE6IyKOe7tMYY4x3ea1wiEgg8C7QA2gM3CIijbM0GwokqGo94HXgpSyvvw7MzOc+jTHGeJE3z3G0AWJVdQeAiEwB+gAbM7XpA4x2L08F3hERUVUVkRuBHcDJfO7TGFNMpaamEh8fT0pKitNR/EpoaCiRkZEEBwd71N6bhaM6sDvT83igbU5tVDVNRBKBCiJyCngCuAZ4NLv2uewTABEZBgwDqFmz5oV/CmOM34iPj6dUqVJERUUhIk7H8QuqypEjR4iPj6d27doebePNcxzZ/V/LOqJiTm2eBV5X1RMXsE/XStVxqhqtqtGVKp13NZkxpghKSUmhQoUKVjTyQUSoUKFCvnpp3uxxxAM1Mj2PBPbm0CZeRIKAMsBRXL2I/iLyMlAWyBCRFGClB/s0xhRjVjTyL7//zbxZOJYD9UWkNrAHGAjcmqVNDDAYWAL0B+apa5z3TmcbiMho4ISqvuMuLnnt0xhjjBd5rXC4z1mMAGYDgcB4Vd0gIs8BK1Q1BvgYmCwisbh6GgMvZJ/e+gzGGGPO59U7x1V1BjAjy7qnMy2nAAPy2MfovPZpTNST0y9q+7gxNxRQEmOKPrtz3BhjCkhcXByXXnopd999N02bNuW2225j7ty5dOjQgfr167Ns2TJOnjzJXXfdRevWrbn88sv5/vvvz23bqVMnWrZsScuWLVm8eDHgGrqma9eu9O/fn0svvZTbbrsNp2duLRZjVRnjLfv3pfPZ54F8/z08ds2TlKvVkNJB4RxPO5n3xsarHnroIdasWVOg+2zRogVvvPFGrm1iY2P5+uuvGTduHK1bt+bzzz9n0aJFxMTE8O9//5vGjRtz1VVXMX78eI4dO0abNm24+uqrqVy5MnPmzCE0NJRt27Zxyy23nBtjb/Xq1WzYsIFq1arRoUMHfvvtNzp27Figny0/rHAYk08pJ1NYFfM9wfETSUxI5dH/zKFt9EnaVvqEykEH6HgXbN1fh/lvPUPHe0cRFJLzr9nFHGKzw2u+qXbt2jRr1gyAJk2a0K1bN0SEZs2aERcXR3x8PDExMbzyyiuA6xLiXbt2Ua1aNUaMGMGaNWsIDAxk69at5/bZpk0bIiMjAVfxiouLs8JhjK9ThbXzV3N85X+5rOxXtA9LZG94JCeCB7Npo3Jpo3Ay0veydeV6vn57NO1q7aNb4+d44x+1uPs/dxER4fQnKH7y6hl4S4kSJc4tBwQEnHseEBBAWloagYGBfPPNN+cNuz969GiqVKnC2rVrycjIIDQ0NNt9BgYGkpaW5uVPkTs7x2FMLjLSlbFjoX59+OTlX2lV4TPWJfRmVdm5VBm2kysfeoFLG7mugQ8IDKBBm+bM2X2U5xeVYHryLB797xA6d4a9e5w9Jm18x3XXXcfbb7997jzF6tWrAUhMTKRq1aoEBAQwefJk0tPTnYyZKyscxuQgIz2DRa/dQ+ysD4iMhNYD7ybjxgN0fGwSLa/vRmBQ7r8+N9x9HTExAaQciePA5HZsW2FXjhsYNWoUqampNG/enKZNmzJq1CgA7rvvPiZOnEi7du3YunUr4eHhDifNmTh9dr4wREdHq03kVLQV9OW46anpLH7jbjpV/4T5h5+my/3P4unNtVknctq0eC3l1vYgNDiZHdW/pWWPqwokt53jON+mTZto1KiR0zH8Unb/7URkpapGZ21rPQ5jskhPTWfJG3e5isaR0XR9wPOikZ1G7S8jvdvvHDoZSdND3Vk4aVLBhTXGAVY4jMkkPU1Z/PqddKw+iflHn6Pr/c8UyH6rN6hJ5dt+Y/3BznQKGszUV7+gGHT2TRFlV1UZk8lT/xQyNjQnvccLdB0xskD3XaZiGZqOmMG3L49hyL978eM60EqCBFoFMf7FehzGuC1fprz6KiRWe7TAi8ZZIaEh9B31NI89FcHXU05y/4mZoBleeS9jvMV6HMYAqadTSZ/bgyHd7mPs2H5efS8ReOYZiK4yixtKv0PC72F8FtbZq+9pTEGyHocxwG8fjaVd1M8MGxZA2bKF857XD+vHD5uv4+lWb1Dn5MHCeVNjCoAVDlPs1Uw+QruI51iypz9t+99YaO8rAcK/TgwlIbkcb9Z6lZB0Z+8GNhcvLi6Opk2bnrf+7rvvZuPGjV597+uvv55jx47l2iYqKorDhw9f9HvZoSpTvGkG/678PqfOhFH3b29f0C6y3ouxf8eRbNdnKySMRzY8weftH2FY/GzeCbR7M4qijz76yOvvMWNG4c02YT0OU6y1PhhPxzpLWB/yKpVrXeJIhsURDXh47vO8Pmkkp+IqOJLBFJy0tDQGDx5M8+bN6d+/P8nJyXTt2vXcSLcRERGMHDmSyy67jHbt2nHgwAEAdu7cSbdu3WjevDndunVj165dAAwZMoThw4dz5ZVXUqdOHX799VfuuusuGjVqxJAhQ869b+bexI033kirVq1o0qQJ48aNK/DPaD0OU2ylHQ/lm8/v5pq4avz0c29Hs3xTrhUHpQzMKEvU0J9JLRHoaJ6i4KGHoIBHVadFC8hr7MQtW7bw8ccf06FDB+666y7ee++9v7x+8uRJ2rVrx4svvsjjjz/Ohx9+yL/+9S9GjBjBoEGDGDx4MOPHj+eBBx5g2rRpACQkJDBv3jxiYmLo1asXv/32Gx999BGtW7dmzZo1tGjR4i/vMX78eMqXL8+pU6do3bo1N910ExUqFNyXEq/2OESku4hsEZFYEXkym9dLiMiX7teXikiUe30bEVnjfqwVkb6ZtokTkXXu12wcEXNBVJWSv5eDjAA2NiyFBFzEreEFICA4g0q9VvLjff14MOB7R7OYi1OjRg06dOgAwO23386iRYv+8npISAg9e/YEoFWrVsTFxQGwZMkSbr31VgDuuOOOv2zXq1evc0OzV6lShWbNmhEQEECTJk3ObZ/ZW2+9da5Hs3v3brZt21agn9FrPQ4RCQTeBa4B4oHlIhKjqpnPEA0FElS1nogMBF4CbgbWA9HuOcarAmtF5AdVPXv28EpVvfgzPKbYuvbEet558CZ6zpjItrKlnI4DQMglJ9iv5bin+adMXtONAyXKOB3Jrzk0qjqSZXyarM+Dg4PPrcttiPTM22Uemj3rsO1Zt58/fz5z585lyZIlhIWF0bVrV1JSUi78A2XDmz2ONkCsqu5Q1TPAFKBPljZ9gInu5alANxERVU3OVCRCAbu11hSYUqmneLHZa6zb14TttXxrBNJ3TvdBVRhRYprTUcwF2rVrF0uWLAHgiy++8HjCpfbt2zNlyhQAPvvsswueqCkxMZFy5coRFhbG5s2b+f333y9oP7nxZuGoDuzO9DzevS7bNu5CkQhUABCRtiKyAVgH3JupkCjwk4isFJFhXsxviqgngr+kXFgCT+y9n4wA37o+5FBoaT5cdzu3Nv2GqORDTscxF6BRo0ZMnDiR5s2bc/ToUYYPH+7Rdm+99RYTJkygefPmTJ48mTfffPOC3r979+6kpaXRvHlzRo0aRbt27S5oP7nx5snx7A4aZ+055NhGVZcCTUSkETBRRGaqagrQQVX3ikhlYI6IbFbVBee9uauoDAOoWbPmxXwOU4S0SdrO7R2mMnb5cLaHVXY6TrbG0Z1BKV8zNHQmoxjkdByTD1FRUdner3F2iH2AEydOnFvu378//fv3P7ftvHnzztv2k08++cv+169fn+1rmc91zJw5M9t82Z0PuRDe/LoVD9TI9DwS2JtTGxEJAsoARzM3UNVNwEmgqfv5XvfPg8B3uA6JnUdVx6lqtKpGV6pU6aI/jPF/GWcCqX0giQ37LuXDoOucjpOjpOCS3PzT+9z71nhS4ss5HceY83izcCwH6otIbREJAQYCMVnaxACD3cv9gXmqqu5tggBEpBbQEIgTkXARKeVeHw5ci+tEujF5Or60Li9/M5Lr177HmUDfvhJ9Y2Q5UoICSVpYl+Iw2ZrxL14rHO5zEiOA2cAm4CtV3SAiz4nI2YvmPwYqiEgs8Ahw9pLdjriupFqDq1dxn/sqqirAIhFZCywDpqvqLG99BlN0RKSm0PTYXsIa7kOqJTsdJ08BIelcft081jzYkfYnY52OY8xfePVrl6rOAGZkWfd0puUUYEA2200GJmezfgdwWcEnNUXdrem/8s+Rb9Pht8/Yg39c5nqsbioBgek8WWs8fQ6/COJbJ/JN8WX/Ek2RF5Z6muFNJvPD5uvYU8o/igZAWkAQL227h8uqr+O6k+ucjmPMOVY4TJF3S/pCyoUf452E/k5HybefwpuyJr4Zj9f/iKAMGz3X+AYrHKZIC007w/BGk5m15Wq2hFd1Ok7+SQBjdg2lbqU/ufq4d4flNhcvMDCQFi1a0LRpUwYMGEBycu7n09q3b19IyQqWFQ5TpNXcnYKgvH3U/3obZy0Jr0uvCV/y0acPkXHGBj/0ZSVLlmTNmjWsX7+ekJAQ3n///VzbL168uJCSFSwrHKbIykgN4Ocf+1HvldVsCI90Os4FExFWXFKd9JMlObmiRt4bGJ/QqVMnYmNdV8S99tprNG3alKZNm/JGpkG0IiIiANi3bx+dO3c+11tZuHAh6enpDBkyhKZNm9KsWTNef/11ANasWUO7du1o3rw5ffv2JSEhAYCuXbvyxBNP0KZNGxo0aMDChQu99tl8+2J2Yy5Cpe2w51QQJdrGOx3lopWofowBvT7h5atH0ztuLMdDwpyO5B/mdj1/Xc2/QYP7IC0Z5l9//ut1hrgeKYdhUZae6tXzPXrbtLQ0Zs6cSffu3Vm5ciUTJkxg6dKlqCpt27alS5cuXH755efaf/7551x33XWMHDmS9PR0kpOTWbNmDXv27Dl3p/jZ2f0GDRrE22+/TZcuXXj66ad59tlnzxWjtLQ0li1bxowZM3j22WeZO3euR3nzy3ocpkgKykgj5pqhjH9oEKGRCU7HKRAHogKoUX43w8RuXfJVp06dokWLFkRHR1OzZk2GDh3KokWL6Nu3L+Hh4URERNCvX7/zegOtW7dmwoQJjB49mnXr1lGqVCnq1KnDjh07uP/++5k1axalS5cmMTGRY8eO0aVLFwAGDx7MggX/G3GpX79+wF+Ha/cG63GYIqlr8iaql9vH7DPREOJ0moKxM6wCn64bwD3NP2XSmm5Ox/EPufUQgsJyfz20osc9jLPOnuPIzJM7/zt37syCBQuYPn06d9xxB4899hiDBg1i7dq1zJ49m3fffZevvvrq3OGqnJwdcj234doLgvU4TJE0+JIYdhyOYklYXaejFKh3z/QmQDIYFPCz01GMhzp37sy0adNITk7m5MmTfPfdd3Tq1OkvbXbu3EnlypW55557GDp0KKtWreLw4cNkZGRw00038fzzz7Nq1SrKlClDuXLlzvVYJk+efK73UZisx2GKnDrJB+nUfAlPL/0HWrJoXYV0sEQZvt3ck0GNvybxyJuUqRDhdCSTh5YtWzJkyBDatHGNx3r33Xf/5fwGuEbPHTt2LMHBwURERDBp0iT27NnDnXfeSUZGBgD/+c9/AJg4cSL33nsvycnJ1KlThwkTJhTuB8IKhymCeslSTp0J5buggp+HwBe8faw/b4x5nJsSInj0UafTmMwyD5me2SOPPMIjjzySY/vBgwczePDg815ftWrVeetatGiR7eRMmYdur1ixolfPcdihKlOkZJwO4vEP3qDbp9+SFFzS6TheEV+6LKvPNOCttyA11ek0pjiywmGKlBPrI0k7Hcqf1X1jHnFvqXrFOkZe838s/Xqq01FMMWSHqkyRoap80+E+poT3Z3bVOk7H8SqtlcQ11X8lZf9KNOMmJCC7yTSN8Q7rcZgiI/pkHFc2WEBwjaJx30auJID4iIdpfMlK1s49b+ZkY7zKCocpMu4oO5MjJ8rzU2hzp6MUiugBgzh8oiKn177mdBRTzFjhMEVCpZQkejaczadbbuJMYLDTcQpFWKmSrEu5j9ZVf+DPP7Y6HccUI3aOwxQJA2QRAZLBF+ldoHjUDQCa9P0740YfZtfykvz7TafT+J6oJ6cX6P7ixtyQr/ajR48mIiKCR3O4bnratGk0aNCAxo0bF0S8QuPVHoeIdBeRLSISKyJPZvN6CRH50v36UhGJcq9vIyJr3I+1ItLX032a4ufMGfhhXj+e+fmf7A8t63ScQlW5RmWWZrzLWx/X4Phxp9OY/Jo2bRobN/rfPCteKxwiEgi8C/QAGgO3iEjWsjoUSFDVesDrwEvu9euBaFVtAXQHPhCRIA/3aYqZb7+FRX9cxbgTvZ2O4ojh9yqXRy5k/pQ5TkcxwIsvvkjDhg25+uqr2bJlCwAffvghrVu35rLLLuOmm24iOTmZxYsXExMTw2OPPUaLFi3Yvn17tu18kTd7HG2AWFXdoapngClAnyxt+gAT3ctTgW4iIqqarKpnR+gKBc6OEubJPk0xc3DhG1xafy2hdQ45HcURrVvDuGEPcWnyP9CMvAfUM96zcuVKpkyZwurVq/n2229Zvnw54Bq1dvny5axdu5ZGjRrx8ccf0759e3r37s3YsWNZs2YNdevWzbadL/Jm4agO7M70PN69Lts27kKRCFQAEJG2IrIBWAfc637dk33i3n6YiKwQkRWHDhXPPyjFwdZl63igw8P07/EZUkxvZZAA4VC54TSovI51vyxxOk6xtnDhQvr27UtYWBilS5emd29XL3j9+vV06tSJZs2a8dlnn7Fhw4Zst/e0ndO8WTiy+zXO+nUoxzaqulRVmwCtgadEJNTDfeLefpyqRqtqdKVKlfIR2/iT/Qve5dSZUKYFtnE6iqNa9h1I4qnSHF/1X6ejFHuSzTeYIUOG8M4777Bu3TqeeeYZUlJSst3W03ZO82bhiAcyz3MZCezNqY2IBAFlgKOZG6jqJuAk0NTDfZpiIvFwIi3Lf8qKw7eQFFI0x6XyVETZCNYcu4Poyl9zZO9hp+MUW507d+a7777j1KlTJCUl8cMPPwCQlJRE1apVSU1N5bPPPjvXvlSpUiQlJZ17nlM7X+PNy3GXA/VFpDawBxgI3JqlTQwwGFgC9Afmqaq6t9mtqmkiUgtoCMQBxzzYpykm1nw7kS4RJ6l4xd/hh/1Ox3Fc1c73krTyS36eupGBD3R2Oo5PyO/lsxerZcuW3HzzzbRo0YJatWqdm3fj+eefp23bttSqVYtmzZqdKxYDBw7knnvu4a233mLq1Kk5tvM1Xisc7j/6I4DZQCAwXlU3iMhzwApVjQE+BiaLSCyunsZA9+YdgSdFJBXIAO5T1cMA2e3TW5/B+C5V+GNFAhH1utLq1lbwQ8Fer++PGrRuSpd/7GXvvmD+NgIC7PZeR4wcOZKRI0eet3748OHnrevQocNfLscdPnx4tu18jVdvAFTVGcCMLOuezrScAgzIZrvJwGRP92mKn3nz4IEPn2HixFG0KoD9FfSNYk4Z9n/B3HFHBgvmHqPrteWdjmOKKPtOYvzS1E+2U7Ei/O1v9k84s5v6KctfuILQtfc4HcUUYfZbZ/zOwZ37efvahnz8+BuEhjqdxreElhSSSnYh+pLv2f9n8bxuRNXuZcmv/P43s8Jh/M7GGZ8SFJhO02u6Ox3FJ9W+ZhhBgelsnu6bN495U2hoKEeOHLHikQ+qypEjRwjNx7cwG+TQ+BXNUCLPTGDdvnY0u/VSp+P4pFpN67Fy5jU0CB9H2pmnCAopPr/mkZGRxMfHYzf95k9oaCiRkZEety8+/6JMkbDpt+U0rrSRBSkfOB3Fp6VGDadaaj+W/jiLtv16Oh2n0AQHB1O7dm2nYxR5VjiMXzm4dDJRFUpyWe+bnY7i01r17sXgG37gSInu/NjP6TSmqLFzHMZvpKTA7a/8h1dXzaJMxTJOx/FpwSWCiOrQkxkzg/jzT6fTmKLGehzGb0ybBnsORHBFH7sr2hP33K0EbhzNpm8jqP2Px5yOY4oQ63EYv1Fyw/08fONErrrK6ST+IbKGcF2btbSOGMuZU6edjmOKECscxi/sjd1Nr4bv0vvK7TaURj4ENhxOpVKHWDHtO6ejmCLEfgWNX9g6exIBAUqdq4c4HcWvtLzhGnYerUPYnvedjmKKECscxudphhKV8Qlr9nahZuM6TsfxKwGBAewI+D9aVPuV7as3OR3HFBFWOIzPWzfvN6IqxHKi8p1OR/FLTXvfyVdLb+bLKU4nMUWFFQ7j82bMCmDWuutp0fsmp6P4pUqRlfju4BRe+bARPjqhnPEzVjiMT0tJgTEftefzPdOJKBvhdBy/NXQoVA7dzPxvVzkdxRQBdh+H8Wm/xGwhlDIMGnSJ01H82lVXKj/983oO760N/Ox0HOPnrHAYn1Z+52Msf+EPqnX9ExCn4/iU/E4+dXdKb/7V5k2ueOhT9oWWK/RpVU3R4dXCISLdgTdxTfP6kaqOyfJ6CWAS0Ao4AtysqnEicg0wBggBzgCPqeo89zbzgarAKfdurlXVg978HMYZh/ccomWVmSw68gg1gqxoXKzvacdTGW/TTxbzLjdc1KyHVnSKN4/PcYhIdRFpLyKdzz7yaB8IvAv0ABoDt4hI4yzNhgIJqloPeB14yb3+MNBLVZsBgzl/GtnbVLWF+2FFo4jaMGMKwUFpRHa6w+koRcKh0FLM2XYVt9T7noCMDKfjGD/mUeEQkZeA34B/AY+5H4/msVkbIFZVd6jqGWAK0CdLmz7ARPfyVKCbiIiqrlbVs9OXbQBC3b0TU4xUSprEpgOXUz+6qdNRiowpR66lQvhRoo4edzqK8WOeHqq6EWioqvkZ8KY6sDvT83igbU5tVDVNRBKBCrh6HGfdBKzO8t4TRCQd+AZ4QbOZ7ktEhgHDAGrWrJmP2MYXxP6xi0srr2TByVdp5HSYImRhyYZUfzSOlApK5b4rnY5j/JSnh6p2AMH53Hd2B6Wz/oHPtY2INMF1+Or/Mr1+m/sQVif3I9vjGKo6TlWjVTW6UqVK+QpunDfhy5pEPbSLRtfbTX8FKT0wkPS6JzgVWxk9YdfGmAvj6b+cZGCNiPwMnPvmr6oP5LJNPFAj0/NIYG8ObeJFJAgoAxwFEJFI4DtgkKpuz/See9w/k0Tkc1yHxCZ5+DmMH8jIgE8/haZtIqlSI+/2Jn8qttjGnJ63MP1wZyZypdNxjB/ytMcRAzwPLAZWZnrkZjlQX0Rqi0gIMNC9n6z7Hexe7g/MU1UVkbLAdOApVf3tbGMRCRKRiu7lYKAnsN7Dz2D8xNo5v/Lfm6/n/27f5XSUIim9XDqpAYHcGvU92RzlNSZPHvU4VHWi+49/A/eqLaqamsc2aSIyApiN63Lc8aq6QUSeA1aoagzwMTBZRGJx9TQGujcfAdQDRonIKPe6a4GTwGx30QgE5gIfevhZjZ84sW4iHRsuIqhnRaejFFmf7+vBu+2f4fIdu1gTXsvpOMbPeFQ4RKQrrquf4nCdl6ghIoNVdUFu26nqDGBGlnVPZ1pOAQZks90LwAs57LaVJ5mNf0o+nsxl5aayNqE/nUqHOR2nyJob2pzE5NIMiPiFNTrE6TjGz3h6qOpVXDfadVHVzsB1uO67MKZArZn+PaVLJhHR1O7d8KbTQcF8vbUXfRtMJyLNRj40+ePpyfFgVd1y9omqbnUfLjKmQAXHT2ZPeA0uu7mL01GKvE+Tr2HJxCtJCKtJcHO7j9Z4ztPCsUJEPuZ/d3DfRt4nx43Jl/374esFV9Gxyw30DrSBm73tz/CKLN7UEQlQqlrhMPngaeEYDvwdeADXOY4FwHveCmWKpy++gLE/PsqdLzudpHgQgSott9BXf+ePE2WJi6jgdCTjJzz6Wqeqp1X1NVXtp6p9VfX1fN5Fbkye/vxtOh3anqKR3SpeaMo03sXLAx/n1pJznY5i/EiuhUNEvnL/XCcif2R9FE5EUxxsW7Get/r15IW7PnI6SrFyIrgk07Zcz98axhCclusV9sack9ehqgfdP3t6O4gp3vYsmkxU+SCaXH+z01GKna9PXMXAsO+4+vQGZga1cDqO8QO59jhUdZ978T5V3Zn5Adzn/XimOEhPTadhyGesPtCdSpGVnY5T7KwMq0XsoToMrDzT6SjGT3h6cvwa4Iks63pks84UcxcyOVDrE9v5uv0e4irYrUGOkAA+39GHTqVWE3AsiIyyaU4nMj4ur3Mcw0VkHXBplvMbfwJ2jsMUiCtkM8eSy3B5z15ORym2xqX14Kp//8yRtXWdjmL8QF5XVX0O9AK+d/88+2ilqrd7OZspBjLOBPLoB2/SbuIsQsNDnY5TbAWVOkPJugcJiwshMD3d6TjGx+V1jiNRVeNwHZLSTI8IEbHZkcxFS952CZoazLGaNpWp0zp1ncWul+rS4dQ2p6MYH+fpOY7puAqGAKFAbWAL0MRLuUwx8UHTZ9hUpjEfVu94QedHTMHZVq4Ch05U5ObyP7HgzKVOxzE+zNMbAJupanP3z/q4Jk9a5N1opqireDqJG5v8SHClJCS7uSBNocoICOCL2Bvp3mAuFU8nOR3H+LALGhBIVVcBrQs4iylmrs9YQWBABjEZ7ZyOYty+Te9AYEAGfXSp01GMD/N0Po5HMj0NAFoCh7ySyBQLqsrNNWewJK41e0qWczqOcdtTshy/bu/AwKgf+GhfN8S6giYbnvY4SmV6lMB1zqOPt+xU9DsAACAASURBVEKZoq/hqQM0qbqJqfuuczqKyeKZbfdzzYtzSdllMzCa7Hk6deyzACJS2vVU7QCouSj7N9blreMj+Kn8ZU5HMVn8Wbk08SeqEbo2mJK1jjgdx/ggj3ocIhLtvhHwD2CdiKwVkTyncBWR7iKyRURiReTJbF4vISJful9fKiJR7vXXiMhK9+CKK0XkqkzbtHKvjxWRt8T60n5H04UNS67gqSX/JCmkpNNxTBYSlEHbLnP5rMd9lDpzyuk4xgd5eqhqPK7xqqJUNQrX3BwTcttARAKBd3ENTdIYuEVEGmdpNhRIUNV6uKaifcm9/jDQS1WbAYP53wRSAP8FhgH13Y/uHn4G4yPqJiTQoeZSIprudjqKyUFYg/3c2CqG61NXOR3F+CBPC0eSqi48+0RVFwF5Ha5qA8Sq6g5VPQNM4fzzIn2Aie7lqUA3ERFVXa2qe93rNwCh7t5JVaC0qi5RVQUmATd6+BmMj7ivypd8+1A/StXdl3dj44jYsCqs2HU5t9eKwfWrZsz/5DVWVUsRaQksE5EPRKSriHQRkfeA+XnsuzqQ+StlvHtdtm1UNQ1IBLJOQ3YTsNo9cVR1935y2+fZ7MNEZIWIrDh0yC4A8xURqSn0rP8T38TeQFqAp/efGid8vqcnzaptoEny3rwbm2Ilr9/cV7M8fybTcl5fQ7I795B1m1zbiEgTXIevrs3HPl0rVccB4wCio6PtK5OPuPrMH4SGnOa7lI4Q7nQak5tZIS14LiWcgeE/8zSDnI5jfEiuhUNVr7yIfccDNTI9jwSyfnU52yZeRIKAMsBRABGJBL4DBqnq9kztI/PYp/FhAy75ic0H6rMprGq23wKM70gOLsFbK+/h8M5IMuoFERBqw60bl1wLh4jcrqqfZrkB8BxVfS2XzZcD9UWkNrAHGAjcmqVNDK6T30uA/sA8VVURKYvrXpGnVPW3TO+3T0SSRKQdsBQYBLyd6yc0PiPseAb1a23nvc13ICWtbPiDt0/2Y//0TpS7egOlW8U5Hcf4iLxOjp89mFAqh0eO3OcsRgCzgU3AV6q6QUSeE5He7mYfAxVEJBZ4BDh7ye4IoB4wSkTWuB9np4YbDnwExALbAZu2zE/sWdWIyAd28Zl2dTqK8VCJS44TUf0AXdPW20lyc05eh6o+cF9We1xV8z09m6rOAGZkWfd0puUUYEA2270AvJDDPlcATfObxTgsQ0neUI3gqERSwy5oiDTjkCE3fsDb7Z+h15L/si68Rt4bmCIvz99gVU0HeufVzpjcdEjeytbnm9C8vQ2q7G/mhDbj+KlS3BIx1+koxkd4+tVvsYi8IyKdzl6i675M1xiP3FxuDiVLnOJQ5WCno5h8SgkKYcrmPvRr+COlUu1OcuN54WiPa9Km53Bdovsq8Iq3QpmipdzpZHrUn8uX23qTHhjodBxzAb48fRWhIafpnbrc6SjGB3h6B9ZQVd2ReYWI1PFCHlME9UxfRnBQGt+kd3Q6irlA28MqsySuNW3CNvCpdnY6jnGYpz2Oqdms+7ogg5iiSVW5NSqGJXFt2FUy66AAxp8MWv4yN770A6d3l3c6inFYXvdxXIrrEFUZEemX6aXSuOYeNyZXZ/aU5fE5rxDcYi9Uzru98V2n65wiIDSVlDVVnY5iHJZXj6Mh0BMoC/TK9GgJ3OPdaKYoSPqjFtPX92BF2VpORzEXKSA4gwF9JrD9wVbsjbWRjYuzvO7j+B74XkSuUNUlhZTJFBGhaWd4otEHTIroQ2pIutNxTAGIu6QUZcOPsfDHd6n20Bin4xiHeHqOo6+IlBaRYBH5WUQOi8jtXk1m/N61p/9gVJ8XadjS5nQoKvaHluHHLddxWfg4ko8nOx3HOMTTwnGtqh7HddgqHmgAPOa1VKZIGHjJTDYfqM+68GxHvjd+6pNjPSkXnsDKbz51OopxiKeF4+xdW9cDX6jqUS/lMUVEzVNHaF97KZ/H9cZm9y1aVoXVZMvBy7jk+Ftoho1fVRx5Wjh+EJHNQDTws4hUAlK8F8v4u5sCF5GWHsiPga2djmIKmIiwvdzb/O31ycz7xb4UFEceFQ5VfRK4AohW1VTgJOdPA2sMAJouBCcF8P3m6zkaEuF0HOMFV93ciT3Jl/Pmm04nMU7I6z6Oq1R1XuZ7OLIcdvjWW8GM/zq1ozJ///Z9Kt20jLB6Nm1vURQaCk+M2EXFfaPYuWEUtZrUczqSKUR5DTnSGZiH694NxTV1a+afVjjMeUr9GczR8BRK1jnsdBTjRbcPCqbcgi9YMrsctZq84XQcU4jyOlSV5J79b32mxwZgnXvZmL8of/oEWx5uxYO3jkEC7MRpUVYlqirLD/yNy8uM5/iR407HMYUor8IRgWumv1a4Zt6rClQD7gUaezea8Ue9M5YRFJjO8pI2BmZxUK7dg5QumcTqbz5xOoopRLkWDlV9VlWfBSoCLVX1UVX9B65CEpnXzkWku4hsEZFYEXkym9dLiMiX7teXikiUe30FEflFRE6IyDtZtpnv3mfWKWWNwzRDuSUqht/+bGcDGhYTjTu1Zt2+K6h1+m0y0jOcjmMKiaeX49YEzmR6fgaIym0D95Sz7wI9cPVObhGRrL2UoUCCqtYDXgdecq9PAUYBj+aw+9tUtYX7cdDDz2C8bN2832hYJZYv9vdwOoopRMeqP8G3S3sxe4ZN8lRceDofx2RgmYh8h+ukeF9gYh7btAFiz87jISJTcF3CuzFTmz7AaPfyVOAdERFVPQksEhG7VMOPJKyayPEKpZhTornTUUwhandTH255tA+NDkCPXk6nMYXB0/s4XgTuBBKAY8CdqvqfPDarDmQeQjPevS7bNqqaBiQCnhzjmOA+TDVKcrgtWUSGicgKEVlx6JBdEuptR4/CzWPGcmvM+5wOsulhi5PgYPj73xXdN5dtq2KdjmMKgaeHqlDVVar6pvux2oNNsvuDnvUyG0/aZHWbqjYDOrkfd2TXSFXHqWq0qkZXqlQpz7Dm4owfDwcSyrKqYp6nvkwR9H9DEvjh0V7sn/ey01FMIfC4cFyAeKBGpueRwN6c2ohIEFAGyHUcLFXd4/6ZBHyO65CYcVBGegbNj/XlsdtmEFIpyek4xgHlq5Zn+eHbia44mYT9R5yOY7zMm4VjOVBfRGqLSAgwEIjJ0iYGGOxe7g/MU9UcexwiEiQiFd3LwbhG67X7SRy28sdZXNt4GjfeYEWjOLukywOUDElh7bcfOR3FeJmnJ8fzTVXTRGQEMBsIBMar6gYReQ5YoaoxwMfAZBGJxdXTGHh2exGJwzVFbYiI3AhcC+wEZruLRiAwF/jQW5/BeGjruxwodQnRQ/rC2jlOpzGFIOrJ6dmunxDeibrl3qHO443JCMj+e2ncmBu8Gc0UAq8VDgBVnQHMyLLu6UzLKcCAHLaNymG3rQoqn7l4OzfsoFXVmSw4OoquoSFOxzEOG3+gDy+Vf5Uqe8+wLzLU6TjGS7x5qMoUA3/+9F8yNIBLew1zOorxAQvDGlBv9AZW/dyNnA86G39nhcNcsORkmPRja76PfZJLatssfwYICKBk9G4yDoVTdr/dSV5UWeEwF2zKFJgw729UvOoFp6MYHxLRNJ5fRnXlnfp53epl/JUVDnNBNEOJ+2UiV7RKoHNnp9MYXyKByvQjnehUZwnNT+5yOo7xAisc5oL88ctinusxhDHDp2BTipusvi3RlkNJFbmvwjdORzFeYIXDXJDTq/7DkRMVaNV/kNNRjA86HRjC+5tup3uDn2lwcp/TcUwBs8Jh8m3rsj9oU3066848SHiZcKfjGB/1ZWBHEpNL00uWOx3FFDArHCbfDv06hqSUCFr8bYTTUYwPOxEcylXffs2D77xH6uEIp+OYAmSFw+TL9th0Thw5ysrEeylbuZzTcYyPO1hPIBCSl9fIu7HxG1Y4TL6MfSWQ3q/PosHfxjgdxfiBwLAz9LtxIluHtyXyVILTcUwBscJhPLZ/52HmfL+bIUOgWvVAp+MYP7EjMoKSIad4JOJLp6OYAmKFw3hs07djWf+fBjz5sA2bbTx3tEQ476y7k35NfqTxiT1OxzEFwAqH8cixgwlEl3mP1YdupPalnkzSaMz/TAy4kgPHK/PPSyagGTaIlb+zwmE8suardykVeoIKnZ90OorxQ6eCS/DKhnvoWOd3FsascDqOuUheHVbdFA3Jx5NpGvImy/dcT+tbL3M6jvFT34a2YdHL8zme2pq1PSHI/vr4LetxmDz99NlCypZMoMTlTzkdxfix9MBAYmtHsHEjfDox2ek45iJY4TC5OnMGHvjPdfT/dBfNr+7odBzj58Ia7OeNe1/j6pRGJB+34uGvrHCYXH31eRK7d8OwB6s5HcUUASLQ9aY2RJbbxdLJbzodx1wgrxYOEekuIltEJFZEzjurKiIlRORL9+tLRSTKvb6CiPwiIidE5J0s27QSkXXubd4SsbFZvSUjPYMrjrdh8sMP06OH02lMUXHZ1R1Zuqc3LUPHcGTvYafjmAvgtdNTIhIIvAtcA8QDy0UkRlU3Zmo2FEhQ1XoiMhB4CbgZSAFGAU3dj8z+CwwDfsc1n3l3YKa3PkdxFPXkdACuPLGRCe03M2brIP711HSHU5mipEK3/xCxsRmLprxAl0fecDqOySdv9jjaALGqukNVzwBTgD5Z2vQBJrqXpwLdRERU9aSqLsJVQM4RkapAaVVdoqoKTAJu9OJnKLZUlb9HfkHsoTr8HJa1dhtzceq1bMzi/XdyeZnx/Ln1uNNxTD55s3BUB3Zneh7vXpdtG1VNAxKB3O4uq+7eT277BEBEhonIChFZcejQoXxGN21O7iC65mrejb2NjAA7FWYKXv3+L9L6mXWMHF3a6Sgmn7z5FyG7cw9Zbxn1pM0FtVfVcaoararRlSpVymWXJjuDy81g77FL+DE02ukopoi6JKoKA4bU4osvYNWyE07HMfngzVtw4oHMYylHAntzaBMvIkFAGeBoHvuMzGOf5iKlxJfj9u+m0LrnLFJr22CGxnsefxxaJN8GC/aj0XORALvWxR94s3AsB+qLSG1gDzAQuDVLmxhgMLAE6A/Mc5+7yJaq7hORJBFpBywFBgFveyN8caUZSuL8BqQGhLKlelkCSHc6kinCSpeGSpe2pWXEg6z4cTbRvbvn2v7shRsXIm7MDRe8rfkrrx2qcp+zGAHMBjYBX6nqBhF5TkR6u5t9DFQQkVjgEeDcJbsiEge8BgwRkXgRaex+aTjwERALbMeuqCpQv3/9NUtHXEWj6xYQEGJFw3jfFYPuZefROpTe8ThpZ9KcjmM84NXRYlR1Bq5LZjOvezrTcgowIIdto3JYv4LzL9E1BeDMqdNUP/QUx4NCOVEnBbDDVMb7QkJD2Fv5Za5I68/88W/R9d5HnI5k8mCXy5hzlkx+n5rld/BC3DA0wIqGKTzt+vdj2Z6eRJ76L3E7Up2OY/Jg41MaABIPHaOZPMfKPVezJLxutpevGeMtEiBU7/chbdqVpPlPwcyY4RqexPgm63EYABZPep+yJROI6DQWG8XFOKF63Ut4YlQZ5s5JZcYXG5yOY3JhhcOwcycMePofjF05h4btWjgdxxRjf/87THvybtokXUXC/tyuzDdOssJhePqfSaRnBHPrI92cjmKKucBAqNfrEcqFHWH9p485HcfkwApHMff7198wtmNd3nh2IzVq5N3eGG9r2PYyFh15lE7VxrNqxs9OxzHZsMJRjB3dd5i6R+/jSEoN7nqwvtNxjDmnzV1Ps+NwQ6rvup1Duw44HcdkYYWjGNv06YOUKZlAYIcJBJcIdjqOMeeElQ4jrd1XxB+J5B8PJpFu96L6FCscxdTSqdPoUP1zFif+iwZtmjsdx5jzNGjTnLVVlzF5Wj1eeMHpNCYzu4+jGDp6FDbPn0WZui3ocN9TTscxJkd33iksWXSSqP33sXLGIKfjGDfrcRRDDz0Ed3/wX1I7/2KHqIxPE4E3XocOl66g5u5bKZ9iw6/7AutxFDO/fbeAxT9VY+TIejRrVdbpOKYYupARbmslP86My+/jjYpvMeT4kza5mMPsv34xcuxgAnX3D+Sbf9zBP5/Kbb4sY3zLzrCKPLryKTrXWcx9aTYgttOscBQj6yc+TMWIg4R0eI+QEjasiPEvMyMu55O1A7mz4ZeE7g5xOk6xZoWjmFg+bTodq09kUcJTNGp/udNxjLkgL6bfTNvXfmXbd11IOx7qdJxiywpHMZB46Bg19g1j28GmXHH3v5yOY8wFSw0K5lSXw2i6cG38VgLtBg9HWOEoBp78ZzBf/d6f1OgJlAgr4XQcYy5KcPmT9LrjE8YPvI8HM35wOk6x5NWrqkSkO/AmrqnkPlLVMVleLwFMAloBR4CbVTXO/dpTwFAgHXhAVWe718cBSe71aaoa7c3P4K/OXrlyantFDk5tS+l29/Paj1vgxwufs9kYX7G2YjXGrb6dB1p9xNrF9fg5wiYFLUxe63GISCDwLtADaAzckmne8LOGAgmqWg94HXjJvW1jYCDQBOgOvOfe31lXqmoLKxq5q3YqgRmXD6dFs98p22Gb03GMKVBjGcCSuDa81/ppWp/Y7nScYsWbh6raALGqukNVzwBTgD5Z2vQBJrqXpwLdxDWLUB9giqqeVtU/gVj3/oyHwlJP83GtF4iquJOSV8YiQRlORzKmQKUGBnLPkcfYdqgeHzQfTfABuy2tsHizcFQHdmd6Hu9el20bVU0DEoEKeWyrwE8islJEhuX05iIyTERWiMiKQ4cOXdQH8Tfpqem8GvFf6lXczrA/nmNvyXJORzLGK04El2TQnqcZ+PFkdnx5JamHI5yOVCx4s3Bkd6NA1rvOcmqT27YdVLUlrkNgfxeRztm9uaqOU9VoVY2uVKmSp5mLhIXv/pMeDefy1IonWBkR5XQcY7zqaIkINrQIQwIyuGLrXqqkHHc6UpHnzcIRD2SeGigS2JtTGxEJAsoAR3PbVlXP/jwIfIcdwvqLTyedQQ8v48M1tzM17Aqn4xhTKILLJVPv1p/54PYRfF57JJGnbNpZb/Jm4VgO1BeR2iISgutkd0yWNjHAYPdyf2Ceqqp7/UARKSEitYH6wDIRCReRUgAiEg5cC6z34mfwKzHfZzD0nhD+veQnxugAp+MYU6hSysMd6/5N6dAkpjV+mGYnd+e9kbkgXisc7nMWI4DZwCbgK1XdICLPiUhvd7OPgQoiEgs8Ajzp3nYD8BWwEZgF/F1V04EqwCIRWQssA6ar6ixvfQZ/8tvnn1NxdSfat07iq6nBpAcG5r2RMUXMxvBI+mx8naOnyvJ1q4e4+sQ6pyMVSV69DEFVZwAzsqx7OtNyCpDtV2NVfRF4Mcu6HcBlBZ/Uvy0YP46OIfeyNqILMT9AKTsXboqxvSXLcdPuMbx/8lUu2ZNG4rE6lG6zw+lYRYpdv+bnfn3/FbqUfozle66n6fCplIwo6XQkYxyXFFySwceeYt+WliRviqROwF7SzqQRFGJ/8gqCDTnip1Rh5htv0qX0YyzZM4DL7v/OioYxmaQFBlKx11pqdV3OnJsHsObNXiQlJDkdq0iw8uujcpvsRhUS5jWiws7ePHjLJt4scT06ek4hpjPGP4gAbQ/y3Lr7eaXNi2yf2InSvadTtU7WW8pMfliPw89IRjrX7d7KiZVRHKkJr4f2QgPsRLgxuZkW3oY1ZaZTrfR2mNWWLUvXOh3Jr1nh8COB6em8HDyecT0f4aabP6Jct42ub1TGmDxF976OvU0WAbDw4/d4803IsJF4LogVDj8RkZrCO2FvM6DJ94xe+gjLa1W3omFMPjVsexkB1y/np0Mv89BDMKTfJnZv+tPpWH7HCocfaJK0h+lRD3NNvV949PeRfFLyKqcjGeO3qkRV5ctvyzB+PNx7+TDKLmnOgvEfkJFu3Q9PWeHwYZouHFvQgKNzmnIqNZRey96zYUSMKQAicOedUOuWz9h2tC2dQ+9l61vRLI+ZiWZkHVLPZGWFw0fVSj5M37j1JC6pz7oSNeix+3U2hduVIMYUpOoNatLi4Z/4LWMi4SEJtD5xPc/c9TULFzqdzLfZ5bg+Jn5zHHE/jGZem8kkJpdhlrQipU7KecMKG2MKRkBgAB1uH8SZlIH88ulkJv58I89PhJF3zeK2uyrRqEMrpyP6HOtx+IiDuw/z66sPUHlZA6IrT+G91XfSdfMHpNRJcTqaMcVCSGgIV949lE1bQnj5ZeXmhk/SaGc0S8YOYPvqTU7H8ylWOBx27BiMHAlXtE2lRZmJ/H5wCAkdY3k1uC/HQ8KcjmdMsRMWBo89JtS8awHzjzxDk/KzidrQlEVjB7Px961Ox/MJVjgcknw8mfnvjWHJS73597+VNl2qcrjDLjo/Oo6qdSOdjmdMsVemYmm63j+a1B47WHjoEVpV+opH742lVSuY8MEREg8nOh3RMXaOo5DFb44jdu5kGgW9R9ey+1l28gbWrjhB81alcM1jZYzxJRWqVaTrw2Np+o+2HKhSncQ1iez/5RWCQ97k69hrmJfYmsUhDUnM5QhB3JgbCjGx91nhKARJSfDNN7Dp55m81ON6IsvDyj1Xc6De17S5taPT8YwxHjgRXJLwy+MJaxHPvJP1qbypLz3rzGVAkxgyMoR527rQf95/KVnnMCGXJBbpG3StcOQht8EGcxKUkUbj5L080vAgZU/P44sFvXj5h0dp3rgT3du+QP1rbqfVrbW8kNYY420isC2iCs9zK8//OZBLkw/QOfAPAo4Hk7ioIYmLLmXOv65iX2pF5p9oxfzgJhf1fhfyN+gsb/V0rHAUAM2AMwfKkLKrPJM6/IOr6i4gvEQyGRnC5oMtad2+LIuegPbtIxAZ6XRcY4qli/kDnCMJYHN4VTZTFUIh8v65pMWVJSGjND3rzuHGgJlUf3ELvF7wb+0kKxz5FJyWSu0zR6ibsY8GIbtpVGoHkip0n/QTAAcaX8LkDQNYktKUN58eSuPK5WjscGZjTOEIDDtDYOODPJZ+DwHbhlI9MYmwK3YC1ZyOVqC8WjhEpDvwJhAIfKSqY7K8XgKYBLQCjgA3q2qc+7WngKFAOvCAqs72ZJ/edGJ9dd6pPYaBV3wJwJm0YLYcbMCyAy2o2GsVoTWP8kJEf1fjCGjx2uLCimaM8TEZAQHsLleG0HJHnY5S4LxWOEQkEHgXuAaIB5aLSIyqbszUbCiQoKr1RGQg8BJws4g0BgYCTXCV6rki0sC9TV779JqQKolM2teLOYsvJzagKvElypIeGAiBEN54X2FEMMb4Ia8cJnOQN3scbYBYVd0BICJTgD5A5j/yfYDR7uWpwDsiIu71U1T1NPCniMS694cH+/SakEon2FipLBspWxhvZ4wxPsmbhaM6sDvT83igbU5tVDVNRBKBCu71v2fZ9uwIf3ntEwARGQYMcz89ISJbLuAzFKSKwGGHM+SXv2X2qbw7X+rpSTOfyuwhy+x9BZJXXrroHNle/unNwpHdVcxZx+rLqU1O67O70z3b8f9UdRwwLreAhUlEVqhqtNM58sPfMvtbXrDMhcXfMvt6Xm8OORIP1Mj0PBLYm1MbEQnCdev00Vy29WSfxhhjvMibhWM5UF9EaotICK6T3TFZ2sQAg93L/YF5qqru9QNFpISI1AbqA8s83Kcxxhgv8tqhKvc5ixHAbFyXzo5X1Q0i8hywQlVjgI+Bye6T30dxFQLc7b7CddI7Dfi7qqYDZLdPb32GAuYzh83ywd8y+1tesMyFxd8y+3RecX3BN8YYYzxjw6obY4zJFyscxhhj8sUKhxeIyMMiskFE1ovIFyIS6j6hv1REtonIl+6T+z5DRB50590gIg+515UXkTnuzHNEpJzDGceLyEERWZ9pXbYZxeUtEYkVkT9EpKUPZR7g/u+cISLRWdo/5c68RUSu85G8Y0Vks/u/43ciUjbTa47mzSXz8+68a0TkJxGp5l7vs/8uMr32qIioiFR0P/eJzH+hqvYowAeuGxX/BEq6n38FDHH/HOhe9z4w3OmsmTI3BdYDYbgumJiL60q2l4En3W2eBF5yOGdnoCWwPtO6bDMC1wMzcd0T1A5Y6kOZGwENgflAdKb1jYG1QAmgNrAdCPSBvNcCQe7llzL9N3Y8by6ZS2dafgB439f/XbjX18B18c9OoKIvZc78sB6HdwQBJd33poQB+4CrcA2rAjARuNGhbNlpBPyuqsmqmgb8CvTFNZzLRHcbxzOr6gJcV99lllPGPsAkdfkdKCsiVQsn6f9kl1lVN6lqdiMZnBtqR1X/BDIPtVMocsj7k/vfBbhGdDg7t7Hjed35sst8PNPTcP53o7DP/rtwex14nL/e2OwTmTOzwlHAVHUP8AqwC1fBSARWAscy/fJlHkLFF6wHOotIBREJw/UNpwZQRVX3Abh/VnYwY05yypjdkDe+9N88O/6Q+S5c337Bx/OKyIsishu4DXjavdpnM4tIb2CPqq7N8pLPZbbCUcDcx9j74Oq6V8P1badHNk195jpoVd2E6xDEHGAWrsMPablu5Ps8GfLG1/h0ZnHNQpYGfHZ2VTbNfCavqo5U1Rq48o5wr/bJzO4vbCP5X4H7y8vZrHM0sxWOgnc18KeqHlLVVOBboD2u7uXZGy59bqgUVf1YVVuqamdcXehtwIGzXWL3z4NOZsxBThn9cXgan80sIoOBnsBt6j7wjg/nzeJz4Cb3sq9mrovry+ZaEYnDlWuViFyCD2a2wlHwdgHtRCRMRATohusO+F9wDasCrmFWvncoX7ZEpLL7Z02gH/AFfx0Sxucyu+WUMQYY5L4ipR2QePaQlg/LaagdR4lr8rQngN6qmpzpJZ/MCyAi9TM97Q1sdi/75L8LVV2nqpVVNUpVo3AVi5aquh9fzOz02fmi+ACexfUPdT0wGddVJ3Vw/VLFAl8DJZzOmSXzQlwFbi3Qzb2uAvAzrt7Hz0B5hzN+geu8zBUTiQAABQRJREFUUSquX6yhOWXE1b1/F9eVPuvIdPWSD2Tu614+DRyA/2/v7EKsLMI4/vsrgqKpiHYjhET05UWKUQgaBGIkRIHSXmSwRBdaGSUERR+UQRIZSUmYSi2URIYVbV24RlqBheJHu1ookQtdFoipRaQ+Xcwcz7tnzx7fd0k7J/4/GM6cd2aeec6cj2femTPPw45C/aezzkeBO9tE359Ia+yHctrYLvq20Hl7/v71A73AzHb/XDSUD1L/V1Vb6FxMdjlijDGmEl6qMsYYUwkbDmOMMZWw4TDGGFMJGw5jjDGVsOEwxhhTCRsO879H0u5GL7SXqJ9HJf0oaevFa1eSu0XSjS3Ku2veX0cpf46kJaNt3yCrR9KynL8s424uPzYcxrSgcNq/DA8BSyLivn9JHgAR8WBE/NCiSjfJvc1omUPyT2ZMKWw4TFsgaVaerW/OsSr6JE3IZRdmrpKmZ5cMtZn2J5J6JR2X9Iik1ZIOSvpO0rRCF8sl7VGKOXJLbj8xx0XYl9vcXZD7oaReoK+JrquznMOqxy7ZSDrk+amkxxvqD5Mn6Yncb7+kFwr6fC7p+yy7q/j6JY3NM/rDkgaU4r4sA24GtirFnpgg6bks+7CkTdmDQU3Oy5L2SjomaaFSXJg1QFdu39Wg+1hJ63J//ZJW5evzJH0lab+kHWrhrbWZ3hf/RJi25r8+gejkFBEAs0gO9Obk59uA5Tm/m3xaFpgODOZ8N+lU8xXADJIn4hW57DXgsUL7zTl/GzkGAvBSoY+pwDGSU8pu0mneYSflgXmk07sTgUnAEWBuLhskn/ZtaDNEHim+xSbSieAxwGdZr6U1PXO9KcXXn/veWSif2jg++fm0Qv5d4K5CvVdzfgnwRUG/DSO8LytJp7Br8TimAeOAPcCMfK0LeDvne4BlZfR26tzkOw7TThyPiEM5v59kTC7Grog4FRG/kgxHb74+0ND+fbgQB2GyUhS7xcCTkg6RfuTGA1fl+jsjolm8hAXAxxFxJiJOk5xYLiyhZ1He4pwOAgeA60l+ngaARfmuYGFEnGyQ8TNwtaQ3sv+o32nO7UrRJgdIcWBmF8o+yo9lx3cRycXIWYD8Gq4jBf/amcfuGeoxOppRVm/TIVRebzXmEvJXIX8OmJDzZ6kvq45v0eZ84fl5hn6+G33rBGnGvzQagipJuhU4M4KOzVxcl6EoT8DaiHhrmHBpHuluYK2kvohYc0HhiBOSbgLuAB4G7iXFxyi2Hw+8SboD+UXS8wwds9r4nKPc918MHzsBRyJifon2pfQ2nYXvOEwnMEha7oC6h+Gq1PYLFpC8i54khehcVdgDmFtCztfAPUrejyeSHBZ+U1GXHcADkiblfmdKujL/M+qPiHiPFAxsSGxppRjUYyJiO/BsofwUabkO6kbityy/zHgV2zfSB6yobernfaOjwAxJ8/O1cZJmj9C+ld6mQ/Edh+kE1gHbJN0PfDlKGSck7QEmU5/tvgisB/qz8RgkxZwYkYg4IKmHuvvwLRFxsIoiEdEn6Qbg22yzTgPLgWuAVySdJ3lNXdnQdCbwjqTahO+p/NgDbJT0JzAf2Exa9hoE9pVQaRf1Jbu1EfFBoWwLcC1pjP4m7cFsyJvyr0uaQvodWU/a72nGSHqbDsXecY0xxlTCS1XGGGMqYcNhjDGmEjYcxhhjKmHDYYwxphI2HMYYYyphw2GMMaYSNhzGGGMq8Q/q4ZHgPfiNtAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 2**G    # number of all cells in an experiment\n",
    "M = N * mu    # expected number of resistant cells\n",
    "x_array = np.linspace(M-30, M+30, 41).astype(int)\n",
    "b_array = st.binom.pmf(x_array, N, mu)    # calculate the binomial probability mass function\n",
    "p_array = st.poisson.pmf(x_array, M)    # calculate Poisson probability mass function with same mean\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(results, bins=20, density=True, label='data')\n",
    "plt.axvline(mean, color='k', label='mean')\n",
    "plt.plot(x_array, b_array, color='b', label='binomial')\n",
    "plt.plot(x_array, p_array, color='orange', linestyle='--', label='Poisson')\n",
    "plt.xlabel('number of resistant cells')\n",
    "plt.ylabel('distribution')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "Thus, in the first scenario where resistant cells emerge after exposure, the number of resistant cells from repeat experiments would show a Poisson distribution, having a variance-to-mean ratio close to 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Scenario 2: resistant cells appear before exposure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "Now consider the second scenario, where resistant mutants arise spontaneously regardless of the phage. This means that resistant cells can appear anytime during the growth of the bacterial population, even before encountering phages. Let us denote by $\\nu$ the probability that a normal cell mutates and becomes resistant during its cell cycle (which can be different from $\\mu$ above). Note that, whenever a cell becomes resistant, *all* its descendants will also be resistant. We assume that the resistant cells grow at the same rate as the normal cells, so the population grows geometrically regardless of its composition. To simulate the growth process in this scenario, we need to keep track of the numbers of both normal and resistant cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "class MutationBeforeExposure:\n",
    "    \"\"\"\n",
    "    simulate the growth of bacterial cells in the scenario where resistant cells can appear anytime before exposure.\n",
    "    \"\"\"\n",
    "    \n",
    "    def __init__(self, nu, N0=1):\n",
    "        \"\"\"\n",
    "        initialize experiment, assuming all cells are normal in the beginning.\n",
    "        inputs:\n",
    "        nu: float, probability that a normal cell becomes resistant *during its generation*.\n",
    "        N0: int, initial number of normal cells.\n",
    "        \"\"\"\n",
    "        self.nu = nu\n",
    "        self.N = N0    # number of normal cells\n",
    "        self.M = 0     # number of mutant cells\n",
    "    \n",
    "    def grow(self, gen=1):\n",
    "        \"\"\"\n",
    "        grow cells for a given number of generations.\n",
    "        inputs:\n",
    "        gen: int, number of generations to grow.\n",
    "        \"\"\"\n",
    "        for g in range(gen):\n",
    "            self.N = self.N * 2    # normal cells divide, some may mutate afterwards\n",
    "            self.M = self.M * 2    # mutant cells divide and all their descendants are resistant\n",
    "            mutate = np.random.binomial(self.N, self.nu)    # draw a random number of normal cells that mutate\n",
    "            self.N = self.N - mutate    # update the number of normal cells\n",
    "            self.M = self.M + mutate    # update the number of mutant cells\n",
    "    \n",
    "    def treat(self):\n",
    "        \"\"\"\n",
    "        treat cells, which kills normal cells but does not affect resistant cells.\n",
    "        \"\"\"\n",
    "        self.N = 0    # normal cells are killed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "Let us simulate experiments in this scenario, growing cells \"overnight\" for $G = 30$ generations to have $N = 2^G \\sim 10^9$ cells as before. It can be shown that, in the second scenario, the expected number of resistant cells will be $\\approx G \\nu N$. Therefore, to match the expected number $\\mu N$ of resistant cells in the previous simulation, we will use a smaller mutation rate, $\\nu = \\mu / G \\approx 0.33 \\times 10^{-8}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "before treatment: number of normal cells = 1073741807, number of resistant cells = 17\n",
      "after treatment: number of normal cells = 0, number of resistant cells = 17\n"
     ]
    }
   ],
   "source": [
    "mbe1 = MutationBeforeExposure(0.33e-8)\n",
    "mbe1.grow(gen=30)\n",
    "print(f'before treatment: number of normal cells = {mbe1.N}, number of resistant cells = {mbe1.M}')\n",
    "mbe1.treat()\n",
    "print(f'after treatment: number of normal cells = {mbe1.N}, number of resistant cells = {mbe1.M}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "In this case, there were already resistant cells even before the treatment, and only those survived after the treatment. However, the number of resistant cells is not so close to what we expected, which should be around $G \\nu N \\sim 10^2$. In fact, if we run the simulation multiple times, we will find that the number of resistant cells fluctuates substantially. To have better statistics on the numbers, we need to repeat the \"experiment\" many more times (which is impossible for real experiments)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "G = 30    # number of generations to grow\n",
    "nu = 0.33e-8    # mutation probability per normal cell per generation\n",
    "\n",
    "repeats = 100000    # number of repeat experiments\n",
    "results = []    # collect numbers of resistant cells from experiments\n",
    "\n",
    "for r in range(repeats):\n",
    "    mbe1 = MutationBeforeExposure(nu)\n",
    "    mbe1.grow(gen=G)\n",
    "    mbe1.treat()\n",
    "    results.append(mbe1.M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean number of resistant cells = 64.65265\n",
      "variance = 2061659.978777977\n",
      "variance-to-mean ratio = 31888.251738760548\n"
     ]
    }
   ],
   "source": [
    "mean = np.mean(results)    # mean number of resistant cells in an experiment\n",
    "var = np.var(results)      # variance of number of resistant cells among experiments\n",
    "ratio = var / mean\n",
    "print(f'mean number of resistant cells = {mean}')\n",
    "print(f'variance = {var}')\n",
    "print(f'variance-to-mean ratio = {ratio}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "You may be amazed by the incredibly large variance-to-mean ratio in this case. Because of the large fluctuations, our calculated mean number of resistant cells is not that accurate. Let us visualize the large fluctuations by plotting the histogram of the simulation results (in log scale!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAXf0lEQVR4nO3df5BX9X3v8ecbgtmoaBQ0V1kjJBAUJQGyoBFxrGEMmqKNOlGqt6EavJJL7rWdZC4d6423TmLS2KkV7ChGgm2N1vgTGjNq2lqTkFsBQQW9KLE07phGxIRSDQ3o5/7xPeCXZffw3WXP+X7P7vMx8539fs/3nM/3/VnYfe05n3M+J1JKSJLUkyHNLkCS1NoMCklSLoNCkpTLoJAk5TIoJEm53tPsAg7EyJEj0+jRo5tdhlrIxo0bARg/fnyTK5Fa15o1a15PKR3V6PqVDIqImA3MHjt2LKtXr252OWohZ555JgBPPPFEU+uQWllE/Gtv1q/koaeU0oqU0pWHH354s0uRpAGvkkEhSSqPQSFJylXJMQpJatTOnTvp7Oxkx44dzS6ldG1tbbS3tzNs2LADasegkDSgdXZ2Mnz4cEaPHk1ENLuc0qSU2Lp1K52dnYwZM+aA2vLQk6QBbceOHYwYMWJQhQRARDBixIh+2ZOqZFBExOyIWLJt27ZmlyKpAgZbSOzWX/2uZFB4eqwklccxChVi9MLv9Wr9zV//dEGVSDpQldyjkCSVx6CQpAJt3ryZE044gc9//vOcfPLJXHrppfzgBz9g+vTpjBs3jqeeeoo333yTyy+/nKlTpzJ58mQefvjhPdvOmDGDKVOmMGXKFFauXAnUpqg588wzueiiizjhhBO49NJLKfJupR56kjRoXH311axbt65f25w0aRI33XRT7jqbNm3iu9/9LkuWLGHq1Kl85zvf4Uc/+hHLly/na1/7GhMmTOCss85i6dKl/OpXv2LatGnMnDmTo48+mscff5y2tjZeeukl5syZs2d+u7Vr17JhwwaOPfZYpk+fzo9//GNOP/30fu3bbgaFJBVszJgxTJw4EYCTTjqJT37yk0QEEydOZPPmzXR2drJ8+XJuvPFGoHZK789+9jOOPfZYFixYwLp16xg6dCgvvvjinjanTZtGe3s7UAurzZs3GxQqnwPSGmj295d/Ud773vfueT5kyJA9r4cMGcKuXbsYOnQo999//z7T41933XV84AMf4JlnnuGdd96hra2t2zaHDh3Krl27Cqu/kkFRP824Bp7eBJThpIHgU5/6FIsWLWLRokVEBGvXrmXy5Mls27aN9vZ2hgwZwp133snbb7/dlPoqOZjtdRSSBpJrr72WnTt38tGPfpSTTz6Za6+9FoAvfOEL3HnnnZx66qm8+OKLHHLIIU2pr5J7FJJUFaNHj2b9+vV7Xi9btqzb92677bZ9th03bhzPPvvsntc33HADULtB1+6bdAEsXry4n6veWyX3KCRJ5TEoJEm5DApJUi6DQpKUy6CQJOUyKCRJuTw9VtKg0tsZB/antxd9XnfddRx66KF86Utf6vb9hx56iI985CNMmDChP8rrF+5RSFILeeihh3j++eebXcZeDApJKthXv/pVxo8fz8yZM9m4cSMAt99+O1OnTuVjH/sYF154IW+99RYrV65k+fLlfPnLX2bSpEn89Kc/7Xa9srVUUETEIRGxJiJ+u9m1SFJ/WLNmDffccw9r167lgQceYNWqVQBccMEFrFq1imeeeYYTTzyRO+64g9NOO43zzjuPb37zm6xbt44Pf/jD3a5XtkKDIiKWRsRrEbG+y/JZEbExIjZFxMK6t/4XcG+RNUlSmX74wx/ymc98hoMPPpjDDjuM8847D4D169czY8YMJk6cyF133cWGDRu63b7R9YpU9B7FMmBW/YKIGArcApwDTADmRMSEiJgJPA/8ouCaJKlUEbHPsrlz57J48WKee+45vvKVr7Bjx45ut210vSIVGhQppSeBN7osngZsSim9nFL6DXAPcD7wW8CpwO8C8yKi29oi4sqIWB0Rq7ds2VJg9ZJ04M444wwefPBBfv3rX7N9+3ZWrFgBwPbt2znmmGPYuXMnd9111571hw8fzvbt2/e87mm9MjXj9NhRwCt1rzuBU1JKCwAiYi7wekrpne42TiktAZYAdHR0FHeTWFXO6IXf499e3rrn+f54L4vBqex/9ylTpnDxxRczadIkjj/+eGbMmAHA9ddfzymnnMLxxx/PxIkT94TDJZdcwrx587j55pu57777elyvTM0Iin33wWDPL/yU0rLySpGk4l1zzTVcc801+yyfP3/+PsumT5++1+mx8+fP73a9MjXjrKdO4Li61+3Aq71pICJmR8SSbdu29WthkqR9NSMoVgHjImJMRBwEXAIs700D3uFOkspT6KGniLgbOBMYGRGdwFdSSndExALgUWAosDSl1Kvzvbxndu/0dsoCj91roEkpdXvm0UCXUv8M4xYaFCmlOT0sfwR45ADaXQGs6OjomNfXNiQNDm1tbWzdupURI0YMqrBIKbF161ba2toOuC0nBZQ0oLW3t9PZ2clgPJ2+ra2N9vb2A26nkkHhoSdJjRo2bBhjxoxpdhmV1lJzPTXKwWxJKk8lg0KSVB6DQpKUq5JB4QV3klSeSgaFYxSSVJ5KBoUkqTwGhSQpVyWDwjEKSSpPJYPCMQpJKk8lr8yW+psTJ0o9q+QehSSpPJUMCscoJKk8lQwKxygkqTyVDApJUnkMCklSLoNCkpTLoJAk5TIoJEm5KhkUnh4rSeWp5JXZKaUVwIqOjo55za6lLF45LKlZKhkUUivpTYgb4KqiSh56kiSVx6CQJOUyKCRJuQwKSVIug0KSlKuSQeF1FJJUnkoGhdOMS1J5KhkUkqTyGBSSpFwGhSQpl0EhScplUEiSchkUkqRcBoUkKZfTjEtN4j1GVBXuUUiSchkUkqRcLRMUEXFiRNwaEfdFxPxm1yNJqik0KCJiaUS8FhHruyyfFREbI2JTRCwESCm9kFK6Cvgs0FFkXZKkxhW9R7EMmFW/ICKGArcA5wATgDkRMSF77zzgR8DfF1yXJKlBhQZFSulJ4I0ui6cBm1JKL6eUfgPcA5yfrb88pXQacGmRdUmSGteM02NHAa/Uve4ETomIM4ELgPcCj/S0cURcCVwJ8MEPfrC4KiVJQHOCIrpZllJKTwBP7G/jlNISYAlAR0dH6tfKJEn7aMZZT53AcXWv24FXe9OAd7iTpPI0IyhWAeMiYkxEHARcAizvTQPe4U6SylP06bF3Az8BxkdEZ0RckVLaBSwAHgVeAO5NKW3oZbvuUUhSSQodo0gpzelh+SPkDFg30O4KYEVHR8e8vrYhSWpMy1yZLUlqTZUMCg89SVJ5KhkUDmZLUnm8H0XJenMPAu8/oJ54LwuVqZJ7FJKk8lQyKByjkKTyVDIoHKOQpPJUMigkSeVpKCgi4n9GxGFRc0dEPB0RZxddnCSp+Rrdo7g8pfTvwNnAUcDvA18vrKr9cIxCksrTaFDsnhr8XODbKaVn6H668FI4RiFJ5Wk0KNZExGPUguLRiBgOvFNcWZKkVtHoBXdXAJOAl1NKb0XECGqHnyRJA1yjexSPp5SeTin9CiCltBX48+LKkiS1itw9iohoAw4GRkbEEbw7LnEYcGzBteXVNRuYPXbs2GaVIEmDxv72KP4bsAY4Ifu6+/EwcEuxpfXMwWxJKk/uHkVK6S+Av4iIL6aUFpVUkySphTQ0mJ1SWhQRpwGj67dJKf1VQXVJKoCzzqovGgqKiPhr4MPAOuDtbHECDApJGuAaPT22A5iQUkpFFiNJaj2Nnh67HvgvRRbSG07hIUnlaXSPYiTwfEQ8Bfzn7oUppfMKqWo/UkorgBUdHR3zmvH5HueVNJg0GhTXFVmEJKl1NXrW0z8VXYgkqTU1etbTdmpnOQEcBAwD3kwpHVZUYZKk1tDoHsXw+tcR8TvAtEIqkiS1lD7dCjWl9BBwVj/XIklqQY0eerqg7uUQatdVeE2FJA0CjZ71NLvu+S5gM3B+v1cjSWo5jY5RtNRNipxmXCqf1w8NXg2NUUREe0Q8GBGvRcQvIuL+iGgvurieOM24JJWn0cHsbwPLqd2saBSwIlsmSRrgGg2Ko1JK304p7coey4CjCqxLktQiGg2K1yPisogYmj0uA7YWWZgkqTU0GhSXA58F/g34OXAR0FID3JKkYjR6euz1wOdSSr8EiIgjgRupBYgkaQBrdI/io7tDAiCl9AYwuZiSJEmtpNGgGBIRR+x+ke1RNLo3IkmqsEZ/2f8ZsDIi7qM2dcdnga8WVpUkqWU0emX2X0XEamoTAQZwQUrp+UIrkyS1hIYPH2XBUGg4ZNOXfxo4GrglpfRYkZ8nSdq/Pk0z3hsRsTSb+mN9l+WzImJjRGyKiIVQm748pTQPmAtcXHRtkqT9KzwogGXArPoFETEUuAU4B5gAzImICXWr/HH2viSpyQoPipTSk8AbXRZPAzallF5OKf0GuAc4P2q+AXw/pfR00bVJkvavjD2K7owCXql73Zkt+yIwE7goIq7qbsOIuDIiVkfE6i1bthRfqSQNcs26FiK6WZZSSjcDN+dtmFJaAiwB6Ojo6PNd9pxbX5Ia06w9ik7guLrX7cCrjW4cEbMjYsm2bdv6vTBJ0t6aFRSrgHERMSYiDgIuoXa/i4Z44yJJKk8Zp8feDfwEGB8RnRFxRUppF7AAeBR4Abg3pbSh6FokSb1X+BhFSmlOD8sfAR7pS5veM1uqlt6MCToe2HqadejpgHjoSZLKU8mgcDBbkspTyaBwj0KSylPJoJAklcegkCTlqmRQOEYhSeWpZFA4RiFJ5alkUEiSymNQSJJyVTIoHKOQpPJUMigco5Ck8lQyKCRJ5TEoJEm5DApJUq5KBoWD2ZJUnkoGhYPZklSewm9cJEnN0JubJYE3TMpTyT0KSVJ5DApJUi6DQpKUq5JB4VlPklSeSgaFZz1JUnkqGRSSpPIYFJKkXAaFJCmXQSFJymVQSJJyGRSSpFzO9SRJXfRmnqjBMEdUJfcovOBOkspTyaDwgjtJKk8lg0KSVB6DQpKUy6CQJOUyKCRJuQwKSVIug0KSlMugkCTlMigkSbkMCklSrpYJioj4UETcERH3NbsWSdK7Cg2KiFgaEa9FxPouy2dFxMaI2BQRCwFSSi+nlK4osh5JUu8VvUexDJhVvyAihgK3AOcAE4A5ETGh4DokSX1UaFCklJ4E3uiyeBqwKduD+A1wD3B+kXVIkvquGWMUo4BX6l53AqMiYkRE3ApMjog/6mnjiLgyIlZHxOotW7YUXaskDXrNuHFRdLMspZS2Alftb+OU0hJgCUBHR0fq59okSV00Y4+iEziu7nU78GpvGvDGRZJUnmYExSpgXESMiYiDgEuA5b1pwBsXSVJ5ij499m7gJ8D4iOiMiCtSSruABcCjwAvAvSmlDUXWIUnqu0LHKFJKc3pY/gjwSF/bjYjZwOyxY8f2tQlJUoNa5srs3vDQkySVp5JB4WC2JJWnkkHhHoUklaeSQSFJKo9BIUnK1Ywrsw+YZz1JakWjF36vV+tv/vqnC6qkf1Vyj8IxCkkqTyWDQpJUHoNCkpSrkkHhdRSSVJ5KBoVjFJJUnkoGhSSpPAaFJCmXQSFJyuUFd5LUAlr5Yr1K7lE4mC1J5alkUEiSymNQSJJyGRSSpFwGhSQpl2c9SWpZrXwm0GBSyT0Kz3qSpPJUMigkSeUxKCRJuQwKSVIug0KSlMugkCTlMigkSbkMCklSrkoGhffMlqTyREqp2TX0WURsAf61xI8cCbxe4uc102Dp62DpJ9jXgaiv/Tw+pXRUoytXOijKFhGrU0odza6jDIOlr4Oln2BfB6Ky+lnJQ0+SpPIYFJKkXAZF7yxpdgElGix9HSz9BPs6EJXST8coJEm53KOQJOUyKCRJuQyKBkTE0oh4LSLWN7uWInTXv4g4MiIej4iXsq9HNLPGvupN36Lm5ojYFBHPRsSU5lXemP7qX0R8Llv/pYj4XDP60lXRfYuIj0fEc9k2N0dEDPS+9fnnOqXkYz8P4AxgCrC+2bWU1T/gT4GF2fOFwDeaXWfRfQPOBb4PBHAq8M/Nrr+M/gFHAi9nX4/Inh8x0PsGPAV8Itvm+8A5A71vff25bvp/9Ko8gNEDNSi66x+wETgme34MsLHZNRbdN+A2YE5367Xy40D7B8wBbqtbvtd6A7Fv2Xv/r275XusN1L719efaQ0/qyQdSSj8HyL4e3eR6+lNPfRsFvFK3Xme2rGp6278q9bu/+jYqe951eTOV0bc+/VwbFNK7ujtGPZDOH++pfwOh373tW5X63PS+GRTqyS8i4hiA7OtrTa6nP/XUt07guLr12oFXS66tP/S2f1Xqd3/1rTN73nV5M5XRtz79XBsU6slyYPdZFJ8DHm5iLf2tp74tB34vO8vkVGDb7t30iult/x4Fzo6II7KzYM7OlrWifulb9t72iDg1OyPo92j+//Ey+ta3n+tmD1ZV4QHcDfwc2Ektra9odk1F9w8YAfw98FL29chm11l036jtst8C/BR4Duhodv1l9Q+4HNiUPX6/2f0qo29AB7A+22Yx2UwVA7lvff25dgoPSVIuDz1JknIZFJKkXAaFJCmXQSFJymVQSJJyGRSqrIh4IiIKv7F8RPyPiHghIu7q53a/FRETct6fGxHHHkD7kyLi3L5u36WtZRFxUfa8lO+7WodBoUEpIt7Ti9W/AJybUrq0n9oDIKX0+ZTS8zmrzAX6HBTAJGozj0oHxKBQoSJidPbX+O0RsSEiHouI92Xv7fnLNCJGRsTm7PnciHgoIlZExL9ExIKI+MOIWBsR/zcijqz7iMsiYmVErI+Iadn2h2Tz/a/Ktjm/rt3vRsQK4LFuav3DrJ31EXF1tuxW4EPA8oj4gy7r79NeRHw5+9xnI+L/1NXzvYh4Jmv74vr+R8TQ7C/29dk9BP4g++u9A7grItZFxPsi4n9nba+PiCXZVbe72/lGRDwVES9GxIyIOAj4E+DibPuLu9Q+NCJuzD7v2Yj4Yrb84xHxTxGxJiIejWy6hx7+bfepe///I1RJzb760sfAflCbSnkXMCl7fS9wWfb8CbKrTIGRwObs+VxqV5oOB44CtgFXZe/9OXB13fa3Z8/PIJuyGfha3We8H3gROCRrt5NurkYFPk7tqtdDgEOBDcDk7L3NwMhuttmrPWpTKCyhdiXtEODvsrou3F1ntt7h9f3PPvvxuvff3/X7k70+su75XwOz69b7s+z5ucAP6upb3MO/y3zgfuA9u9sGhgErgaOyZRcDS7Pny4CLGqnbx8B7uEehMvxLSmld9nwNtfDYn39MKW1PKW2hFhQrsuXPddn+boCU0pPAYRHxfmq/sBdGxDpqv9TagA9m6z+eUnqjm887HXgwpfRmSuk/gAeAGQ3UWd/e2dljLfA0cAIwLqt5ZvZX/4yU0rYubbwMfCgiFkXELODfe/is34qIf46I54CzgJPq3nsg+9ro93cmcGtKaRdA1ofxwMnA49n37o/Ze3K5rhqtWxXX6+OqUh/8Z93zt4H3Zc938e7hz7acbd6pe/0Oe/+/7ToHze5pli9MKW2sfyMiTgHe7KHGvt4Gs769AG5IKd22T+MRH6f21/4NEfFYSulP9hSc0i8j4mPAp4D/DnyW2hw+9du3AX9JbQ/jlYi4jr2/Z7u/P2/T2M91sO/3LoANKaVPNLB9Q3VrYHCPQs20mdrhC4CL+tjG7uP9p1ObVXMbtVk1v1h3DH9yA+08CfxORBwcEYcAnwF+2MtaHgUuj4hDs88dFRFHZ2cuvZVS+hvgRmq3wNwjIkYCQ1JK9wPX1r2/ndrhN3g3FF7P2m/k+1W/fVePAVftHoTPxn02AkdFxCeyZcMi4qQets+rWwOMexRqphuBeyPivwL/0Mc2fhkRK4HDePev2euBm4Bns7DYDPx2XiMppacjYhm1ew0DfCultLY3haSUHouIE4GfZBn1H8BlwFjgmxHxDrXZQud32XQU8O2I2P2H2x9lX5cBt0bEr6nd//h2aoexNgOrGijpH3n3ENwNKaW/rXvvW8BHqH2PdlIbQ1mcDaLfHBGHU/v9cBO18Zru9FS3Bhhnj5Uk5fLQkyQpl0EhScplUEiSchkUkqRcBoUkKZdBIUnKZVBIknL9f+sROE8u6PAzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "results = np.array(results)\n",
    "results_nz = results[np.nonzero(results)[0]]    # remove zeros from results so that we can take the log\n",
    "num_zero = repeats - len(results_nz)    # number of zeros in the original results\n",
    "bins = np.arange(-0.125, 5, 0.25)    # specify bin edges for histogram\n",
    "hist = np.histogram(np.log10(results_nz), bins)[0]    # calculate histogram in log scale\n",
    "hist = np.hstack(([num_zero], hist))    # add back number of zeros into the histogram\n",
    "\n",
    "xticks = [1, 10, 100, 1000, 10000, 100000]\n",
    "plt.figure()\n",
    "plt.bar(bins, hist, width=0.2, label='data')\n",
    "plt.axvline(np.log10(mean), color='k', label='mean')\n",
    "plt.xticks(np.log10(xticks), xticks)\n",
    "plt.xlabel('number of resistant cells')\n",
    "plt.yscale('log')\n",
    "plt.ylabel('counts')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "The long tail (again, in log scale!) of the distribution shows that, in some experimental trials, we get a huge number of resistant cells compared to the mean, which contributes to the large variance that we saw above. The exact shape of the distribution is not easy to calculate theoretically. However, it suffices to use the very large variance-to-mean ratio to distinguish this scenario from the first one. Indeed, it was by observing such a large variance-to-mean ratio in their real experiments that Luria & Delbruck (1943) were able to reject the first scenario (mutation after exposure) and establish the second scenario (mutation before exposure). Remember that this was before people knew the molecular structrue of DNA and the cause of mutations!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}