mt_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling microtubule interactions\n",
    "\n",
    "\\- Tamas Nagy\n",
    "\n",
    "Building on [Sumino et al](https://www.nature.com/nature/journal/v483/n7390/full/nature10874.html#supplementary-information)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib osx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "import scipy.spatial.distance as scidist\n",
    "from IPython.display import HTML\n",
    "from math import sqrt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameter setup\n",
    "\n",
    "These are the parameters that I extracted from the paper, the main ones to changes are $A$, the grid size, and $\\lambda$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Particle number\n",
    "N = 1500\n",
    "\n",
    "# particle speed (μm/s)\n",
    "v0 = 10 #8.75\n",
    "\n",
    "# trajectory curvature (1/μm)\n",
    "k0 = -7.1e-04 # mean\n",
    "k_sd = 1.8e-3 # std dev\n",
    "\n",
    "# relaxation time (s)\n",
    "τ = 62.0\n",
    "\n",
    "# microtubule length (μm)\n",
    "l = 40 #20\n",
    "\n",
    "# non-dimensionalized relaxation time \n",
    "λ = 100 #v0*τ/l\n",
    "\n",
    "# preferred angular velocity\n",
    "ω0 = v0 * k0\n",
    "\n",
    "# initial angular velocities\n",
    "ω = v0 * np.random.normal(k0, k_sd, N)\n",
    "\n",
    "# preferred non-dimensionalized angular velocity\n",
    "Ω0 = ω0 * l / v0\n",
    "\n",
    "# initial non-dimensionalized angular velocities\n",
    "Ω = ω * l / v0\n",
    "\n",
    "# white Gaussian noise parameters\n",
    "ξ_μ = 0\n",
    "ξ_σ = sqrt(2 * v0**2 * k_sd**2 / τ)*8\n",
    "\n",
    "# size of microtubule grid in microns\n",
    "grid_size = 100\n",
    "\n",
    "# weight of neighbors' angular velocity versus self\n",
    "A = 0.8 #0.1\n",
    "\n",
    "# random initial orientation\n",
    "θ = 2 * np.pi * np.random.rand(N)\n",
    "\n",
    "# random x,y positions\n",
    "X = grid_size * np.random.rand(N, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update (edit this!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def update(X, Ω, θ, λ, Ω0, A, ξ_μ, ξ_σ, grid_size):\n",
    "    \"\"\"\n",
    "    Updates positions, angles, and angular velocities of all microtubules\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    X: Nx2 matrix\n",
    "        x, y positions\n",
    "    Ω: vector of length N\n",
    "        the angular velocities\n",
    "    θ: vector of length N\n",
    "        the angles\n",
    "    λ, Ω0, A, ξ_μ, ξ_σ, grid_size:\n",
    "        parameters from above\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    Updated X, Ω, and θ\n",
    "    \"\"\"\n",
    "    D = scidist.squareform(scidist.pdist(X))\n",
    "    # boolean mask of distance matrix\n",
    "    # true if (i,j) pair of microtubules are closer than l\n",
    "    # NxN size\n",
    "    D = (D <= l)\n",
    "    \n",
    "    # count neighbors and exclude self\n",
    "    num_neighbors = (np.sum(D, 1) - 1).astype(float)\n",
    "    # 1/inf = 0 so set particles with no neighbors to infinite neighbors\n",
    "    num_neighbors[num_neighbors == 0] = np.inf\n",
    "    # num_neighbors is a vector of N elements that contains the number of \n",
    "    \n",
    "    Ω = Ω - 1 / λ * (Ω - Ω0) + np.random.normal(ξ_μ, ξ_σ, N)\n",
    "\n",
    "    θ = θ + Ω + A / num_neighbors * np.sum(2*np.sin(np.subtract.outer(θ, θ)*D), 1)\n",
    "    \n",
    "    X = np.mod(X+np.transpose([np.cos(θ), np.sin(θ)]), grid_size)\n",
    "    \n",
    "    return (X, Ω, θ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1,1);\n",
    "# initial quiver plot\n",
    "Q = ax.quiver(X[:, 0], X[:, 1], np.cos(θ), np.sin(θ),\n",
    "              pivot='mid', scale=1/(45*grid_size/N), scale_units='x', color='c')\n",
    "\n",
    "ax.set_xlim(-1, grid_size+1)\n",
    "ax.set_ylim(-1, grid_size+1)\n",
    "\n",
    "# display plot\n",
    "plt.show()\n",
    "\n",
    "for T in range(0, 10000):\n",
    "    X, Ω, θ = update(X, Ω, θ, λ, Ω0, A, ξ_μ, ξ_σ, grid_size)\n",
    "    Q.set_offsets(X)\n",
    "    Q.set_UVC(np.cos(θ), np.sin(θ))\n",
    "    plt.title(\"T={}\".format(T+1))\n",
    "    plt.draw()\n",
    "    plt.pause(0.0001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:test-env]",
   "language": "python",
   "name": "conda-env-test-env-py"
  },
  "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}