{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Using Symbolic Math to Understand Dynamics of Natural Selection Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Table of contents**\n", " * [Example code for codominant selection](#section1)\n", " * [Exercises](#section2)\n", " * [$A_{2}$ dominant](#type1)\n", " * [$A_{2}$ advantageous](#type1a)\n", " * [$A_{2}$ deleterious](#type1b)\n", " * [$A_{2}$ recessive](#type2)\n", " * [$A_{2}$ advantageous](#type2a)\n", " * [$A_{2}$ deleterious](#type2b)\n", " * [Overdominance](#type3)\n", " * [Underdominance](#type4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example code for codominant selection\n", "\n", "For this lab, we will use symbolic math with the Python package [SymPy](http://docs.sympy.org/dev/index.html) to find the fixed points and their stability for the model\n", "\n", "$$q_{t+1} = F(q_{t}) = \\frac{w_{12}q_{t}(1-q_{t}) + w_{22}q_{t}^2}{w_{11}(1-q_{t})^2 + 2w_{12}q_{t}(1-q_{t}) + w_{22}q_{t}^2}$$\n", "\n", "To begin, we need to first load the SymPy package into the workspace." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sympy import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next step is to use the `var` command to create symbolic variables. We will need $q, w_{11}, w_{12}, w_{22}, s$ and $t$. After this command, we will specify the function $F$." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(q**2*w22 + q*w12*(-q + 1))/(q**2*w22 + 2*q*w12*(-q + 1) + w11*(-q + 1)**2)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "q,w12,w22,w11,s,t = var('q,w12,w22,w11,s,t')\n", "F = ((w12*(1-q)*q + w22*q*q)/((1-q)*(1-q)*w11 + 2*w12*(1-q)*q + w22*q*q))\n", "F" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we will create a dictionary relating the weights as functions of $s$ and possibly $t$. We'll start with **codominance** for this walkthrough." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [], "source": [ "selection_type = {w11: 1, w12: 1+s, w22: 1+2*s}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To make things easier, let's go ahead and simplify $F$ by substituting the correct values for the weights (given our choice of selection). We do that using the `subs` command as follows:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(q**2*(2*s + 1) + q*(-q + 1)*(s + 1))/(q**2*(2*s + 1) + 2*q*(-q + 1)*(s + 1) + (-q + 1)**2)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F = F.subs(selection_type)\n", "F" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that `**` means _exponent_, i.e. `q**2` is the same as $q^2$. Let's also try and simplify this expression, which we can do using the `simplify` command." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "q*(q*s + s + 1)/(2*q*s + 1)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "F = F.simplify()\n", "F" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's better. To compute fixed points for this system, we need to `solve` for $q$ in the equation\n", "\n", "$$F(q) = q$$\n", "\n", "which is equivalent to \n", "\n", "$$F(q)-q = 0.$$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "These are the fixed points:\n", "[0, 1]\n" ] } ], "source": [ "fps = solve(F-q,q)\n", "print \"These are the fixed points:\"\n", "print fps" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, so we've got two fixed points: $q_{\\infty} = 0$ and $q_{\\infty} = 1$. Our job now is to compute the stability of these fixed points. First thing we need to do is compute the (symbolic) derivative $F^{\\prime}(q)$. We do that using the `diff` function. The first argument to the function is $F$, and the second argument specifies wich variable to take a derivative with respect to. We only have one variable $q$, so we'll specify that one." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "q*s/(2*q*s + 1) - 2*q*s*(q*s + s + 1)/(2*q*s + 1)**2 + (q*s + s + 1)/(2*q*s + 1)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dF = diff(F,q)\n", "dF" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Holy smokes! That is complicated. Let's use the simplify command again." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(2*q**2*s**2 + 2*q*s + s + 1)/(4*q**2*s**2 + 4*q*s + 1)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dF = dF.simplify()\n", "dF" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Uh, I guess that's better? It's time now to substitute specific values for $q$ for the two fixed points." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "dF/dq at q=0: s + 1\n", "dF/dq at q=1: (s + 1)/(2*s + 1)\n" ] } ], "source": [ "dF_q0 = dF.subs(q,0).simplify()\n", "dF_q1 = dF.subs(q,1).simplify()\n", "print \"dF/dq at q=0: \", dF_q0\n", "print \"dF/dq at q=1: \", dF_q1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So our question now is: If allele $A_{2}$ has a selective advantage (i.e. $s > 0$), what is the stability of the fixed points? Remember we have the following stability criteria:\n", "\n", "$$\n", "\\begin{align}\n", "\\textrm{stable}: & \\ |F^{\\prime}(q_{\\infty})| < 1 \\\\\n", "\\textrm{unstable}: & \\ |F^{\\prime}(q_{\\infty})| > 1 \\\\\n", "\\textrm{neutral}: & \\ |F^{\\prime}(q_{\\infty})| = 1\n", "\\end{align}\n", "$$\n", "\n", "From our results above, we see that\n", "\n", "$$\n", "\\begin{align}\n", "|F^{\\prime}(0)| = |s+1| > 1 & \\Rightarrow \\textrm{$q_{\\infty} = 0$ is unstable} \\\\\n", "|F^{\\prime}(1)| = \\left|\\frac{s+1}{2s+1}\\right| < 1 & \\Rightarrow \\textrm{$q_{\\infty} = 1$ is stable}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In case the first derivative is 1, in which case we are uncertain of stability, we can compute the second derivative and substitute the appropriate $q$ value." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Here is the second derivative, in case the first derivative equals 1: -2*s/(4*s**2 + 4*s + 1)\n" ] } ], "source": [ "ddF = diff(dF,q)\n", "ddF_q0 = ddF.subs(q,1).simplify()\n", "print \"Here is the second derivative, in case the first derivative equals 1:\", ddF_q0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises\n", "\n", "Repeat for all selection types under each condition, listed again here for convenience:\n", "\n", "Type of selection | $w_{11}$ | $w_{12}$ | $w_{22}$ | Conditions\n", "----------------- | -------- | -------- | -------- | ----------\n", "$A_{2}$ dominant | 1 | 1+s | 1+s | $s > 0$ (advantageous), $s < 0$ (deleterious)\n", "$A_{2}$ recessive | 1 | 1 | 1+s | $s > 0$ (advantageous), $s < 0$ (deleterious)\n", "Overdominance | 1 | 1+s | 1+t | $s > 0$ and $s > t$\n", "Underdominance | 1 | 1+s | 1+t | $s < 0$ and $s < t$\n", "\n", "For each of the types listed above, it may be helpful to use the visualization at http://mathinsight.org/cobwebbing_graphical_solution to gain insight into the system and verify your results." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $A_{2}$ dominant" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $A_{2}$ advantageous" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $A_{2}$ deleterious" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### $A_{2}$ recessive" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $A_{2}$ advantageous" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### $A_{2}$ deleterious" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Overdominance" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Underdominance" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2 (SageMath)", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 0 }