The reaction-diffusion solver solves the following system of 1-D PDEs describing the dynamics of the concentrations of two species, $a$ and $b$.
\begin{align} \frac{\partial a}{\partial t} &= D_a \,\frac{\partial^2a}{\partial x^2} + r_a \\ \frac{\partial b}{\partial t} &= D_b \,\frac{\partial^2b}{\partial x^2} + r_b, \end{align}
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).
For the Gierer-Meinhardt system, we call $a$ the "activator" and the second species $h$ the "inhibitor." The reaction terms are
\begin{align} r_a &= \rho_a\,\frac{a^2}{h} - \mu_a a \\[2mm] r_h &= \rho_h a^2 - \mu_h h. \end{align}
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
\begin{align} \frac{\partial a}{\partial t} &= d\, \frac{\partial^2a}{\partial x^2} +\frac{a^2}{h} - a \\[2mm] \frac{\partial h}{\partial t} &= \frac{\partial^2h}{\partial x^2} + \rho(a^2 - h), \end{align}
where $a$, $h$, are now properly rescaled and $d = D_a / D_h$. We therefore only have two parameters, $d$ and $\rho$ to vary.
The Gierer-Meinhardt system has a unique homogenious steady state of $a_0 = h_0 = 1$.
For the ASDM, we call $a$ the "activator" and the second species $s$ the "substrate." The reaction terms are
\begin{align} r_a &= \rho_a a^2 s - \mu_a a \\[2mm] r_s &= -\rho_s a^2 s + \sigma_s. \end{align}
Upon nondimensionalizing, we get
\begin{align} \frac{\partial a}{\partial t} &= d\, \frac{\partial^2a}{\partial x^2} + a^2s - a \\[2mm] \frac{\partial h}{\partial t} &= \frac{\partial^2h}{\partial x^2} + \mu(1 - a^2 s), \end{align}
giving the two parameters $\mu$ and $d = D_a/D_s$.
The ASDM has a unique homogeneous steady state of $a_0 = s_0 = 1$.
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.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
%matplotlib inline
import rd_solver as rd
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).
rxn_fun = rd.gierer_meinhardt_rxn
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$.
# Set up steady state (using 500 grid points)
a_0 = np.ones(500)
h_0 = np.ones(500)
# Seed the random number generator for reproducibility
np.random.seed(42)
# Make a small perturbation to a_0
a_0 += 0.001 * np.random.rand(500)
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.
# Time points
t = np.linspace(0.0, 100.0, 100)
# Physical length of system
len_x = 20.0
# Diffusion coefficient
diff_coeff = (0.05, 1.0)
# Reaction parameter (must be a tuple of params, even though only 1 for G-M)
rxn_params = (1.4,)
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.
a, h = rd.solve(a_0, h_0, t, len_x, diff_coeff, rxn_fun, rxn_params,
banded=True)
We can now plot the final profile.
# Didn't need x to compute, need it to plot (remember: it's dimensionless)
x = np.linspace(0.0, len_x, len(a_0))
# Plot profiles
plt.plot(x, a[-1,:], '-', color='#1F78B4', lw=2, label='activator')
plt.plot(x, h[-1,:], '-', color='#33A02C', lw=2, label='inhibitor')
plt.xlabel('$x$', fontsize=18)
plt.ylabel('concentration', fontsize=18)
plt.legend(framealpha=0.5)
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.
The embedded video may not work with some browsers. I tested it with Chrome 38 and Safari 8.
# Import modules for encoding and embedding movie
from IPython.display import HTML
from base64 import b64encode
# Get bounds on results
max_conc = max(a.max(), h.max())
# Set up figure and set axis bounds
fig = plt.figure()
ax = plt.axes(xlim=(x[0], x[-1]), ylim=(0, max_conc * 1.02))
lines = [ax.plot([], [], lw=2, color='#1F78B4', label='activator')[0],
ax.plot([], [], lw=2, color='#33A02C', label='inhibitor')[0]]
ax.set_xlabel('$x$', fontsize=18)
ax.set_ylabel('concentration', fontsize=18)
ax.legend(framealpha=0.5)
# Initialize plotting
def init():
lines[0].set_data([],[])
lines[1].set_data([],[])
return lines
# This paints in the animation
def animate(i):
lines[0].set_data(x, a[i,:])
lines[1].set_data(x, h[i,:])
return lines
# call the animator.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=len(t), interval=20, blit=True)
# Save the animation
anim.save('gierer_meinhardt.mp4', fps=10, dpi=100,
extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
# Close the figure
plt.close()
# Encode/embed the animation
video = open("gierer_meinhardt.mp4", "rb").read()
video_encoded = b64encode(video)
video_tag = """
<video controls>
<source alt="G-M video" src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support MP4 H.264 video.
</video>
""".format(video_encoded)
HTML(data=video_tag)
We can also plot the result as kymographs. We will use the blue channel for activator and the green channel for inhibitor.
# RGB values for kymgraph image
rgb_a = np.zeros((len(t), len(x), 3))
rgb_h = np.zeros((len(t), len(x), 3))
# Values for the activator
rgb_a[:,:,2] = a / a.max()
# Values for the inhibitor
rgb_h[:,:,1] = h / h.max()
# Make kymograph
fig, ax = plt.subplots(1, 2, sharey=True)
ax[0].imshow(rgb_a, extent=(t[0], t[-1]/2, t[-1], t[0]))
ax[1].imshow(rgb_h, extent=(t[0], t[-1]/2, t[-1], t[0]))
ax[0].set_ylabel('t (s)')
ax[0].set_title('activator')
ax[1].set_title('inhibitor')
for i in [0, 1]:
ax[i].set_xlim((x[0], x[-1]))
ax[i].set_xticks([t[0], t[-1]/2])
ax[i].set_xticklabels([x[0], x[-1]])
ax[i].set_xlabel(u'x (µm)')