lecture-04.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lecture 04: Combining natural selection with other forces of evolution\n",
    "\n",
    "### Natural selection and genetic drift\n",
    "\n",
    "We saw that natural selection and genetic drift can interact in complicated ways. When an allele is very rare (for example, a new mutation), its initial fate is determined almost entirely by drift. Once the allele grows larger than an establishment frequency ($\\sim 1/Ns$ following [Desai and Fisher](https://academic.oup.com/genetics/article/176/3/1759/6062236), but keep in mind that this is a fuzzy boundary rather than a sharp one) its future change is determined more strongly by selection. In particular, alleles for which $2Ns\\gg 1$ will grow rapidly.\n",
    "\n",
    "Here, we'll explore the interaction between drift and selection in simulations. To do this, we'll use the Wright-Fisher model that we considered previously, but adapted to take natural selection into account. To make our analysis simpler, we'll also assume that we're working with a haploid population rather than diploids, allowing us to avoid complications due to dominance.\n",
    "\n",
    "As before, let's call two types of alleles A and B, with frequencies $p$ and $1-p$. Under the Wright-Fisher model **without selection**, the probability of getting $n$ alleles of type A in the next generation is binomial,\n",
    "\n",
    "$$ P(n) = \\binom{N}{n} p^n \\left(1-p\\right)^{N-n}\\,. $$\n",
    "\n",
    "Since we're considering haploids, the total number of alleles is equal to the population size $N$ instead of $2N$. \n",
    "\n",
    "To account for selection, we simply adjust the probability $p$ that appears in the binomial distribution above. Instead of depending only on the allele frequency, we now weight the probability by fitness. If the fitness of the A allele is $1+s$ and the fitness of the B allele is $1$, then the probability to select an A allele $p_A$ is \n",
    "\n",
    "$$ p_A = \\frac{p (1 + s)}{\\bar{w}} = \\frac{p (1 + s)}{p (1 + s) + (1 - p)} \\,. $$\n",
    "\n",
    "With this adjustment, the probability of obtaining $n$ A alleles in the next generation is \n",
    "\n",
    "$$ P(n) = \\binom{N}{n} p_A^n \\left(1-p_A\\right)^{N-n}\\,. $$\n",
    "\n",
    "Let's use this expression to simulate the number of A alleles over time in a population. We'll run the simulation multiple times and plot the results to get a sense of the randomness involved in this process.\n",
    "\n",
    "(**Note:** We will not consider this right now, but mutation could also be incorporated into the Wright-Fisher model in a similar way as we have done above for selection by modifying $p_A$.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np          # here we import numpy\n",
    "import numpy.random as rng  # and here we import the random number generation (sub-)library\n",
    "\n",
    "\n",
    "# Set the starting parameters for the simulations\n",
    "\n",
    "N   = 100   # population size\n",
    "s   = 0.01  # selection coefficient\n",
    "n_A = 1     # starting number of A alleles\n",
    "n_g = 100   # number of generations of evolution\n",
    "n_r = 50    # number of replicate simulations\n",
    "n_lost = 0  # number of simulations in which the A allele is lost\n",
    "n_t = np.zeros((n_r, n_g))  # a blank matrix of replicates x allele numbers over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Next, let's fill in the vector of allele frequencies over time by simulating the WF model\n",
    "\n",
    "n_lost = 0\n",
    "\n",
    "for i in range(n_r):\n",
    "\n",
    "    n_t[i][0] = n_A  # set the starting frequency to p\n",
    "\n",
    "    for j in range(n_g-1):\n",
    "        p_A = n_t[i][j] * (1 + s) / (n_t[i][j] * (1 + s) + N - n_t[i][j])  # probability to select an A allele\n",
    "        n_t[i][j+1] = rng.binomial(N, p_A)  # get As in next generation\n",
    "        \n",
    "    if n_t[i][-1]==0:\n",
    "        n_lost += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finally, let's make a plot\n",
    "\n",
    "import seaborn as sns            # import seaborn\n",
    "import matplotlib.pyplot as plt  # and matplotlib\n",
    "\n",
    "\n",
    "# If the establishment number is nonnegative and < N, plot it here\n",
    "\n",
    "n_estab = 1 / s\n",
    "if n_estab > 0 and n_estab < N:\n",
    "    plt.hlines(y=n_estab, xmin=0, xmax=n_g, label='establishment', color='black', ls='--')\n",
    "\n",
    "\n",
    "# Plot individual simulations\n",
    "\n",
    "for i in range(n_r):\n",
    "    sns.lineplot(x=np.arange(n_g), y=n_t[i], alpha=0.5)\n",
    "    \n",
    "\n",
    "plt.xlabel('Time (generations)')\n",
    "plt.ylabel('Number of A alleles')\n",
    "plt.yscale('log')\n",
    "plt.ylim(1, 1.1*N);\n",
    "\n",
    "\n",
    "# Let's also write out the fraction of simulations in which A alleles were lost\n",
    "\n",
    "print('At the final time point, there are no A alleles in %d out of %d simulations.' % (n_lost, n_r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Regimes of selection.** Recall that $Ns$ (or $2Ns$ for diploids) is a critical parameter for undertsanding the interaction between selection and drift. \n",
    "\n",
    "- If $Ns\\ll 1$, then **drift is dominant**.\n",
    "\n",
    "- If $Ns \\sim 1$, then **drift and selection are comparable**, even at high allele frequencies.\n",
    "\n",
    "- If $Ns \\gg 1$, then **selection is dominant** at high allele frequencies.\n",
    "\n",
    "To reinforce our understanding, let's return to the simulation above. Try different values for $N$ and $s$ that fit within these regimes. Do they behave in the way that you expect? \n",
    "\n",
    "In addition, explore simulations for deleterious alleles with $s<0$. Do they appear to follow rules that are similar to those for beneficial alleles?\n",
    "\n",
    "You may also want to consider how the evolution changes when you change the starting number of alleles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.9"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}