2 Gradients and segments

2.1 Morphogen Gradients and Patterning

In this practical, we a going to look at how an organism can form segments along its body axis . The mathematical model that we will implement and study is an implementation of the so-called French flag conceptual model first proposed by Lewis Wolpert (Wolpert 1969). It assumes the spatially graded expression of a morphogen “M” that influences the expression of some downstream genes A, B and C. Their expression is often visualized by red, white and blue and the arising pattern resembles the French flag, hence the name (see the power of visualization).

One of the most well-studied organisms when it comes to body axis segmentation (although its segmentation mechanism is evolutionary derived and a-typical!) is the development of the fruit fly Drosophila melanogaster. Supporting the ideas of Wolpert, it was found that through tethering maternal Bicoid mRNA to one side of the embryo, Bicoid protein can form a gradient extending along the anterior-posterior axis, with so called gap genes as a first tier in the segmentation hierarchy differentially responding to different Bicoid protein levels (Driever and Nüsslein-Volhard 1988). Later it was found that often at least two opposing morphogen gradients drive downstream gene expression and genes typically respond to multiple inputs.

2.2 Mathematical modeling - integrating multiple signals

Promotors/enhancers driving gene expression frequently make use of so called OR and AND gates to integrate inputs from different transcription factors. An OR gate can be implemented mathematically with a sum of the effect of the transcription factors, while an AND can be implemented mathematically with a product. Some examples:

\(\frac{dX}{dt} = a(\text{tf1}) + b(\text{tf2})\): Gene X is induced by transcription factor 1 and 2 in an OR fashion, either \(a(\text{tf1})\) or \(b(\text{tf2})\) needs to be high to give high transcription of X.
\(\frac{dY}{dt} = a(\text{tf1})\cdot b(\text{tf2})\): Gene Y is induced by transcription factor 1 and 2 in an AND fashion, both \(a(\text{tf1})\) and \(b(\text{tf2})\), which are being multiplied, need to be high to give high transcription of X.

Note that the shape of \(a(\text{tf1})\) and \(b(\text{tf2})\) (increasing or decreasing function of the transcription factor) determines whether tf1 and tf2 are repressing or activating.

2.3 Python code

The code below is also in 01_morphogengradient_to_segments.py.

Starting code for this practical

import numpy as np
import matplotlib.pyplot as plt

# Parameters
Lx = 40.0  # Length of the domain in x in microm
Ly = 10.0  # Length of the domain in y in microm
T = 200  # Total time in seconds
dx = 0.5  # Grid spacing in x
dt = 0.1  # Time step
nx = int(Lx/dx)+2  # Number of grid points in x + padding grid points
ny = int(Ly/dx)+2  # Number of grid points in y + padding grid points
# Padding grid points to account for boundary conditions
nt = int(T/dt)  # Number of time steps
D = 0.4  # Diffusion coefficient in mm^2/s
decayM =0.01 # Decay rate in 1/s


# Parameters for A, B, C
... # TODO create parameters for A, B, C as needed in Q5

# Stability criterion
if D * dt / dx**2 > 0.5:
    raise ValueError("Stability criterion not met")

# A, B and C are required for later exercises.
A = np.zeros((nx, ny))
B = np.zeros((nx, ny))
C = np.zeros((nx, ny))

# Initial condition
u = np.zeros((nx, ny))
u[0, :] = 100

# Reaction-diffusion equation
def reaction_diffusion_step(u, D, dt, dx, decay):
    un = u.copy()
    u[1:-1, 1:-1] = un[1:-1, 1:-1] +  D *dt / dx**2 * (un[2:, 1:-1] + un[:-2, 1:-1] + \
                    un[1:-1, 2:]  + un[1:-1, :-2] - 4 * un[1:-1, 1:-1]) - \
                    decay * un[1:-1, 1:-1] * dt
    ## for loop version to understand the equation
    # for i in range(1, nx-1):
    #     for j in range(1, ny-1):
    #         u[i, j] = (un[i, j] +
    #                    D * dt / dx**2 * (un[i+1, j] + un[i-1, j] - 2 * un[i, j] +
    #                    un[i, j+1] + un[i, j-1] - 4 * un[i, j]) - decay * un[i, j] * dt)
    #boundary conditions
    u[-1, :] = (u[-2, :]/u[-3, :])*u[-2, :]  if sum(u[-3, :]) != 0 else np.zeros(ny)
    #to understand this line:
    #if sum(u[-3, :]) != 0:
    #    u[-1, :] = (u[-2, :]/u[-3, :])*u[-2, :]#extrapolate from third to last row
    #else:
    #    u[-1, :] = np.zeros(ny) #if already zero in third to last row, also zero in last row
    u[:, 0] = u[:, 1]
    u[:, -1] = u[:, -2]

    return u

def reaction_diffusion_gradient(t, u, D, dx, decay, switch_time = None, noise = False):
    '''
    Function to create a gradient in the u array that could decay after a certain time.
    t: current time step
    u: array to create the gradient in
    D: diffusion coefficient
    dx: grid spacing
    decay: decay rate
    switch_time: time step after which the gradient decays. If no switch is desired, set to None
    noise: whether to add noise to the gradient
    '''
    # TODO for student: write code for the noise and the switch.
    added_noise = np.zeros_like(u)  # Initialize noise array
    if noise:
        ...  # TODO: add noise generation code here for Q10
    
    if switch_time is None or t <= switch_time:
        # define a exponential decay gradient over the array in the x direction with numpy array operations using the index
        for i in range(u.shape[0]):
            u[i, :] = np.maximum(100 * np.exp(-i*dx/np.sqrt(D/decay))+added_noise[i, :], 0)
        return u
    if t > switch_time:
        ...# TODO Q7: implement a gradient that decays over time, otherwise return the original u array
        return u
    # In all other cases, return the original u array        
    return u

def hill(x, Km, pow):
    """Hill function for the reaction kinetics."""
    return (x**pow) / (Km**pow + x**pow) 

def ihill(y, Km, pow):
    """Inverse Hill function for the reaction kinetics."""
    return( (Km**pow) / (y **pow  + Km**pow))

# TODO for student: write update functions for A, B, C as needed in Q5


# initilize figure and axes for plotting
# TODO for student: Add a new axis for the ABC flag visualization as suggested in Q5
fig, (ax_M, ax_lines) = plt.subplots(2, figsize=(10, 8), gridspec_kw={'height_ratios': [3, 1]})  # Make the first graph 3 times the height of the second

# Time-stepping simulation loop
for n in range(nt):
    # Update all variables
    u = reaction_diffusion_step(u, D, dt, dx, decayM)
    # TODO for student: use precomputed gradient, update A, B, C as needed in Q5
    
    if n == 0:  # Initial plot
        imshow_M_plot = ax_M.imshow(u.T, cmap='viridis', origin='lower', aspect='auto')
        ax_M.set_title(f"Time: {n*dt}")
        ax_M.set_xlabel('x direction')
        ax_M.set_ylabel('y direction')
        ax_M.set_xticks([])
        ax_M.set_yticks([])

        # Plot the concentration at a specific y index (e.g., y=2)    
        line_plot = ax_lines.plot([x*dx for x in range(nx)], u[:, 2], label='M', color='green')
        # TODO: Add lines for A, B, C as needed in Q5

        
        ax_lines.legend(loc='upper right')
        ax_lines.set_ylim(0, 100)
        ax_lines.tick_params(axis='y')
        ax_lines.set_xlim(0, dx*nx)
        ax_lines.set_xlabel('x')
        ax_lines.set_ylabel('Concentration at y=2')
        ax_lines.tick_params(axis='x')

    if n % 20 == 0:  # Update plot every so many time steps
        #update the imshow M plot with the new data
        imshow_M_plot.set_data(u.T)
        ax_M.set_title(f"Time: {n*dt}")
            
        # Update the line plots with new data
        line_plot[0].set_ydata(u[:, 2])
        # TODO: Update A, B, C line plots as needed in Q5

    plt.pause(0.001)  # Pause to refresh the plot

# And keep the last plot open
# plt.show()

# Or close the plot window when done
plt.close(fig)

2.4 Questions

Exercise 2.1 (Algorithmic thinking) Have a look at the reaction-diffusion equation for the morphogen and the implementation of it in the file 01_morphogengradient_to_segments.py in the reaction_diffusion_step function, as well as at the initialization of the array \(u\). How is this different from the gradient formation modeling we discussed during the lecture? What is done at the terminal boundary in the length direction and why does this make sense? (hint: outcommented we provided code doing essentially the same but not using numpy arrays and hence written in a less compact matter to help you understand what is happening)

Exercise 2.2 (Important concept) Play with the morphogen \(M\) diffusion rate and the morphogen \(M\) decay rate and describe what happens. What happens in terms of dynamics and steady state if you change both, but the ratio stays the same? Hint: it might help to draw a horizontal line at a certain height to ease comparison.

Exercise 2.3 (Biology & mathematical thinking) Next, we are going to introduce the genes \(A\), \(B\) and \(C\) in the model. We want these genes to be expressed dependent on \(M\), and in the head, trunk and tail respectively, so \(A\) on the left, \(B\) in the middle and \(C\) on the right. Think of what conditions in terms of \(M\) should lead to expression of \(A/B/C\). What Hill functions (normal/inverse/any combination) corresponds to those conditions? Write down (pen and paper, not in code) full equations for the genes, do we need any other terms than just Hill functions, would we need specific parameter conditions?

Exercise 2.4 (Biology & algorithmic/mathematical thinking) From hereon we assume that the morphogen \(M\) gradient reaches steady state very quickly, and no longer use the numerical implementation of the diffusion equation and instead work with a superimposed morphogen profile defined by reaction_diffusion_gradient to save time. (Hint: to not call the function any more use # in front of where it is called). Change the simulation loop such that it computes the morphogen \(M\) gradient once. What type of function is the superimposed morphogen profile and how does this relate to question 2.

Exercise 2.5 (Biology & algorithmic/mathematical thinking)

Now create functions to update \(A\), \(B\) and \(C\) according to your equations from the previous question. You may use the predefined hill and ihill (inverse) functions that are provided in the code. For simplicity, you may keep most of the parameters the same across genes, but some have to be different to ensure the right location of the genes (see your reasoning to the previous question).
Also make sure that the levels of \(A\), \(B\) and \(C\) are updated in the simulation loop. (Hint: using array properties to update \(A\),\(B\) and \(C\), like in the reaction_diffusion_step function, makes your code run a lot faster than using for-loops)
Next, ensure \(A\), \(B\) and \(C\) are visualized in the bottom plot axis (copy and adapt the code for the visualization of \(M\)). You can add an extra third axis to the plot to visualize the (French) flag pattern, by getting which of the three genes is maximal at each location with np.argmax(np.array(\[A, B, C\]), axis=0) and turning that into an array of RGB colors of choice.
Do you get the expected “French flag” pattern? If not, think of why not and improve your equations from previous question, parameters or your code.

Exercise 2.6 (Biology) At some point in development, the morphogen gradient \(M\) will disappear. For example, in the case of Drosophila embryos the Bicoid gradient disappear because the maternal mRNA is degraded. Predict what will happen to the expression of \(A\), \(B\) and \(C\) (and hence the French flag pattern) if the morphogen \(M\) gradient disappears over time and assume \(A\), \(B\), and \(C\) are regulated by our equations (first try think about this without actually simulating this).

Exercise 2.7 (Algorithmic/mathematical thinking) Now write the code in reaction_diffusion_gradient that updates the morphogen \(M\) concentration, such that after the time point A, B and C have gotten close to equilibrium, the gradient gradually disappears. Adapt the simulation loop where necessary to implement this, and run your code:

Was your prediction on \(A\), \(B\), \(C\) from the previous question correct? Why/why not? Hint: you might need to increase the duration of your simulation, especially if it takes a long time for \(A\), \(B\) and \(C\) to reach equilibrium (or you can change the parameters to speed things up by using same production/degradation rate ratio yet higher absolute values of the individual parameters).

Exercise 2.8 (Biological & mathematical thinking) \(A\), \(B\) and \(C\) cannot remain in a stable pattern if they are only influenced by \(M\). How can we stabilize the pattern in absence of M? Test your ideas by creating new update functions for \(A\), \(B\) and \(C\) and let these new ‘rules’ kick in at the same time when \(M\) starts to decline. Again, you may find Hill functions useful and perhaps also the before/after switch time structure used in reaction_diffusion_gradient for \(M\). (Hint: think about how the genes should affect each other, and assume that in this new phase, when genes have already been initialized in absence of repression the genes will be expressed). Can you implement changes that maintain the expression domains of \(A\), \(B\) and \(C\) and does their shape change?

Exercise 2.9 (Biology & algorithmic/mathematical thinking) A sudden switch from phase 1 (stable \(M\) gradient) to phase 2 (decaying \(M\) gradient) resulting in the genes following different differential equations is biologically implausible. In reality, genes have complex promoters and enhancers integrating different inputs that arise in different developmental stages. Try to come up with one integrated expression for each gene, incorporating simultaneous input from the morphogen \(M\) gradient and the other genes. Create new update functions for \(A\),\(B\),\(C\) and test your ideas. Can you get a stable pattern before and after the decline in \(M\)? A couple of things you could consider:

Think of how positively and negatively \(M\) regulated expression behave once the gradient starts declining: what is the best way to combine (\(AND/OR\)) that with the regulation by the genes?
Consider splitting up the \(M\) regulated expression of middle gene \(B\) into a positive and negatively regulated morphogen part before integrating the other genes inputs.
It is also possible to give certain genes a bit of a constant boost to prevent their takeover by other genes due to timing issues

Exercise 2.10 (Biology & algorithmic/mathematical thinking) Write code to get noise in the morphogen \(M\) gradient, both in its steady state and decaying phase. During which phase do the expression domains of \(A\), \(B\) and \(C\) suffer more from the noise. Explain why? How could we make the system more robust in a manner that is also likely occurring in nature?

2.5 Extra questions

If you’re done early or a master student, you can make these extra questions. These questions need not be made in the order they are provided, you can choose what you would like to investigate.

Exercise 2.11 (Biological thinking) Play with the size of the domain to see what happens to the gene expression domains. What would this mean for an organism?

Exercise 2.12 (Biological & algorithmic thinking) Various mechanisms have been proposed to ensure that the domains of the morphogen controlled genes scale with the size of the domain. One proposed mechanism suggests the existence of an also diffusible “expander” molecule which expression is repressed by the morphogen but which itself either represses the degradation of the morphogen or enhances its diffusion (see e.g. https://www.pnas.org/doi/full/10.1073/pnas.0912734107 and https://onlinelibrary.wiley.com/doi/10.1111/dgd.12337).

For the expander we can write:

\[ \frac{∂E}{∂t}=p\frac{K^h}{K^h+M^h}-dE+D_E ∆E \]

Assume \(E\) reduces degradation of \(M\) (easier to implement than enhancement of diffusion) and study the effect of scaling. (Hint, vary the size of the domain in the length direction but plot domains as a function of relative instead of absolute domain size to compare domain sizes).

Exercise 2.13 (Algorithmic thinking) In b, we can also implement other boundary conditions, especially for the right boundary. What happens to the profile if we make a no-flux boundary by copying the value at \(n-2\) to \(n-1\)? And what if we set it to a sink (force concentration to zero)?

Exercise 2.14 (Algorithmic thinking) From question 4 onwards, how do things change if we do not assume a quasi-steady-state for the morphogen gradient (i.e. keep the morphogen dynamics instead of replacing it by the superimposed exponential)? How would you implement a disappearing gradient and how does it shape change the outcomes of the flag?

Exercise 2.15 (Algorithmic thinking) Implement your solution to make the system more robust for noise from question 10. Did it work?

2.6 Relevant literature / further reading

Dagmar Iber group, scaling from non-steady state dynamics and uniform growth: Fried and Iber (2014) https://www.nature.com/articles/ncomms6077

Cheung et al. (2011) https://pmc.ncbi.nlm.nih.gov/articles/PMC3109599/

He et al. (2015) https://www.nature.com/articles/ncomms7679

Bicoid gradient: larger eggs get more Bicoid mRNA so higher production rate, when it scales with volume this helps scale the gradient: Inomata (2017) https://onlinelibrary.wiley.com/doi/10.1111/dgd.12337

Expansion repression model. Players: Chordin, Bmp, Sizzled. Chordin represses Bmp which induces Sizzled (which has low decay rate), yet Sizzled reduces Chordin decay. Chordin is morphogen. Sizzled is expander (by reducing decay of morphogen) Ben-Zvi and Barkai (2010) https://www.pnas.org/doi/abs/10.1073/pnas.0912734107