male-genotype-fitness.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import ipywidgets as widgets\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import integrate\n",
    "import sympy as sym\n",
    "\n",
    "import models\n",
    "import selection_functions\n",
    "import symbolics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.style.use(\"ggplot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Monomorphic $\\gamma$ model\n",
    "\n",
    "## Only G allele present \n",
    "\n",
    "Set $x_2=1-x_1$, $x_3=x4=0$ so that only the G allele of the $\\gamma$ gene is present"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x1, x2, x3 = sym.symbols(\"x1, x2, x3\", real=True, nonnegative=True)\n",
    "x4 = 1 - x1 - x2 - x3\n",
    "T, R, P, S = sym.symbols('T, R, P, S', real=True, positive=True)\n",
    "M, m = sym.symbols(\"M, m\", real=True, nonnegative=True)\n",
    "epsilon = sym.symbols(\"epsilon\", real=True, positive=True)\n",
    "\n",
    "UGA = symbolics.UGA\n",
    "UgA = symbolics.UgA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x = np.array([[x1], [x2], [x3], [x4]])\n",
    "payoff_kernel = np.array([[R, S], [T, P]])\n",
    "W = models.generalized_sexual_selection(x, UGA, UgA, payoff_kernel, M, m, epsilon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "male_genotype_fitness = W.sum(axis=0).sum(axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "male_GA_fitness = sym.together(sym.expand(male_genotype_fitness[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.0*x1*(-M*UGA(x1 + x3) + R*UGA(x1 + x3)**2 + R*UgA(x1 + x3)**2 - S*UGA(x1 + x3)**2 + S*UGA(x1 + x3) - S*UgA(x1 + x3)**2 + S*UgA(x1 + x3) - m*UgA(x1 + x3))/(x1 + x3)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "male_GA_fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitness of males carrying the $GA$ genotype.\n",
    "\n",
    "$$ 4\\big((R - S)U_{GA}(x_1)^2 + (S - M)U_{GA}(x_1)\\big) + 4\\big((R - S)U_{gA}(x_1)^2 + (S - m)U_{gA}(x_1)\\big) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "male_Ga_fitness = sym.together(sym.expand(male_genotype_fitness[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x2*(4.0*M*UGA(x1 + x3) - 4.0*M + 4.0*P*UGA(x1 + x3)**2 - 8.0*P*UGA(x1 + x3) + 4.0*P*UgA(x1 + x3)**2 - 8.0*P*UgA(x1 + x3) + 8.0*P - 4.0*T*UGA(x1 + x3)**2 + 4.0*T*UGA(x1 + x3) - 4.0*T*UgA(x1 + x3)**2 + 4.0*T*UgA(x1 + x3) + 4.0*m*UgA(x1 + x3) - 4.0*m)/(-x1 - x3 + 1)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "male_Ga_fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitness of males carrying the $Ga$ genotype.\n",
    "\n",
    "$$ 4\\big(2P - (M + m)\\big) + 4\\big((T + M - 2P)U_{GA}(x_1) - (T - P)U_{GA}(x_1)^2\\big) + 4\\big((T + m - 2P)U_{gA}(x_1) - (T - P)U_{gA}(x_1)^2\\big) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.0*x1*(-M*UGA(x1 + x3) + R*UGA(x1 + x3)**2 + R*UgA(x1 + x3)**2 - S*UGA(x1 + x3)**2 + S*UGA(x1 + x3) - S*UgA(x1 + x3)**2 + S*UgA(x1 + x3) - m*UgA(x1 + x3))/(x1 + x3) - x2*(4.0*M*UGA(x1 + x3) - 4.0*M + 4.0*P*UGA(x1 + x3)**2 - 8.0*P*UGA(x1 + x3) + 4.0*P*UgA(x1 + x3)**2 - 8.0*P*UgA(x1 + x3) + 8.0*P - 4.0*T*UGA(x1 + x3)**2 + 4.0*T*UGA(x1 + x3) - 4.0*T*UgA(x1 + x3)**2 + 4.0*T*UgA(x1 + x3) + 4.0*m*UgA(x1 + x3) - 4.0*m)/(-x1 - x3 + 1)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.collect(male_GA_fitness - male_Ga_fitness, (UGA(x1), UgA(x1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fitness differential is\n",
    "\n",
    "$$ 4\\big((T + R - P - S)U_{GA}(x_1)^2 - (T + 2M - 2P - S)U_{GA}(x_1)\\big) + 4\\big((T + R - P - S)U_{gA}(x_1)^2 - (T + 2m - 2P - S)U_{gA}(x_1)\\big) - 4\\big(2P - (M + m)\\big) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_male_genotype_fitness(x1, x2, x3, selection_function, d1, d3,\n",
    "                               T, R, P, S, M, m, epsilon, max_time):\n",
    "\n",
    "    fig, axes = plt.subplots(1,2, figsize=(15,8), sharex=True)\n",
    "\n",
    "    # prepare the axes\n",
    "    axes[0].set_ylim((0, 1.05))\n",
    "    axes[0].set_xlabel(r\"Time, $t$\", fontsize=15)\n",
    "    axes[0].set_ylabel(r\"Offspring genotype shares, $x_i$\", fontsize=15)\n",
    "\n",
    "    axes[1].set_xlabel(r\"Time, $t$\", fontsize=15)\n",
    "    axes[1].set_ylabel(r\"Male genotype fitness\", fontsize=15)\n",
    "\n",
    "    # create the initial condition\n",
    "    x4 = 1 - x1 - x2 - x3\n",
    "    y0 = np.array([x1, x2, x3, x4])\n",
    "    assert y0.sum() <= 1\n",
    "\n",
    "    # create the payoff kernel\n",
    "    assert (T > R) and (R > P) and (R > S), \"Payoffs must satisfy either Prisoner's Dilemma or Stag Hunt constraints! T={}, R={}, P={}, S={}\".format(T, R, P, S)\n",
    "    payoff_kernel = np.array([[R, S], [T, P]])\n",
    "\n",
    "    # create the selection functions\n",
    "    if selection_function == \"kirkpatrick\":\n",
    "        UGA = lambda x_A: selection_functions.kirkpatrick_selection(x_A, d1)\n",
    "        UgA = lambda x_A: selection_functions.kirkpatrick_selection(x_A, d3)\n",
    "    elif selection_function == \"seger\":\n",
    "        UGA = lambda x_A: selection_functions.seger_selection(x_A, d1)\n",
    "        UgA = lambda x_A: selection_functions.seger_selection(x_A, d3)\n",
    "    elif selection_function == \"wright\":\n",
    "        UGA = lambda x_A: selection_functions.wright_selection(x_A, d1)\n",
    "        UgA = lambda x_A: selection_functions.wright_selection(x_A, d3)\n",
    "    else:\n",
    "        valid_funcs = (\"kirkpatrick\", \"seger\", \"wright\")\n",
    "        msg = \"Selection_function must be one of {}, {}, or {}.\".format(*valid_funcs)\n",
    "        raise ValueError(msg)\n",
    "\n",
    "    # simulate the model starting from a random initial condition\n",
    "    def f(t, y):\n",
    "        W = models.generalized_sexual_selection(y, UGA, UgA, payoff_kernel, M, m, epsilon)\n",
    "        y_dot = models.offspring_genotypes_evolution(W, y)\n",
    "        return y_dot\n",
    "\n",
    "    solution = integrate.solve_ivp(f, t_span=(0, max_time), y0=y0, method=\"RK45\",\n",
    "                                   rtol=1e-9, atol=1e-12, dense_output=True, vectorized=True)\n",
    "\n",
    "    axes[0].plot(solution.t, solution.y[0], label=\"GA\")\n",
    "    axes[0].plot(solution.t, solution.y[1], label=\"Ga\")\n",
    "    axes[0].plot(solution.t, solution.y[2], label=\"gA\")\n",
    "    axes[0].plot(solution.t, solution.y[3], label=\"ga\")\n",
    "    axes[0].legend()\n",
    "\n",
    "    def _male_genotype_fitness(yt):\n",
    "        W = models.generalized_sexual_selection(yt, UGA, UgA, payoff_kernel, M, m, epsilon)\n",
    "        male_genotype_fitness = W.sum(axis=0).sum(axis=0).reshape(-1, 1)\n",
    "        return male_genotype_fitness\n",
    "\n",
    "    def male_genotype_fitness(y):\n",
    "        _, T = y.shape\n",
    "        ws = []\n",
    "        for t in range(T):\n",
    "            yt = y[:,[t]]\n",
    "            w = _male_genotype_fitness(yt)\n",
    "            ws.append(w)\n",
    "        return np.hstack(ws)\n",
    "\n",
    "    fitness = male_genotype_fitness(solution.y)\n",
    "    axes[1].plot(solution.t, fitness[0], label=\"GA\")\n",
    "    axes[1].plot(solution.t, fitness[1], label=\"Ga\")\n",
    "    axes[1].plot(solution.t, fitness[2], label=\"gA\")\n",
    "    axes[1].plot(solution.t, fitness[3], label=\"ga\")\n",
    "    axes[1].legend()\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "    return (solution, fitness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5ea781fc92ae4b9d9e9874dc83df147f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sliders used to control the initial condition\n",
    "x1_slider = widgets.FloatSlider(value=0.25, min=0.0, max=1.0, step=1e-3, description=r\"$x_1$\", readout_format=\".3f\")\n",
    "x2_slider = widgets.FloatSlider(value=0.25, min=0.0, max=1.0, step=1e-3, description=r\"$x_2$\", readout_format=\".3f\")\n",
    "x3_slider = widgets.FloatSlider(value=0.25, min=0.0, max=1.0, step=1e-3, description=r\"$x_3$\", readout_format=\".3f\")\n",
    "\n",
    "# sliders used to control the Prisoner's Dilemma Payoffs\n",
    "T_slider = widgets.FloatSlider(value=10, min=0, max=100, step=0.1, description=r\"$T$\")\n",
    "R_slider = widgets.FloatSlider(value=8, min=0, max=100, step=0.1, description=r\"$R$\")\n",
    "P_slider = widgets.FloatSlider(value=6, min=0, max=100, step=0.1, description=r\"$P$\")\n",
    "S_slider = widgets.FloatSlider(value=4, min=0, max=100, step=0.1, description=r\"$S$\")\n",
    "\n",
    "# sliders used to control the metabolic costs\n",
    "M_slider = widgets.FloatSlider(value=1, min=0, max=100, step=0.1, description=r\"$M_G$\")\n",
    "m_slider = widgets.FloatSlider(value=0, min=0, max=100, step=0.1, description=r\"$m_g$\")\n",
    "\n",
    "# slider used to control which selection function is being used\n",
    "U_slider = widgets.Dropdown(options=[\"kirkpatrick\", \"seger\", \"wright\"], index=0, description=r\"$U_{\\gamma(j)A}$\")\n",
    "\n",
    "# slider that controls the parameters of the selection function\n",
    "d1_slider = widgets.FloatSlider(value=1, min=0.0, max=10, step=0.05, description=r\"$d_1$\")\n",
    "d3_slider = widgets.FloatSlider(value=1, min=0.0, max=10, step=0.05, description=r\"$d_3$\")\n",
    "\n",
    "# slider used to control the mutation rate\n",
    "e_slider = widgets.FloatSlider(value=0.0, min=0.0, max=1.0, step=1e-3, description=r\"$\\epsilon$\", readout_format=\".3f\")\n",
    "\n",
    "# slider that controls max simulation time\n",
    "max_time_slider = widgets.IntSlider(value=25, min=1, max=100000, description=r\"$\\max t$\")\n",
    "\n",
    "w = widgets.interactive(plot_male_genotype_fitness, x1=x1_slider, x2=x2_slider, x3=x3_slider,\n",
    "                        selection_function=U_slider, d1=d1_slider, d3=d3_slider, \n",
    "                        T=T_slider, R=R_slider, P=P_slider, S=S_slider,\n",
    "                        M=M_slider, m=m_slider, epsilon=e_slider,\n",
    "                        max_time=max_time_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}