![](\"caltech_rosen.png\")
\n",
"\n",
"\n",
"\\n\"+\n", " \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n", " \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n", " \"
\\n\"+\n", " \"\\n\"+\n",
" \"from bokeh.resources import INLINE\\n\"+\n",
" \"output_notebook(resources=INLINE)\\n\"+\n",
" \"\\n\"+\n",
" \"\\n\"+\n \"BokehJS does not appear to have successfully loaded. If loading BokehJS from CDN, this \\n\"+\n \"may be due to a slow or bad network connection. Possible fixes:\\n\"+\n \"
\\n\"+\n \"\\n\"+\n \"from bokeh.resources import INLINE\\n\"+\n \"output_notebook(resources=INLINE)\\n\"+\n \"\\n\"+\n \"![](\"expander.png\")
\n",
"\n",
"The morphogen and expander are thus in a negative feedback loop. How might this help us with the issue that the morphogen gradient should scale?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Numerical solution to R-D equations\n",
"\n",
"To solve R-D systems, I wrote the two functions below, which are the main engines for computation. The first computes the time derivative of concentrations at each of the spatial grid points. The second is a wrapper around `scipy.integrate.odeint()` to do the time stepping. Their doc strings describe them in more detail."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def dc_dt(c, t, x, derivs_0, derivs_L, diff_coeff_fun, diff_coeff_params, \n",
" rxn_fun, rxn_params, n_species, h):\n",
" \"\"\"\n",
" Time derivative of concentrations in an R-D system \n",
" for constant flux BCs.\n",
" \n",
" Parameters\n",
" ----------\n",
" c : ndarray, shape (n_species * n_gridpoints)\n",
" The concentration of the chemical species interleaved in a\n",
" a NumPy array. The interleaving allows us to take advantage \n",
" of the banded structure of the Jacobian when using the \n",
" Hindmarsh algorithm for integrating in time.\n",
" t : float\n",
" Time.\n",
" derivs_0 : ndarray, shape (n_species)\n",
" derivs_0[i] is the value of the diffusive flux,\n",
" D dc_i/dx, at x = 0, the leftmost boundary of the domain of x.\n",
" derivs_L : ndarray, shape (n_species)\n",
" derivs_0[i] is the value of the diffusive flux,\n",
" D dc_i/dx, at x = L, the rightmost boundary of the domain of x.\n",
" diff_coeff_fun : function\n",
" Function of the form diff_coeff_fun(c_tuple, t, x, *diff_coeff_params).\n",
" Returns an tuple where entry i is a NumPy array containing\n",
" the diffusion coefficient of species i at the grid points.\n",
" c_tuple[i] is a NumPy array containing the concentrations of\n",
" species i at the grid poitns.\n",
" diff_coeff_params : arbitrary\n",
" Tuple of parameters to be passed into diff_coeff_fun.\n",
" rxn_fun : function\n",
" Function of the form rxn_fun(c_tuple, t, *rxn_params).\n",
" Returns an tuple where entry i is a NumPy array containing\n",
" the net rate of production of species i by chemical reaction\n",
" at the grid points. c_tuple[i] is a NumPy array containing \n",
" the concentrations of species i at the grid poitns.\n",
" rxn_params : arbitrary\n",
" Tuple of parameters to be passed into rxn_fun.\n",
" n_species : int\n",
" Number of chemical species.\n",
" h : float\n",
" Grid spacing (assumed to be constant)\n",
" \n",
" Returns\n",
" -------\n",
" dc_dt : ndarray, shape (n_species * n_gridpoints)\n",
" The time derivatives of the concentrations of the chemical\n",
" species at the grid points interleaved in a NumPy array.\n",
" \"\"\"\n",
" # Tuple of concentrations\n",
" c_tuple = tuple([c[i::n_species] for i in range(n_species)])\n",
" \n",
" # Compute diffusion coefficients\n",
" D_tuple = diff_coeff_fun(c_tuple, t, x, *diff_coeff_params)\n",
" \n",
" # Compute reaction terms\n",
" rxn_tuple = rxn_fun(c_tuple, t, *rxn_params)\n",
"\n",
" # Return array\n",
" conc_deriv = np.empty_like(c)\n",
" \n",
" # Convenient array for storing concentrations\n",
" da_dt = np.empty(len(c_tuple[0]))\n",
"\n",
" # Useful to have square of grid spacing around\n",
" h2 = h**2\n",
" \n",
" # Compute diffusion terms (central differencing w/ Neumann BCs)\n",
" for i in range(n_species):\n",
" # View of concentrations and diffusion coeff. for convenience\n",
" a = np.copy(c_tuple[i])\n",
" D = np.copy(D_tuple[i])\n",
" \n",
" # Time derivative at left boundary\n",
" da_dt[0] = D[0] / h2 * 2 * (a[1] - a[0] - h * derivs_0[i])\n",
" \n",
" # First derivatives of D and a\n",
" dD_dx = (D[2:] - D[:-2]) / (2*h)\n",
" da_dx = (a[2:] - a[:-2]) / (2*h)\n",
"\n",
" # Time derivative for middle grid points\n",
" da_dt[1:-1] = D[1:-1] * np.diff(a, 2) / h2 + dD_dx*da_dx\n",
" \n",
" # Time derivative at left boundary\n",
" da_dt[-1] = D[-1] / h2 * 2 * (a[-2] - a[-1] + h * derivs_L[i])\n",
" \n",
" # Store in output array with reaction terms\n",
" conc_deriv[i::n_species] = da_dt + rxn_tuple[i]\n",
" \n",
" return conc_deriv\n",
"\n",
"\n",
"def rd_solve(c_0_tuple, t, L=1, derivs_0=0, derivs_L=0, diff_coeff_fun=None, \n",
" diff_coeff_params=(), rxn_fun=None, rxn_params=(),\n",
" rtol=1.49012e-8, atol=1.49012e-8):\n",
" \"\"\"\n",
" Parameters\n",
" ----------\n",
" c_0_tuple : tuple\n",
" c_0_tuple[i] is a NumPy array of length n_gridpoints with the \n",
" initial concentrations of chemical species i at the grid points.\n",
" t : ndarray\n",
" An array of time points for which the solution is desired.\n",
" L : float\n",
" Total length of the x-domain.\n",
" derivs_0 : ndarray, shape (n_species)\n",
" derivs_0[i] is the value of dc_i/dx at x = 0.\n",
" derivs_L : ndarray, shape (n_species)\n",
" derivs_L[i] is the value of dc_i/dx at x = L, the rightmost\n",
" boundary of the domain of x.\n",
" diff_coeff_fun : function\n",
" Function of the form diff_coeff_fun(c_tuple, t, *diff_coeff_params).\n",
" Returns an tuple where entry i is a NumPy array containing\n",
" the diffusion coefficient of species i at the grid points.\n",
" c_tuple[i] is a NumPy array containing the concentrations of\n",
" species i at the grid poitns.\n",
" diff_coeff_params : arbitrary\n",
" Tuple of parameters to be passed into diff_coeff_fun.\n",
" rxn_fun : function\n",
" Function of the form rxn_fun(c_tuple, t, *rxn_params).\n",
" Returns an tuple where entry i is a NumPy array containing\n",
" the net rate of production of species i by chemical reaction\n",
" at the grid points. c_tuple[i] is a NumPy array containing \n",
" the concentrations of species i at the grid poitns.\n",
" rxn_params : arbitrary\n",
" Tuple of parameters to be passed into rxn_fun.\n",
" rtol : float\n",
" Relative tolerance for solver. Default os odeint's default.\n",
" atol : float\n",
" Absolute tolerance for solver. Default os odeint's default.\n",
" \n",
" Returns\n",
" -------\n",
" c_tuple : tuple\n",
" c_tuple[i] is a NumPy array of shape (len(t), n_gridpoints)\n",
" with the initial concentrations of chemical species i at \n",
" the grid points over time.\n",
" \n",
" Notes\n",
" -----\n",
" .. When intergrating for long times near a steady state, you\n",
" may need to lower the absolute tolerance (atol) because the\n",
" solution does not change much over time and it may be difficult\n",
" for the solver to maintain tight tolerances.\n",
" \"\"\"\n",
" # Number of grid points\n",
" n_gridpoints = len(c_0_tuple[0])\n",
" \n",
" # Number of chemical species\n",
" n_species = len(c_0_tuple)\n",
"\n",
" # Grid spacing\n",
" h = (L / (n_gridpoints - 1))\n",
" \n",
" # Grid points\n",
" x = np.linspace(0, L, n_gridpoints)\n",
" \n",
" # Set up boundary conditions\n",
" if np.isscalar(derivs_0):\n",
" derivs_0 = np.array(n_species * [derivs_0])\n",
" if np.isscalar(derivs_L):\n",
" derivs_L = np.array(n_species * [derivs_L])\n",
"\n",
" # Set up parameters to be passed in to dc_dt\n",
" params = (x, derivs_0, derivs_L, diff_coeff_fun, diff_coeff_params, \n",
" rxn_fun, rxn_params, n_species, h)\n",
"\n",
" # Set up initial condition\n",
" c0 = np.empty(n_species * n_gridpoints)\n",
" for i in range(n_species):\n",
" c0[i::n_species] = c_0_tuple[i]\n",
"\n",
" # Solve using odeint, taking advantage of banded structure\n",
" c = scipy.integrate.odeint(\n",
" dc_dt, c0, t, args=params, ml=n_species, mu=n_species,\n",
" rtol=rtol, atol=atol)\n",
" \n",
" return tuple([c[:,i::n_species] for i in range(n_species)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### A sample calculation: the activator-substrate depletion model (ASDM)\n",
"\n",
"We will use the R-D solver to solve the ASDM model, which exhibits Turing patterns. 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$.\n",
"\n",
"We will now write our functions to compute the diffusion coefficients and reaction terms for the ASDM, using the required I/O formats for the functions."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def constant_diff_coeffs(c_tuple, t, x, diff_coeffs):\n",
" n = len(c_tuple[0])\n",
" return tuple([diff_coeffs[i] * np.ones(n) for i in range(len(c_tuple))])\n",
"\n",
"\n",
"def asdm_rxn(as_tuple, t, mu):\n",
" \"\"\"\n",
" Reaction expression for activator-substrate depletion model.\n",
"\n",
" Returns the rate of production of activator and substrate, respectively.\n",
"\n",
" r_a = a**2 * s - a\n",
" r_s = mu * (1 - a**2 * s)\n",
" \"\"\"\n",
" # Unpack concentrations\n",
" a, s = as_tuple\n",
" \n",
" # Compute and return reaction rates\n",
" a2s = a**2 * s\n",
" return (a2s - a, mu * (1.0 - a2s))\n",
"\n",
"\n",
"# Set up steady state (using 500 grid points)\n",
"a_0 = np.ones(500)\n",
"s_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(len(a_0))\n",
"\n",
"# Time points\n",
"t = np.linspace(0.0, 100.0, 100)\n",
"\n",
"# Physical length of system\n",
"L = 20.0\n",
"\n",
"# Diffusion coefficient\n",
"diff_coeffs = (0.05, 1.0)\n",
"\n",
"# Reaction parameter (must be a tuple of params, even though only 1 for ASDM)\n",
"rxn_params = (1.4,)\n",
"\n",
"# Solve\n",
"a, s = rd_solve((a_0, s_0), t, L=L, derivs_0=0, derivs_L=0, \n",
" diff_coeff_fun=constant_diff_coeffs, \n",
" diff_coeff_params=(diff_coeffs,), rxn_fun=asdm_rxn, \n",
" rxn_params=rxn_params)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that we've done the solution, let's plot the result with a slider so that we can change the values of the time points."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"