{ "metadata": { "name": "", "signature": "sha256:1ce33c8e48cfac71628e86512053691932a544c53e70c1fe3c46fc60cd380586" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "heading", "level": 1, "metadata": {}, "source": [ " 1D Turing patterns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*This document was generated from an IPython notebook. You can download the notebook [here](turing.ipynb). The module you need for the calculation, `rd_solver.py`, may be downloaded [here](rd_solver.py).*" ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Turing patterns from two reactants" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The reaction-diffusion solver solves the following system of 1-D PDEs describing the dynamics of the concentrations of two species, $a$ and $b$.\n", "\n", "\\begin{align}\n", "\\frac{\\partial a}{\\partial t} &= D_a \\,\\frac{\\partial^2a}{\\partial x^2} + r_a \\\\\n", "\\frac{\\partial b}{\\partial t} &= D_b \\,\\frac{\\partial^2b}{\\partial x^2} + r_b,\n", "\\end{align}\n", "\n", "where $r_a$ and $r_b$ are respectively the rates of production of $a$ and $b$ by chemical reaction. While the user can specify arbitrary $r_a$ and $r_b$, two classic models are included. They are the Gierer-Meinhardt (G-M) and the activator-substrate depletion model (ASDM)." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Gierer-Meinhardt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the Gierer-Meinhardt system, we call $a$ the \"activator\" and the second species $h$ the \"inhibitor.\" The reaction terms are\n", "\n", "\\begin{align}\n", "r_a &= \\rho_a\\,\\frac{a^2}{h} - \\mu_a a \\\\[2mm]\n", "r_h &= \\rho_h a^2 - \\mu_h h.\n", "\\end{align}\n", "\n", "When doing numerical calculations, it is always a good idea to work in dimensionless units. It is possible to nondimensionalize the governing equations to give\n", "\n", "\\begin{align}\n", "\\frac{\\partial a}{\\partial t} &= d\\, \\frac{\\partial^2a}{\\partial x^2} +\\frac{a^2}{h} - a \\\\[2mm]\n", "\\frac{\\partial h}{\\partial t} &= \\frac{\\partial^2h}{\\partial x^2} + \\rho(a^2 - h),\n", "\\end{align}\n", "\n", "where $a$, $h$, are now properly rescaled and $d = D_a / D_h$. We therefore only have two parameters, $d$ and $\\rho$ to vary.\n", "\n", "The Gierer-Meinhardt system has a unique homogenious steady state of $a_0 = h_0 = 1$." ] }, { "cell_type": "heading", "level": 3, "metadata": {}, "source": [ "Activator-substrate depletion model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the ASDM, we call $a$ the \"activator\" and the second species $s$ the \"substrate.\" The reaction terms are\n", "\n", "\\begin{align}\n", "r_a &= \\rho_a a^2 s - \\mu_a a \\\\[2mm]\n", "r_s &= -\\rho_s a^2 s + \\sigma_s.\n", "\\end{align}\n", "\n", "Upon nondimensionalizing, we get\n", "\n", "\\begin{align}\n", "\\frac{\\partial a}{\\partial t} &= d\\, \\frac{\\partial^2a}{\\partial x^2} + a^2s - a \\\\[2mm]\n", "\\frac{\\partial h}{\\partial t} &= \\frac{\\partial^2h}{\\partial x^2} + \\mu(1 - a^2 s),\n", "\\end{align}\n", "\n", "giving the two parameters $\\mu$ and $d = D_a/D_s$.\n", "\n", "The ASDM has a unique homogeneous steady state of $a_0 = s_0 = 1$." ] }, { "cell_type": "heading", "level": 2, "metadata": {}, "source": [ "Performing a calculation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We step through an example calculation using the Gierer-Meinhardt system. We start by importing the `rd_solver` module, along with `numpy`, whose utilities we will use. We also import some `matplotlib` tools for plotting." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from matplotlib import animation\n", "%matplotlib inline\n", "\n", "import rd_solver as rd" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we specify the function describing the chemical reactions. This is the pre-made `gierer-meinhardt_rxn` function. You can write your own with arbitrary parameters, if you please. They just have to have the functional form `rxn_name(a, b, *params)`." ] }, { "cell_type": "code", "collapsed": false, "input": [ "rxn_fun = rd.gierer_meinhardt_rxn" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to define our parameters for the calculation. We start by defining our initial condition. This will also define the number of grid points we will use. We will make a small random perturbation about the homogeneous steady state of $a_0 = h_0 = 1$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Set up steady state (using 500 grid points)\n", "a_0 = np.ones(500)\n", "h_0 = np.ones(500)\n", "\n", "# Seed the random number generator for reproducibility\n", "np.random.seed(42)\n", "\n", "# Make a small perturbation to a_0\n", "a_0 += 0.001 * np.random.rand(500)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now define our system parameters. Specifically, we choose the physical length of our system, the time points we want to get results for, and the parameters $d$ and $\\rho$. We have to package them properly for what the functions requires. We will use $d = 0.05$ and $\\rho = 1.4$ on a system of length $20$. Note that all of these parameters are dimensionless." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Time points\n", "t = np.linspace(0.0, 100.0, 100)\n", "\n", "# Physical length of system\n", "len_x = 20.0\n", "\n", "# Diffusion coefficient\n", "diff_coeff = (0.05, 1.0)\n", "\n", "# Reaction parameter (must be a tuple of params, even though only 1 for G-M)\n", "rxn_params = (1.4,)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 21 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are now ready to solve the system of PDEs! Under the hood, the solver used the method of lines. It does second order central differencing to compute the derivatives, and then uses `scipy.odeint` to do the time stepping." ] }, { "cell_type": "code", "collapsed": false, "input": [ "a, h = rd.solve(a_0, h_0, t, len_x, diff_coeff, rxn_fun, rxn_params, \n", " banded=True)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 22 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now plot the final profile." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Didn't need x to compute, need it to plot (remember: it's dimensionless)\n", "x = np.linspace(0.0, len_x, len(a_0))\n", "\n", "# Plot profiles\n", "plt.plot(x, a[-1,:], '-', color='#1F78B4', lw=2, label='activator')\n", "plt.plot(x, h[-1,:], '-', color='#33A02C', lw=2, label='inhibitor')\n", "plt.xlabel('$x$', fontsize=18)\n", "plt.ylabel('concentration', fontsize=18)\n", "plt.legend(framealpha=0.5)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 23, "text": [ "" ] }, { "metadata": {}, "output_type": "display_data", "svg": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text": [ "" ] } ], "prompt_number": 23 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have all the time points, so we can animate the solution if we like. There are lots of ways to do this, such as looping through and making plots, saving as PNGs, and then stiching them together into a movie with some other program. We will use some of `matplotlib`'s nifty animators, which rely on `ffmpeg` under the hood.\n", "\n", "The embedded video may not work with some browsers. I tested it with Chrome 38 and Safari 8." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# Import modules for encoding and embedding movie\n", "from IPython.display import HTML\n", "from base64 import b64encode\n", "\n", "# Get bounds on results\n", "max_conc = max(a.max(), h.max())\n", "\n", "# Set up figure and set axis bounds\n", "fig = plt.figure()\n", "ax = plt.axes(xlim=(x[0], x[-1]), ylim=(0, max_conc * 1.02))\n", "lines = [ax.plot([], [], lw=2, color='#1F78B4', label='activator')[0],\n", " ax.plot([], [], lw=2, color='#33A02C', label='inhibitor')[0]]\n", "ax.set_xlabel('$x$', fontsize=18)\n", "ax.set_ylabel('concentration', fontsize=18)\n", "ax.legend(framealpha=0.5)\n", "\n", "# Initialize plotting\n", "def init():\n", " lines[0].set_data([],[])\n", " lines[1].set_data([],[])\n", " return lines\n", "\n", "# This paints in the animation\n", "def animate(i):\n", " lines[0].set_data(x, a[i,:])\n", " lines[1].set_data(x, h[i,:])\n", " return lines\n", "\n", "# call the animator.\n", "anim = animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=len(t), interval=20, blit=True)\n", "\n", "# Save the animation\n", "anim.save('gierer_meinhardt.mp4', fps=10, dpi=100, \n", " extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])\n", "\n", "# Close the figure\n", "plt.close()\n", "\n", "# Encode/embed the animation\n", "video = open(\"gierer_meinhardt.mp4\", \"rb\").read()\n", "video_encoded = b64encode(video)\n", "video_tag = \"\"\"\n", "\n", "\"\"\".format(video_encoded)\n", "HTML(data=video_tag)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "\n", "\n" ], "metadata": {}, "output_type": "pyout", "prompt_number": 24, "text": [ "" ] } ], "prompt_number": 24 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can also plot the result as kymographs. We will use the blue channel for activator and the green channel for inhibitor." ] }, { "cell_type": "code", "collapsed": false, "input": [ "# RGB values for kymgraph image\n", "rgb_a = np.zeros((len(t), len(x), 3))\n", "rgb_h = np.zeros((len(t), len(x), 3))\n", "\n", "# Values for the activator\n", "rgb_a[:,:,2] = a / a.max()\n", "\n", "# Values for the inhibitor\n", "rgb_h[:,:,1] = h / h.max()\n", "\n", "# Make kymograph\n", "fig, ax = plt.subplots(1, 2, sharey=True)\n", "ax[0].imshow(rgb_a, extent=(t[0], t[-1]/2, t[-1], t[0]))\n", "ax[1].imshow(rgb_h, extent=(t[0], t[-1]/2, t[-1], t[0]))\n", "ax[0].set_ylabel('t (s)')\n", "ax[0].set_title('activator')\n", "ax[1].set_title('inhibitor')\n", "for i in [0, 1]:\n", " ax[i].set_xlim((x[0], x[-1]))\n", " ax[i].set_xticks([t[0], t[-1]/2])\n", " ax[i].set_xticklabels([x[0], x[-1]])\n", " ax[i].set_xlabel(u'x (\u00b5m)')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "svg": [ "\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text": [ "" ] } ], "prompt_number": 44 } ], "metadata": {} } ] }