{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.set_printoptions(suppress=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What Bruce Levin says is [@Turrientes2013] S1 File:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p0 = 1e-5\n",
      "s = 0.001\n",
      "pe = 0.1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "S = array([[1,0],[0,1-s]])\n",
      "p = array([p0, 1-p0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for t in range(100000):\n",
      "    if p[0] > pe:\n",
      "        break\n",
      "    p = S.dot(p)\n",
      "    p /= p.sum()\n",
      "print t "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "9312\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "What I say:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 500 # max number of mutations\n",
      "mu = 0.0004 # Wielgoss2011\n",
      "tau = 35 # mutator rate increase\n",
      "s = 0.001  # selection per site\n",
      "print \"mutation rate:\",mu\n",
      "print \"selection coef:\",s\n",
      "print \"lowest fitness:\",(1-s)**n\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "mutation rate: 0.0004\n",
        "selection coef: 0.001\n",
        "lowest fitness: 0.606378944861\n"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "considering only single deleterious mutations the mutation matrix is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "zeros_matrix = lambda n: matrix(zeros((n,n)))\n",
      "mutation_matrix = lambda n,mu: diag([1-mu] * (n+1)) + diag([mu] * n,-1)\n",
      "A = mutation_matrix(n,mu)\n",
      "B = zeros_matrix(n+1)\n",
      "C = zeros_matrix(n+1)\n",
      "D = mutation_matrix(n,tau * mu)\n",
      "M = np.bmat([[A, B], [C, D]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "the selection matrix is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "selection_matrix =  lambda n,s: diag([(1-s)**i for i in range(n+1)])\n",
      "A = selection_matrix(n,s)\n",
      "B = zeros_matrix(n+1)\n",
      "C = zeros_matrix(n+1)\n",
      "D = selection_matrix(n,s)\n",
      "S = np.bmat([[A, B], [C, D]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "the evolution matrix is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "E = array(M*S)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "now for the evolution:"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "starting from mutation-free:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p = zeros((n+1)*2)\n",
      "p[0] = p0\n",
      "p[n+1] = 1-p0\n",
      "print p"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.00001  0.       0.      ...,  0.       0.       0.     ]\n"
       ]
      }
     ],
     "prompt_number": 37
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "starting from mutation-selection balance:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import poisson"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p = array([p0*poisson(mu/s).pmf(i) for i in range(n+1)] + [(1-p0)*poisson(tau*mu/s).pmf(i) for i in range(n+1)])\n",
      "print p"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.0000067   0.00000268  0.00000054 ...,  0.          0.          0.        ]\n"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "mean fitness at the mutation-selection balance:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sv = array([(1-s)**i for i in range(n+1)])\n",
      "w_nm = sv.dot(p[:n+1])/p0\n",
      "w_cm = sv.dot(p[n+1:])/(1-p0)\n",
      "print w_nm, w_cm"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.999600079989 0.986097544263\n"
       ]
      }
     ],
     "prompt_number": 40
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "evolve until NM reach `pe`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for t in range(1000000):\n",
      "    if p[0] > pe:\n",
      "        break\n",
      "    p = E.dot(p)\n",
      "    p /= p.sum()\n",
      "print t "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "718\n"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "this is an order of magnitude lower than what Levin the result of Levin (~9000)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Disclaimer\n",
      "\n",
      "This notebook is intended for internal use. The author doesn't stand behind anything written here and no disrespect is intended."
     ]
    }
   ],
   "metadata": {}
  }
 ]
}