Skip to content

Instantly share code, notes, and snippets.

@tlnagy
Last active September 7, 2017 19:29
Show Gist options
  • Save tlnagy/cba938ffd5c98236e90bfd1dc3d23d11 to your computer and use it in GitHub Desktop.
Save tlnagy/cba938ffd5c98236e90bfd1dc3d23d11 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"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
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment