2014-01-08-Reaction_Diffusion_Accelerating_Domain.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Reaction Diffusion System on an Accelerating Domain"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In a [previous blog post](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain) we derived equations\n",
      "that describe the time behaviour of a reaction-diffusion system on a growing space domain.\n",
      "\n",
      "In that work we assumed that the velocity of domain growth is an increasing function of\n",
      "one of the two unknown variables (here, protein concentrations) described by our\n",
      "reaction diffusion system.\n",
      "\n",
      "Let us add another layer to this problem and assume that **not the growth velocity** is\n",
      "dependend on the unknown protein concentration but that **the growth acceleration** is\n",
      "a function of this unknown concentration.\n",
      "\n",
      "We derived the equations that describe the time dynamics of our previous (velocity) problem\n",
      "in material coordinates \n",
      "[here](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain#specification_of_the_local_growth_velocity)\n",
      "and reproduce them for clarity:\n",
      "\n",
      "$$\\frac{\\partial u}{\\partial t} = \\frac{D_u}{g^2} \\frac{\\partial^2 u}{\\partial X^2} - \\left( D_u \\frac{g_X}{g^3} + \\frac{a}{g} \\right) \\frac{\\partial u}{\\partial X} + f(u,v) - \\frac{a_X}{g} u,$$\n",
      "\n",
      "$$\\frac{\\partial v}{\\partial t} = \\frac{D_v}{g^2} \\frac{\\partial^2 v}{\\partial X^2} - \\left( D_v \\frac{g_X}{g^3} + \\frac{a}{g} \\right) \\frac{\\partial v}{\\partial X} - f(u,v) - \\frac{a_X}{g} v,$$\n",
      "\n",
      "$$\\frac{\\partial g}{\\partial t} = \\frac{\\partial a}{\\partial X},$$\n",
      "\n",
      "taken with initial and boundary conditions for $u$, $v$, and $g$ as\n",
      "[described earlier](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain).\n",
      "\n",
      "Let us now introduce a fourth equation to this system to describe the time behaviour of\n",
      "the growth velocity $a$ as a function of its acceleration:\n",
      "\n",
      "$$\\frac{\\partial a}{\\partial t} = s(X,t),$$\n",
      "\n",
      "where $s(X,t)$ is the local growth acceleration.\n",
      "\n",
      "In this expanded system of partial differential equations the time behaviour of the, now, unknown variable $a$\n",
      "is described by an uncoupled partial differential equation.\n",
      "Since there is no spatial coupling in this partial differential equation, we do not need to impose\n",
      "any boundary conditions.\n",
      "As initial condition, we will assume that the system is at rest at first:\n",
      "\n",
      "$$a(X,0) = 0.$$\n",
      "\n",
      "Let us now modify our \n",
      "[previous code](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain#python_code)\n",
      "to simulate this expanded system.\n",
      "\n",
      "Most of this code is just a copy-and-paste of our\n",
      "[previous code](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain#python_code)\n",
      "so please refer to our comments there if anything is unclear."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "from scipy import integrate\n",
      "from matplotlib import pyplot\n",
      "numpy.set_printoptions(precision=3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L = 1.\n",
      "J = 100\n",
      "dX = float(L)/float(J)\n",
      "X_grid = numpy.array([float(dX)/2. + j*dX for j in range(J)])\n",
      "x_grid = numpy.array([j*dX for j in range(J+1)])\n",
      "\n",
      "T = 200\n",
      "N = 1000\n",
      "dt = float(T)/float(N-1)\n",
      "t_grid = numpy.array([n*dt for n in range(N)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "D_v = float(10.)/float(100.)\n",
      "D_u = 0.01 * D_v"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "k0 = 0.067\n",
      "f_vec = lambda U, V: numpy.multiply(dt, numpy.subtract(numpy.multiply(V, \n",
      "                     numpy.add(k0, numpy.divide(numpy.multiply(U,U), numpy.add(1., numpy.multiply(U,U))))), U))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The initial condition for the growth velocity `a` is set to zero everywhere so that the space domain is at rest at first.\n",
      "\n",
      "The expression for the growth acceleration `s` is the exact same expression we used\n",
      "[earlier](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain#python_code)\n",
      "for the growth velocity."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "total_protein = 2.26\n",
      "\n",
      "no_high = 10\n",
      "U = numpy.array([0.1 for i in range(no_high,J)] + [2. for i in range(0,no_high)])\n",
      "V = numpy.array([float(total_protein-dX*sum(U))/float(J*dX) for i in range(0,J)])\n",
      "g = numpy.array([1. for j in range(J)])\n",
      "a = numpy.array([0. for j in range(J)])\n",
      "s = lambda U: numpy.array([0.001*X_grid[j]*U[j] for j in range(J)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is the *polarized* initial condition we choose for our concentration system `U` and `V`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pyplot.ylim((0., 2.1))\n",
      "pyplot.xlabel('X'); pyplot.ylabel('concentration')\n",
      "pyplot.plot(X_grid, U)\n",
      "pyplot.plot(X_grid, V)\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As before, the initial condition for the slope `g` of the trajectories $G$ is `1.` everywhere:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pyplot.plot(X_grid, g)\n",
      "pyplot.xlabel('X'); pyplot.ylabel('g')\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Initially, `s` looks the same as `a` did in our [previous code](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain#python_code),\n",
      "while in this version of our code `a` is set initially to `0.` everywhere:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pyplot.plot(X_grid, a)\n",
      "pyplot.plot(X_grid, s(U))\n",
      "pyplot.xlabel('X'); pyplot.ylabel('g')\n",
      "pyplot.ylim((-.0001, 0.0021))\n",
      "pyplot.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# syntax to get one element in an array: http://stackoverflow.com/a/7332880/3078529\n",
      "a_left = lambda a: numpy.concatenate((a[1:J], a[J-1:J]))\n",
      "a_right = lambda a: numpy.concatenate((a[0:1], a[0:J-1]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "f_vec_u = lambda U, V, g, a: numpy.subtract(f_vec(U, V),\n",
      "                                         numpy.multiply(numpy.subtract(a_left(a), a_right(a)),\n",
      "                                                        numpy.multiply(float(dt)/(float(2.*dX)),\n",
      "                                                                       numpy.divide(U,g))))\n",
      "\n",
      "f_vec_v = lambda U, V, g, a: numpy.subtract(numpy.multiply(-1., f_vec(U, V)),\n",
      "                                         numpy.multiply(numpy.subtract(a_left(a), a_right(a)),\n",
      "                                                        numpy.multiply(float(dt)/(float(2.*dX)),\n",
      "                                                                       numpy.divide(V,g))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sigma_u_func = lambda g: numpy.divide(float(D_u*dt)/float(2.*dX*dX),\n",
      "                                      numpy.multiply(g, g))\n",
      "sigma_v_func = lambda g: numpy.divide(float(D_v*dt)/float(2.*dX*dX),\n",
      "                                      numpy.multiply(g, g))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# syntax to get one element in an array: http://stackoverflow.com/a/7332880/3078529\n",
      "g_left = lambda g: numpy.concatenate((g[1:J], g[J-1:J]))\n",
      "g_right = lambda g: numpy.concatenate((g[0:1], g[0:J-1]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "rho_u_func = lambda g, a: numpy.multiply(float(-dt)/float(4.*dX),\n",
      "                                      numpy.add(numpy.divide(a, g),\n",
      "                                                numpy.multiply(numpy.subtract(g_left(g), g_right(g)),\n",
      "                                                               numpy.divide(float(D_u)/(2.*dX),\n",
      "                                                                            numpy.power(g, 3)))))\n",
      "rho_v_func = lambda g, a: numpy.multiply(float(-dt)/float(4.*dX),\n",
      "                                      numpy.add(numpy.divide(a, g),\n",
      "                                                numpy.multiply(numpy.subtract(g_left(g), g_right(g)),\n",
      "                                                               numpy.divide(float(D_v)/(2.*dX),\n",
      "                                                                            numpy.power(g, 3)))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that here we use slightly more meaningful variable names for the ODE steppers for `g` and `x_grid` than\n",
      "in our [previous code](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain#python_code).\n",
      "\n",
      "The ODE stepper for `a` is of the same form except that the right-hand side is governed by the\n",
      "growth acceleration `s`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "g_rhs = lambda t, g, a: numpy.divide(numpy.subtract(a_left(a), a_right(a)),\n",
      "                                     numpy.array([dX]+[2.*dX for j in range(J-2)]+[dX]))\n",
      "\n",
      "g_stepper = integrate.ode(g_rhs)\n",
      "g_stepper = g_stepper.set_integrator('dopri5', nsteps=10, max_step=dt)\n",
      "g_stepper = g_stepper.set_initial_value(g, 0.)\n",
      "g_stepper = g_stepper.set_f_params(a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a_rhs = lambda t, a, s: s\n",
      "\n",
      "a_stepper = integrate.ode(a_rhs)\n",
      "a_stepper = a_stepper.set_integrator('dopri5', nsteps=10, max_step=dt)\n",
      "a_stepper = a_stepper.set_initial_value(a, 0.)\n",
      "a_stepper = a_stepper.set_f_params(s(U))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x_grid_rhs = lambda t, x_grid, a: numpy.concatenate(([0.],\n",
      "                                                     numpy.divide(numpy.add(a[0:J-1], a[1:J]), 2.),\n",
      "                                                     a[J-1:J]))\n",
      "\n",
      "x_stepper = integrate.ode(x_grid_rhs)\n",
      "x_stepper = x_stepper.set_integrator('dopri5', nsteps=10, max_step=dt)\n",
      "x_stepper = x_stepper.set_initial_value(x_grid, 0.)\n",
      "x_stepper = x_stepper.set_f_params(a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U_record = []\n",
      "V_record = []\n",
      "g_record = []\n",
      "a_record = []\n",
      "x_record = []\n",
      "\n",
      "U_record.append(U)\n",
      "V_record.append(V)\n",
      "g_record.append(g)\n",
      "a_record.append(a)\n",
      "x_record.append(x_grid)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us now integrate this system for a total of `N` time points.\n",
      "\n",
      "Notice that we update the ODE stepper of `g`, `a` and `x_grid` every time step with the current state\n",
      "of their corresponding right-hand-side expressions.\n",
      "This is done with a call to `.set_f_params()` on the corresponding objects."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for ti in range(1,N):\n",
      "    sigma_u = sigma_u_func(g)\n",
      "    sigma_v = sigma_v_func(g)\n",
      "    rho_u = rho_u_func(g, a)\n",
      "    rho_v = rho_v_func(g, a)\n",
      "    \n",
      "    A_u = numpy.diagflat([-sigma_u[j]+rho_u[j] for j in range(1,J)], -1) +\\\n",
      "          numpy.diagflat([1.+sigma_u[0]+rho_u[0]]+[1.+2.*sigma_u[j] for j in range(1,J-1)]+[1.+sigma_u[J-1]-rho_u[J-1]]) +\\\n",
      "          numpy.diagflat([-(sigma_u[j]+rho_u[j]) for j in range(0,J-1)], 1)\n",
      "        \n",
      "    B_u = numpy.diagflat([sigma_u[j]-rho_u[j] for j in range(1,J)], -1) +\\\n",
      "          numpy.diagflat([1.-sigma_u[0]-rho_u[0]]+[1.-2.*sigma_u[j] for j in range(1,J-1)]+[1.-sigma_u[J-1]+rho_u[J-1]]) +\\\n",
      "          numpy.diagflat([sigma_u[j]+rho_u[j] for j in range(0,J-1)], 1)\n",
      "        \n",
      "    A_v = numpy.diagflat([-sigma_v[j]+rho_v[j] for j in range(1,J)], -1) +\\\n",
      "          numpy.diagflat([1.+sigma_v[0]+rho_v[0]]+[1.+2.*sigma_v[j] for j in range(1,J-1)]+[1.+sigma_v[J-1]-rho_v[J-1]]) +\\\n",
      "          numpy.diagflat([-(sigma_v[j]+rho_v[j]) for j in range(0,J-1)], 1)\n",
      "        \n",
      "    B_v = numpy.diagflat([sigma_v[j]-rho_v[j] for j in range(1,J)], -1) +\\\n",
      "          numpy.diagflat([1.-sigma_v[0]-rho_v[0]]+[1.-2.*sigma_v[j] for j in range(1,J-1)]+[1.-sigma_v[J-1]+rho_v[J-1]]) +\\\n",
      "          numpy.diagflat([sigma_v[j]+rho_v[j] for j in range(0,J-1)], 1)\n",
      "            \n",
      "    U_new = numpy.linalg.solve(A_u, B_u.dot(U) + f_vec_u(U, V, g, a))\n",
      "    V_new = numpy.linalg.solve(A_v, B_v.dot(V) + f_vec_v(U, V, g, a))\n",
      "    \n",
      "    while a_stepper.successful() and a_stepper.t + dt < ti*dt:\n",
      "        a_stepper.integrate(a_stepper.t + dt)\n",
      "    while g_stepper.successful() and g_stepper.t + dt < ti*dt:\n",
      "        g_stepper.integrate(g_stepper.t + dt)\n",
      "    while x_stepper.successful() and x_stepper.t + dt < ti*dt:\n",
      "        x_stepper.integrate(x_stepper.t + dt)\n",
      "        \n",
      "    g_stepper = g_stepper.set_f_params(a)\n",
      "    x_stepper = x_stepper.set_f_params(a)\n",
      "    a_stepper = a_stepper.set_f_params(s(U))\n",
      "    \n",
      "    U = U_new\n",
      "    V = V_new\n",
      "    \n",
      "    # these are the correct \"y\" values to save for the current time step since\n",
      "    # we integrate only up to t < ti*dt\n",
      "    g = g_stepper.y\n",
      "    a = a_stepper.y\n",
      "    x_grid = x_stepper.y\n",
      "    \n",
      "    U_record.append(U)\n",
      "    V_record.append(V)\n",
      "    g_record.append(g)\n",
      "    a_record.append(a)\n",
      "    x_record.append(x_grid)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As in our\n",
      "[previous code](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain#python_code)\n",
      "protein mass is conserved by definition of our boundary conditions for `U` and `V`.\n",
      "\n",
      "To make certain that our numerical integration does conserve protein mass (at least to\n",
      "some sensible degree) we work out the initial total mass and final total mass respectively:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum(numpy.multiply(numpy.diff(x_record[0]),U_record[0]) + numpy.multiply(numpy.diff(x_record[0]),V_record[0]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sum(numpy.multiply(numpy.diff(x_record[-1]),U_record[-1]) + numpy.multiply(numpy.diff(x_record[-1]),V_record[-1]))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These numbers look good and we therefore assume that our numerical integration is somewhat stable."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A look at a kymograph of the concentration of `U` reveals that the initial condition on `U` is\n",
      "lost quickly - which will be due to reasons we discussed\n",
      "[earlier](http://georg.io/2013/12/Reaction_Diffusion_Growing_Domain)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U_record = numpy.array(U_record)\n",
      "V_record = numpy.array(V_record)\n",
      "\n",
      "fig, ax = pyplot.subplots()\n",
      "pyplot.xlabel('X'); pyplot.ylabel('t')\n",
      "heatmap = ax.pcolormesh(x_record[0], t_grid, U_record, vmin=0., vmax=1.2)\n",
      "colorbar = pyplot.colorbar(heatmap)\n",
      "colorbar.set_label('concentration U')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since with the above setup **we only introduce acceleration and no deceleration**, we are not surprised\n",
      "to observe that the growth velocity simply increases over time:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a_record = numpy.array(a_record)\n",
      "\n",
      "fig, ax = pyplot.subplots()\n",
      "pyplot.xlabel('X'); pyplot.ylabel('t')\n",
      "heatmap = ax.pcolormesh(X_grid, t_grid, a_record)\n",
      "colorbar = pyplot.colorbar(heatmap)\n",
      "colorbar.set_label('local growth velocity a(X,t)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This continual increase in growth velocity is only due to an initial\n",
      "acceleration due to the initial condition we imposed on `U`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = pyplot.subplots()\n",
      "pyplot.xlabel('X'); pyplot.ylabel('t')\n",
      "s_record = numpy.empty(shape=(N, J))\n",
      "for ti in range(N):\n",
      "    s_record[ti] = s(U_record[ti])\n",
      "heatmap = ax.pcolormesh(X_grid, t_grid, s_record)\n",
      "colorbar = pyplot.colorbar(heatmap)\n",
      "colorbar.set_label('local growth acceleration s')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}