21. Scaling reaction-diffusion patterns

Design principle

  • Expansion repression allows for scaled profiles.


  • Scaling in morphogenesis.

import numpy as np
import scipy.integrate

import biocircuits

import bokeh.layouts
import bokeh.plotting
import bokeh.io

Loading BokehJS ...

Reaction-diffusion models for morphogen patterning

In the last lesson, we discussed morphogens, chemical species that determine cell fate in a concentration-dependent manner. Morphogens are spatially distributed throughout a developing organism, thereby imparting spatial variation in cell types. They present us with a great opportunity to connect our genetic network dynamics to spatial variation.

We specifically discussed two-component systems that give what is now known as Turing patterns. More generally, we can consider a set of biochemical species, \(\mathbf{M} = \{M_1, M_2, \ldots\}\). The set of concentrations of these biochemical species is \(\mathbf{c} = \{c_1, c_2, \ldots\}\). We can write down a system of partial differential equations describing the dynamics of morphogen concentrations (which is more or less a statement of conservation of mass) as

\begin{align} \frac{\partial c_i}{\partial t} &= \frac{\partial}{\partial x}\,D_i\,\frac{\partial c_i}{\partial x} + r_i(\mathbf{c}), \end{align}

valid for all \(i\). Here, \(r_i(\mathbf{c})\) describes the net rate of production of species \(M_i\), dependent on the concentrations of all other species. The diffusion coefficient \(D_i\) parametrizes the diffusion of species \(i\). For our purposes today, we are restricting ourselves to one-dimensional systems. We also have implicitly assumed that the diffusivity tensor is diagonal. Finally, for all of our analyses in this lecture, we will assume \(x \in [0,L]\) with constant flux (usually no flux) boundary conditions at the ends of the domain. Recall that the flux of morphogen is given by Fick’s first law,

\begin{align} j_i &= -D_i\,\frac{\partial c_i}{\partial x}. \end{align}

A gradient from a point source

We will begin by analyzing the simplest case of morphogen patterning. We have a single morphogen, \(M\), which has a source at \(x = 0\). We model the source as a constant flux, \(\eta\) at \(x = 0\).

Simple decay of the morphogen

We begin with the case where our morphogen only undergoes a single reaction, degradation which we would represent as \(M \to \emptyset\), and has a constant diffusion coefficient. Our PDE is

\begin{align} \frac{\partial c}{\partial t} &= D \,\frac{\partial^2 c}{\partial x^2} - \gamma c, \\[1em] \left.\frac{\partial c}{\partial x}\right|_{x=0} &= -\frac{\eta}{D}\\[1em] \left.\frac{\partial c}{\partial x}\right|_{x=L} &= 0. \end{align}

We can readily solve for the unique steady state in this case as

\begin{align} c(x) = \frac{\eta}{\sqrt{D\gamma}}\,\frac{\mathrm{e}^{-x/\lambda} + \mathrm{e}^{(x-2L)/\lambda}}{1 - \mathrm{e}^{-2L/\lambda}}, \end{align}

where \(\lambda \equiv \sqrt{D/\gamma}\). So, we have an exponentially decaying morphogen profile with decay length \(\lambda\). Typically, \(L \gg \lambda\), so we have

\begin{align} c(x) &\approx \frac{\eta}{\sqrt{D\gamma}}\,\mathrm{e}^{-x/\lambda}. \end{align}

Two obvious problems with this expression

Considering that morphogens determine cell fate in a concentration-dependent manner, there are two obvious problems with this expression. 1. The total level of the morphogen is set by the strength of the source, \(\eta\), assuming \(D\) and \(\gamma\) are constant. This would mean that cell fate would not be robust to the strength of the source. 2. The decay length of the concentration profile is independent of size. This means that a larger organism would have all of its morphogen highly concentrated near \(x/L = 0\).

Issue 2 is most clearly seen if we plot two profiles together.

# x-values for plot
L_1 = 100
L_2 = 200
x_1 = np.linspace(0, L_1, 200)
x_2 = np.linspace(0, L_2, 200)

# Value of lambda
lam = 25

# Make plots
p1 = bokeh.plotting.figure(
    height=250, width=350, x_axis_label="x/L", y_axis_label="conc"
p2 = bokeh.plotting.figure(
    height=250, width=350, x_axis_label="x/L", y_axis_type="log"
p1.line(x_1 / L_1, np.exp(-x_1 / lam), legend_label="L=100", line_width=2)
    x_2 / L_2,
    np.exp(-x_2 / lam),
p2.line(x_1 / L_1, np.exp(-x_1 / lam), line_width=2)
p2.line(x_2 / L_2, np.exp(-x_2 / lam), color="orange", line_width=2)

bokeh.io.show(bokeh.layouts.gridplot([[p1, p2]]))

We will consider a solution to these issues proposed by Ben-Zvi and Barkai (PNAS, 2010).

The concept of expansion repression

Consider another molecule, called an expander, that facilitates expansion of the morphogen. Specifically, the expander enhances diffusion of the morphogen and also prevents its degradation. However—and this is important because we need (you guessed it) feedback—the morphogen represses production of the expander. Our adjusted differential equations are now

\begin{align} \frac{\partial m}{\partial t} &= \frac{\partial}{\partial x}\,D_m(m, e)\,\frac{\partial m}{\partial x} - \gamma_m(m, e)\, m \\[1em] \frac{\partial e}{\partial t} &= D_e \,\frac{\partial^2 e}{\partial x^2} + \frac{\beta_e}{1 + (m/k)^n} - \gamma_e\, e. \end{align}

Here, \(m\) is the concentration of morphogen and \(e\) is the concentration of expander,both of which are functions of both time \(t\) and position \(x\). We will consider the case where \(D_m\) is constant, where the morphogen enhances its own degradation, and the expander represses degradation, shown schematically below.


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? If we think about the construction of the gradient starting from no morphogen being present, the mechanism becomes clear. Initially, we have no morphogen present, so the morphogen concentration grows throughout the entire tissue. As morphogen is introduced from the end of the tissue, the expander gets locally depleted. This leads to the morphogen spreading more toward the distal end of the tissue, where there is more expander, which inhibits the morphogen’s degradation. However, as the morphogen levels toward the distal end rise, the expander can no longer limit the degradation of the morphogen because its production is inhibited by the elevated morphogen levels. Eventually, the expander level at the distal end is set by the balance of expander repression and morphogen degradation. This sets the level at the distal end, thereby setting the morphogen profile across the entire tissue.

In the absence of an expander, there is nothing to pin the value of the morphogen at the distal end, which means that the morphogen cannot “sense” how long the tissue is.

Numerical solution to R-D equations for simple degradation

We will first solve the system with simple degradation before treating the expander. To solve R-D systems, we can use the rd_solve() function from the biocircuits package. The function is general for 1D RD systems with a diagonal mobility tensor, but allows for nonconstant diffusion coefficients. We therefore need to provide the solver with a function to compute the reaction terms and the diffusion coefficients. For the oft used scenario where the diffusion coefficients are constant, we can use biocircuits.constant_diff_coeffs() for the latter.

We additionally need to provide the values of the derivatives at the boundaries to enforce the Neumann conditions. In this case, the flux is zero at the \(x = L\) boundary, but given by \(-\left.\eta\middle/ \sqrt{D}\right.\) at the \(x = 0\) boundary to give the constant flux of morphogen.

def degradation_rxn(c_tuple, t, gamma):
    """Reaction expression for simple degradation."""
    return (-gamma * c_tuple[0],)

# Time points
t = np.linspace(0.0, 10.0, 100)

# Physical length of system
L = 5.0

# x for plotting
x = np.linspace(0, L, 500)

# Diffusion coefficient
diff_coeffs = (1.0,)

# Gamma
rxn_params = (1.0,)

# Initial concentration (500 grid points)
c_0_tuple = (np.zeros(500),)

# Production rate
eta = 1.0
derivs_0 = -eta / np.sqrt(diff_coeffs[0])

# Solve
conc = biocircuits.rd_solve(

And let’s make a plot. We will not plot the steady state because we already solved for that analytically and plotted it above, so we will show the dynamics with an interactive plot using the rd_plot() function of the biocircuits package. (The plot is not visible in the HTML rendering of this notebook.)

bokeh.io.show(biocircuits.xyt_plot(x, conc, t))