{
"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
}