4  Clock-and-wavefront patterning

4.1 Goal of the tutorial:

In this tutorial you will look at gradient formation and patterning using different mechanisms than you have seen previously. We will model the so-called clock-and-wavefront pattern, which stems from oscillations in gene expression and gene products and results in a regularly striped segmentation pattern in a growing tissue. Note that as compared to the Drosophila French Flag type case, this is thought to be the evolutionary ancestral model of segmenting animal body axes.

4.2 The model system

The clock-and-wavefront model is an important model in describing somitogenesis. In this process early in embryo development, the somites, a precursor tissue for the vertebrae and other tissues later in development, are formed from the pre-somitic mesoderm. This mesodorm extends on the posterior end by growth, and the somites bud off periodically at the anterior end in the order of a couple of weeks (in humans).

Key players in this model system are the protein FGF (fibroblast growth factor), and many genes and gene products that have an oscillatory pattern that will determine cell fate. However, for this tutorial, we simplify this to a single pair of gene mRNA and product, denoted with \(m\) and \(p\) for short.

4.3 Programming with classes

Because we are going to make a more complex model with a tissue existing of multiple cells, and each cell having its own concentrations of FGF, \(m\) and \(p\), we are going to use Classes in our code. You have been using Classes already: the data types such as int, str and bool have their own class, and the str class has many methods (=functions working on that class) defined, such as "hello world".upper(), but it is also possible to create custom classes. With classes, you can easily make objects, which is part of the object-oriented programming paradigm.

Today, we are going to use classes for the different levels of our model: 1) tissue, 2) cell, 3) \(m\) & \(p\) clock and 4) the plotting. By using classes, we can separate things that happen on a tissue/cellular/clock scale, and separate the model from the visualization of it. You will first work with 1 & 2 & 4, then 3 on its own and then combine all four yourself into one model.

Important concepts when working with classes are

  • Class versus instance
  • Defining the __init__ method and other methods
  • Class attributes and using the keyword self

This tutorial should be doable without an extensive knowledge of classes as there are plenty of examples to copy-paste from, but feel free to read up on these concepts here at this online tutorial

4.4 Questions

Exercise 4.1 (Biology - An alternative gradient forming mechanism) In practical 1 we saw how gradients can be created through local production, diffusion and decay. However, other mechanisms for gradient formation are possible, such as cell lineage transport. Here we work with a model for the FGF gradient (01-fgfgradientfromgrowth.py) where only the rightmost/posterior cell produces FGF and grows, and in which cells upon division inherit this FGF from their mother cell. First study the code to see how it uses Classes and see what happens. Next, play with the model by varying the decay rate of FGF and the division rate of the cells. How does this affect the gradient? What happens if divisions are only allowed during the first half of the simulation (put divisiontime to 0.5 instead of 1).

Exercise 4.2 (Mathematics) Let us now move to the other half of the clock and wavefront model, the clock part (02_clock.py), in which we implemented one of the earliest models for the somitogenesis oscillator from Lewis (2003) which models a gene that codes for a mRNA (\(m\)) that encodes a protein (\(p\)) that acts as a repressive transcription factor of this same gene. In class we discussed how for oscillations negative feedbacks, delays and non-linearity are important. Examine the code to find the differential equations governing this model and determine the negative feedback, delay and non-linearities in them.

Exercise 4.3 (Biology) Play with the parameters of the model. How does the delay (\(\tau\)) affect oscillations?

Exercise 4.4 (Algorithmic thinking) In the file 03_rolling_clock.py, there is a different implementation of the clock. Compare the two files and find out how they differ. What benefits for studying the model does 02_clock.py have over 03_rolling_clock.py and vice versa? Ignore the added functions __copy__ and set_tau in this comparison. Some differences become clearer when you run the code too.

Exercise 4.5 Next, we will combine the clock and FGF wavefront into a single model.

Read and interpret existing code: Read the code of 04_clockplusfgf.py and try to understand the assumptions from this model implementation by answering the following questions:

  1. Each cell has its own clock. What clock states is the tissue initialized with? And what clock states do newly divided cells get?

  2. The FGF wavefront affects \(\tau\) : what function is used for that? What are your expectations for the effect of the FGF wavefront on the cells’ clocks?

  3. Does growth affect the clock state?

  4. What is your opinion on the assumptions discussed in the three questions?

Exercise 4.6 (Biology) Describe how the model behaves and why. Reverse the effect of FGF and see how that affects the patterning? Do we get stable somites?

Exercise 4.7 (Biology & Programming) What is still missing is a means to transform the temporal oscillations in the posterior of the tissue into a spatial pattern in the anterior. In the French Flag morphogen gradient lecture we discussed that memory mechanisms are important to stabilize spatial patterns once the start up signal that broke the symmetry and initialized patterning has gone. It turns out that such a memorization/stabilization mechanism is also essential to convert oscillations to stripes.

  1. In the code 05_clockplusfgf.py (available in the Modeling Life Teams channel, and also here in the dropdown text below) we added an extra ‘memory’ molecule \(M\) as a new attribute to the Class Cell, which value is inherited from the mother cell upon division. Write your own code to perform its updating inside the function run_clocks. Remember that to create memory we need bistability, which can be easily achieved by having the memory molecule having a non-linear saturating positive feedback on itself. Study what happens (you can see what happens since we already incorporated visualisation of \(M\) into the python code). What happens for different initial values of \(M\)?
from rolling_clock import Clock # import the Clock class from rolling_clock.py
from graphics import * # import the visualization functions from graphics.py

# Similation parameters
totaltime=50*3600
divisiontime=1 # fraction of simulation time divisions are allowed
dt=1
time_steps = int(totaltime / dt)
time_visualization = 600 # How often to visualize the tissue growth in seconds

growth_rate =128/(10*3600)
cell_width=4
doubling_threshold = 2*cell_width

#when modifying the degradation rate, modify the production rate
#similarly, that way p/d stays constant and you have same max level
#at the right end of the tissue if no growth takes place
pfgf=60*0.001
dfgf=4*0.00001

# Parameters for the clock
alpha = 1.0/60
beta = 1.0/60
K = 1.0
n=4
tau=20*60
mu=0.03/60
v=0.03/60


class Cell:
    def __init__(self, fgf=0, m=0, p=0, memory=0):
        self.fgf = fgf
        self.xmin = 0
        self.xmax = 0
        self.xheight = 0
        self.clock = Clock(alpha, beta, K, n, tau, mu, v, m0=m, p0=p)  # Initialize the clock for each cell
        self.memory = memory  # Initialize memory value

    def grow_cell(self, dt, growth_rate):
        self.fgf = self.fgf * self.xheight / (self.xheight + dt * growth_rate)  # Dilute content with volume
        self.xmax += dt * growth_rate
        self.xheight = self.xmax - self.xmin

    def update_clock_tau(self):
        self.clock.set_tau((1+0.5*(100-self.fgf)/100)*tau) # Update tau based on FGF level and global tau

class Tissue:
    def __init__(self, doubling_threshold=doubling_threshold, time_visualization=time_visualization, growth_rate=growth_rate, dt=dt):
        self.cells = []
        self.num_cells = 0
        self.xmax_tissue = 0
        self.doubling_threshold = doubling_threshold
        self.time_visualization = time_visualization
        self.growth_rate = growth_rate
        # self.dt=dt

    # some functions to update the tissue properties
    def update_xmax_tissue(self):
        if self.cells:
            self.xmax_tissue = max(cell.xmax for cell in self.cells)
        
    # Functions to add, insert, and divide cells
    def add_cell(self, cell):
        self.cells.append(cell)
        self.num_cells += 1
        self.update_xmax_tissue()
    
    def insert_cell(self, cell, index):
        # Insert a cell at the specified index
        self.cells.insert(index, cell)
        self.num_cells += 1
        self.update_xmax_tissue()

    def divide_cell(self, cell_index):
        # Divide the cell at the specified index
        cell = self.cells[cell_index]
        mid_point = (cell.xmin + cell.xmax) / 2
        new_cell = Cell(fgf=cell.fgf)  # Create a new cell with the same FGF
        new_cell.clock=cell.clock.__copy__()  # Copy the clock state to the new cell

        new_cell.xmin = mid_point
        new_cell.xmax = cell.xmax
        new_cell.xheight = new_cell.xmax - new_cell.xmin
        # Adjust the original cell
        cell.xmax = mid_point
        cell.xheight = cell.xmax - cell.xmin
        # Add the new cell to the list of cells
        self.insert_cell(new_cell, cell_index + 1)  # Insert after the original cell
        return new_cell


    # functions to initialize and grow the tissue 
    def initialize_regular_tissue(self, num_cells=1, initial_fgf=100, cell_width=cell_width):
        # Initialize a regular tissue with a specified number of cells of which the first has initial FGF
        # and the rest have FGF=0
        for i in range(num_cells):
            if i <1:
                cell = Cell(fgf=initial_fgf)  # First cell has initial FGF
            else:
                cell = Cell(fgf=0)  # Rest have FGF=0
            cell.xmin = i * cell_width
            cell.xmax = (i + 1) * cell_width
            cell.xheight = cell_width
            self.add_cell(cell)
        self.update_xmax_tissue()

    def grow_tissue(self, dt, growth_rate, doubling_threshold,n=1):
        for cell_index,cell in enumerate(self.cells[-n:],-n):  # Apply growth to the n last cells only
            # Grow the cell
            cell.grow_cell(dt, growth_rate)
            # Update xmin and xmax of all other cells above it "pushing" them
            if cell_index != - 1:  # Avoid index out of range
                for other_cell in self.cells[cell_index+1:]:
                    if other_cell.xmin > cell.xmin and other_cell.xmin < cell.xmax:
                        other_cell.xmin += dt*growth_rate
                        other_cell.xmax += dt*growth_rate
                    elif other_cell.xmin >= cell.xmax:
                        other_cell.xmin += dt*growth_rate
                        other_cell.xmax += dt*growth_rate
            self.update_xmax_tissue()


        # After growing the tissue, check for cell division
        for cell_index, cell in enumerate(self.cells):
            # Check if the cell has reached the doubling threshold
            if cell.xheight >= doubling_threshold:
                # Print which cell is being divided
                # print(f"Dividing cell at index: {cell_index}")
                # Create a new cell for division
                self.divide_cell(cell_index)
                # Add the new cell to the list of new cells

            # print(f"Current number of cells: {self.num_cells}")


    # Simulate the morphogen gradient over time
    # The last cell produces fgf, the others only degrade it
    def simulate_morphogen(self, time_steps, pfgf, dfgf):
        for i, cell in enumerate(self.cells):
            if i == len(self.cells) - 1:  # only last cell produces fgf
                dtfgf = dt * (pfgf  - (dfgf * cell.fgf))
            else: # other cells degrade fgf
                dtfgf = dt * (- (dfgf * cell.fgf))
            cell.fgf += dtfgf

    def run_clocks(self, t, dt):
        # Update the clocks of all cells
        for cell in self.cells:
            cell.update_clock_tau()
            cell.clock.simulate(t, dt)
            #write a method for updating the memory 
            #cell.memory +=

# Modify simulate_development to include plotting
def simulate_development(time_steps):
    # Initialize a tissue
    tissue = Tissue()
    # Initialize the tissue with a single cell
    tissue.initialize_regular_tissue(num_cells=1, initial_fgf=100, cell_width=cell_width)

    plt.ion()
    tissueplot=TissuePlot(tissue, time_visualization, dt,  n_axis=3)
    tissueplot.initialize_axis_cell_data(axis_index=0, colormap='viridis', attribute='fgf',update_direction={"max_up":1,"max_down":0.75 }, label='FGF Level')
    tissueplot.initialize_axis_cell_data(axis_index=1, colormap='Greys', attribute='clock.p_values[-1]',update_direction={"max_up":1.3,"max_down":0.6} , label='P Level')
    tissueplot.initialize_axis_cell_data(axis_index=2, colormap='Greys', attribute='memory' , update_direction={"max_up":1,"max_down":0}, label='Memory')
    for t in range(time_steps):
        if t / time_steps < divisiontime:  # Allow divisions only in the first fraction of the simulation time
            tissue.grow_tissue(dt, growth_rate, doubling_threshold)
        tissue.simulate_morphogen(time_steps, pfgf, dfgf)   
        tissue.run_clocks(t * dt, dt)     
        if t *dt % time_visualization == 0:  # Plot every something timesteps
            tissueplot.update_xmax_display() # Update the maximum x-axis limit for display purposes
            tissueplot.update_plot_cell_data(t,axis_index=0)
            tissueplot.update_plot_cell_data(t,axis_index=1)
            tissueplot.update_plot_cell_data(t,axis_index=2)
    print("Simulation finished")
    # Finalize and show the plot
    plt.ioff()
    # plt.show()

# Call the function to simulate tissue growth and plot
simulate_development(time_steps)
  1. Thus, we need to make \(M\) dependent on the clock. To achieve this, we can start with low values of \(M\) that prevent autoactivation, and then have \(M\) activated by either \(p\) (part of the clock) or \(M\) itself using the following function: \[\frac{dM}{dt}=c\max\left(\frac{p^4}{h_p^4+p^4}, \frac{M^4}{h_M^4+M^4}\right)-\delta M.\]

What happens now?

  1. Since oscillations occur in the PSM and stripe formation occurs only more anteriorly we should constrain memorization to occur only below a certain \(FGF\) value. Implement this in the code. Did this help?

  2. So in addition to have memorization occur only below a certain \(FGF\) value it should also occur above another, lower \(FGF\) value, constraining it to occur in a limited temporal window that enables cells passing through there with a high \(p\) state to induce a high \(M\) state while cells passing through with a low \(p\) state to not induce a high \(M\) state. Implement this and see if you get stable somite formation. What would be another way of achieving this?

Exercise 4.8 (Biology) Both zebrafish and mice are common model organisms, so we know a lot about the biological parameters of their somitogenesis. See the following list:

Parameter Zebrafish Reference Mouse Reference
Duration of somitogenesis 18 h (a) 5 days (b)
Number of somites ~30 (a) 65 (c)
Somite size 50 \(\mu\) / 30 \(\mu\) (d) (a) 120 \(\mu\) (c)
Cells per somite ~5 cell in diameter (a) 5-11 (estimated from total cell size in 3D, ranging 1 order of magnitude) (c)
Clock period 25 min /30 min (d) / (e) / (f) 2-3 h (e)

From these parameters, you can derive a number of desired model inputs/outcomes:

  • The total size of the tissue at the end of somitogenesis
  • The size of a cell
  • Speed of division

Use the model and try to find suitable model parameters to recreate the development of both zebrafish and mice: is this model able to describe both of these processes? I.e., is this model able to deal with the scale differences between zebrafish and mice?

A couple of notes:

  • You don’t have to exactly recreate the biological parameter with 100% precision, except for the number of segments/somites, although it can be a fun challenge to get a complete match.
  • You might want to adapt your plotting timestep to have sufficient but not too many plot updates in one simulation.
  • You can test the clock parameters separately with the clock.py/rolling_clock.py script. Don’t forget that FGF has an effect on \(\tau\) !

Exercise 4.9 (Questions for master students) The model currently has a number of assumptions that we can question. Feel free to explore any of these further open questions and study how it effects the outcome of the model:

  • What if the relationship between FGF and \(\tau\) is shaped differently? For instance, if the clock runs faster on the left/anterior than on the right/posterior?
  • What if the clock state of a daughter cell is started fresh rather than being a copy from the mother cell?
  • What if the clock \(m\) and \(p\) are diluted by growth?
  • What if \(\tau\) is unaffected by FGF: can we still get somites fixed in place?

4.5 Relevant literature

If you want to know more about the model system and previous models, have a look at the following (after the tutorial):

Lewis (2003) Autoinhibition with transcriptional delay: a simple mechanism for the zebrafish somitogenesis oscillator.

Hester (2012) A Multi-cell, Multi-scale Model of Vertebrate Segmentation and Somite Formation.

Herrgen et al. (2010) Intercellular coupling regulates the period of the segmentation clock.

Soroldoni et al. (2014) Genetic oscillations. A Doppler effect in embryonic pattern formation.

Sonnen et al. (2018) Modulation of Phase Shift between Wnt and Notch Signaling Oscillations Controls Mesoderm Segmentation.

Bulusu et al. (2017) Spatiotemporal Analysis of a Glycolytic Activity Gradient Linked to Mouse Embryo Mesoderm Development.

Oostrom et al. (2025) Coupling of cell proliferation to the segmentation clock ensures robust somite scaling.

4.6 References for biological parameters:

  1. Stickney, Barresi, and Devoto (2000)

  2. Saga (2012)

  3. Tam (1981)

  4. Ishimatsu et al. (2018)

  5. Carraco, Martins-Jesus, and Andrade (2022)

  6. Tomka, Iber, and Boareto (2018)