22. Scaling reaction-diffusion patterns

Design principle

  • Cooperative autodegradation of a morphogen retarded by a diffusing expander allows for scaled morphogen profiles.

Concepts

  • Scaling in morphogenesis.

[1]:

# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
    cmd = "pip install --upgrade biocircuits watermark"
    process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    stdout, stderr = process.communicate()
# ------------------------------

import numpy as np
import scipy.integrate

import biocircuits

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

bokeh.io.output_notebook()
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 versus the fraction length along the embryo, \(x/L\).

[2]:

# 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 plot
p = bokeh.plotting.figure(
    height=250, width=350, x_axis_label="x/L", y_axis_label="conc"
)
p.line(x_1 / L_1, np.exp(-x_1 / lam), legend_label="L=100", line_width=2)
p.line(
    x_2 / L_2,
    np.exp(-x_2 / lam),
    color="orange",
    legend_label="L=200",
    line_width=2,
)

bokeh.io.show(p)

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.

expander

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 autodegradation is minimal and the morphogen concentration is free to grow throughout the entire tissue. As morphogen is introduced from the proximal (left, in our plots) 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 encountered 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.

[3]:

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(
    c_0_tuple,
    t,
    L=L,
    derivs_0=derivs_0,
    derivs_L=0,
    diff_coeff_fun=biocircuits.constant_diff_coeffs,
    diff_coeff_params=(diff_coeffs,),
    rxn_fun=degradation_rxn,
    rxn_params=rxn_params,
)

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.

[4]:

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

Numerical solution of degradation with expander

To study the dynamics of the expander system, we will consider constant diffusion coefficients, and we take

\begin{align} \gamma_m(m, e) = (\gamma_{m,1} + \gamma_{m,2}m)\,\frac{m}{1+e}. \end{align}

This means that in addition to natural degradation of the morphogen, we also have enhanced autodegradation, both of which are retarded by the expander. Thus, our PDEs are

\begin{align} &\frac{\partial m}{\partial t} = D_m \,\frac{\partial^2 m}{\partial x^2} - (\gamma_{m,1} + \gamma_{m,2}m)\,\frac{m}{1+e}, \\[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}

Note the \(m^2\) term in the degradation of the morphogen. This serves to make the shape of the gradient more robust to the source strength. Momentarily, we will see that this robustness is lost when ultrasensitive degradation is eliminated.

We will perform the calculation for two values of \(L\) and see how the gradient scales.

[5]:

def expander_gradient_rxn(
    c_tuple, t, gamma_m1, gamma_m2, gamma_e, beta_e, k, n
):
    """
    Reaction expression for simple degradation.
    """
    # Unpack
    m, e = c_tuple

    # Morphogen reaction rate
    m_rate = -(gamma_m1 + gamma_m2 * m) * m / (1 + e)

    # Expander production rate
    e_rate = beta_e / (1 + (m / k) ** n) - gamma_e * e

    return (m_rate, e_rate)


# Time points
t = np.linspace(0.0, 5e5, 300)

# Physical length of system and grid points
L_1 = 100.0
L_2 = 200.0
n_gridpoints = 100

# Diffusion coefficients
diff_coeffs = (1, 10)

# Reaction params
gamma_m1 = 0
gamma_m2 = 100
gamma_e = 1e-6
beta_e = 1e-2
k = 1e-3
n = 4
rxn_params = (gamma_m1, gamma_m2, gamma_e, beta_e, k, n)

# Initial concentration
c_0_tuple = (np.zeros(n_gridpoints), np.zeros(n_gridpoints))

# Production rate
eta = 10.0

# Solve
c_tuple_L_1 = biocircuits.rd_solve(
    c_0_tuple,
    t,
    L=L_1,
    derivs_0=(-eta, 0),
    derivs_L=0,
    diff_coeff_fun=biocircuits.constant_diff_coeffs,
    diff_coeff_params=(diff_coeffs,),
    rxn_fun=expander_gradient_rxn,
    rxn_params=rxn_params,
    atol=1e-7,
)

c_tuple_L_2 = biocircuits.rd_solve(
    c_0_tuple,
    t,
    L=L_2,
    derivs_0=(-eta, 0),
    derivs_L=0,
    diff_coeff_fun=biocircuits.constant_diff_coeffs,
    diff_coeff_params=(diff_coeffs,),
    rxn_fun=expander_gradient_rxn,
    rxn_params=rxn_params,
)

With the results in hand, we can compute the morphogen levels and see if we have scaling.

[6]:

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

# Make plots
p1 = bokeh.plotting.figure(
    height=250, width=350, x_axis_label="x/L", y_axis_type="log"
)
p2 = bokeh.plotting.figure(
    height=250, width=350, x_axis_label="x", y_axis_type="log"
)
p1.line(x_1 / L_1, c_tuple_L_1[0][-1, :], legend_label="L=100", line_width=2)
p1.line(
    x_2 / L_2,
    c_tuple_L_2[0][-1, :],
    color="orange",
    legend_label="L=200",
    line_width=2,
)
p2.line(x_1, c_tuple_L_1[0][-1, :], line_width=2)
p2.line(x_2, c_tuple_L_2[0][-1, :], color="orange", line_width=2)

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

The slope of the gradient is similar for both values of \(L\) in scaled (\(x/L\)) units. Furthermore, the absolute levels remain close to each other. So, we have achieved (approximate) scaling.

However, the fall-off of the morphogen concentration from \(x = 0\) is very rapid, much more rapid that we would expect for an exponential profile. This can be seen more clearly if we plot the profiles on a linear scale.

[7]:

p1 = bokeh.plotting.figure(height=250, width=350, x_axis_label="x/L")
p2 = bokeh.plotting.figure(height=250, width=350, x_axis_label="x")
p1.line(x_1 / L_1, c_tuple_L_1[0][-1, :], legend_label="L=100", line_width=2)
p1.line(
    x_2 / L_2,
    c_tuple_L_2[0][-1, :],
    color="orange",
    legend_label="L=200",
    line_width=2,
)
p2.line(x_1, c_tuple_L_1[0][-1, :], line_width=2)
p2.line(x_2, c_tuple_L_2[0][-1, :], color="orange", line_width=2)

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

This also makes more clear the ≈30% difference in morphogen concentration at \(x = 0\) for the two different embryo lengths, though the concentrations become very close even for small nonzero \(x\).

Let’s make a movie to explore the dynamics.

[8]:

# Normalize concentrations
c_tuple_L_1_norm = tuple([c / c.max() for c in c_tuple_L_1])

bokeh.io.show(
    biocircuits.xyt_plot(
        x_1,
        c_tuple_L_1_norm,
        t,
        x_axis_label='x',
        y_axis_label='normalized concentration',
        legend_names=["morphogen", "expander"],
    )
)

As can be seen in the movie, the pattern is established rather quickly, with the expander moving toward its steady state (uniform) concentration.

Numerical solution for expander with non-cooperative degradation

We now consider the case where we do not have cooperative degradation. I.e., we have

\begin{align} \gamma_m(m, e) = \gamma_{m,1}\,\frac{m}{1+e}. \end{align}

Thus, our PDEs are

\begin{align} &\frac{\partial m}{\partial t} = D_m \,\frac{\partial^2 m}{\partial x^2} - \gamma_{m,1}\frac{m}{1+e} \\[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}

We can reuse the functions we already wrote with \(\gamma_{m,1} >0\) and \(\gamma_{m,2} = 0\).

[9]:

# Reaction params
gamma_m1 = 0.1
gamma_m2 = 0
gamma_e = 1e-6
beta_e = 1e-2
k = 1e-3
n = 4
rxn_params = (gamma_m1, gamma_m2, gamma_e, beta_e, k, n)

# Initial concentration
c_0_tuple = (np.zeros(n_gridpoints), np.zeros(n_gridpoints))

# Production rate
eta = 0.1

# Solve
c_tuple_L_1 = biocircuits.rd_solve(
    c_0_tuple,
    t,
    L=L_1,
    derivs_0=(-eta, 0),
    derivs_L=0,
    diff_coeff_fun=biocircuits.constant_diff_coeffs,
    diff_coeff_params=(diff_coeffs,),
    rxn_fun=expander_gradient_rxn,
    rxn_params=rxn_params,
    atol=1e-6,
)

c_tuple_L_2 = biocircuits.rd_solve(
    c_0_tuple,
    t,
    L=L_2,
    derivs_0=(-eta, 0),
    derivs_L=0,
    diff_coeff_fun=biocircuits.constant_diff_coeffs,
    diff_coeff_params=(diff_coeffs,),
    rxn_fun=expander_gradient_rxn,
    rxn_params=rxn_params,
    atol=1e-6,
)

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

# Make plots
p1 = bokeh.plotting.figure(
    height=250, width=350, x_axis_label="x/L", y_axis_type="log"
)
p2 = bokeh.plotting.figure(
    height=250, width=350, x_axis_label="x", y_axis_type="log"
)
p1.line(x_1 / L_1, c_tuple_L_1[0][-1, :], legend_label="L=100", line_width=2)
p1.line(
    x_2 / L_2,
    c_tuple_L_2[0][-1, :],
    color="orange",
    legend_label="L=200",
    line_width=2,
)
p2.line(x_1, c_tuple_L_1[0][-1, :], line_width=2)
p2.line(x_2, c_tuple_L_2[0][-1, :], color="orange", line_width=2)

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

Without ultrasensitive degradation, the profiles are exponential (as seen by the linear profiles in the semi-log plot). Their slopes are similar, but clearly different. The absolute level of morphogen also varies throughout the embryo. That is, the profiles when plotted against \(x/L\) do not rapidly come together. This is evident when plotted on a linear scale.

[10]:

p1 = bokeh.plotting.figure(height=250, width=350, x_axis_label="x/L")
p2 = bokeh.plotting.figure(height=250, width=350, x_axis_label="x")
p1.line(x_1 / L_1, c_tuple_L_1[0][-1, :], legend_label="L=100", line_width=2)
p1.line(
    x_2 / L_2,
    c_tuple_L_2[0][-1, :],
    color="orange",
    legend_label="L=200",
    line_width=2,
)
p2.line(x_1, c_tuple_L_1[0][-1, :], line_width=2)
p2.line(x_2, c_tuple_L_2[0][-1, :], color="orange", line_width=2)

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

Therefore, we have arrived at a design principle: Cooperative autodegradation of a morphogen that is retarded by a diffusing expander can give scaling of a morphogen gradient.

Computing environment

[11]:

%load_ext watermark
%watermark -v -p numpy,scipy,numba,bokeh,biocircuits,jupyterlab

Python implementation: CPython
Python version       : 3.8.10
IPython version      : 7.22.0

numpy      : 1.20.2
scipy      : 1.6.2
numba      : 0.53.1
bokeh      : 2.3.2
biocircuits: 0.1.5
jupyterlab : 3.0.14