colorlinks Department of Biochemistry and Molecular Medicine
Université de Montréal
MILA Québec

Department of Physics
McGill University

Montréal, Québec

<EMAIL_ADDRESS>
# Waves, patterns and bifurcations:
a tutorial review on the vertebrate
segmentation clock

Paul François and Victoria Mochulska

###### Contents

  1. 0 Verterbrate segmentation for theorists: why?
  2. 1 Characterizing vertebrate segmentation : clock, waves, morphogens
    1. 1 Early concepts
    2. 2 Statics and dynamics of metazoan segmentation
    3. 3 The segmentation clock paradigm
  3. 2 Early models
    1. 1 The clock and wavefront framework
    2. 2 Meinhardt’s model
    3. 3 Cell Cycle model
  4. 3 Phase models
    1. 1 From chemical equations to phase
    2. 2 Clock and Unclock
    3. 3 The Lewis Phase Model (LPM)
    4. 4 Flow-based phase models
    5. 5 Delayed coupled models
    6. 6 Doppler period shift
    7. 7 Wavefront as a phase shock emerging from coupling
    8. 8 Phase-amplitude coupling and excitability for oscillation arrest
  5. 4 From systems biology to geometric models
    1. 1 Simple delayed oscillators
    2. 2 Molecular (clock and) Wavefront models
    3. 3 Somite AP patterning: Inverse problem approach
    4. 4 Landscape geometry of segment formation
    5. 5 Clock and switch modulated
    6. 6 From geometric back to phase models ?
  6. 5 Hacking the segmentation clock
    1. 1 Monolayer PSM cell cultures : the $\alpha$ two-oscillator model
    2. 2 Entraining the segmentation clock
    3. 3 Inferring coupling rules : walkie talkie model
    4. 4 Exploring cell communications/coupling with optogenetics
    5. 5 In vitro mouse segmentation clock
    6. 6 Zebrafish cultures
    7. 7 Stem-cell systems
    8. 8 Randomization
  7. 6 Theoretical challenges and future insights
    1. 1 Categorizing models: primary waves and bifurcations
    2. 2 Reconciling models with experimental observations
    3. 3 Towards a universal model of the cellular dynamics during segmentation ?
    4. 4 Unknown and Known limitations, blind spots
    5. 5 Theory, learning and Evolution
  8. 7 Supplementary Materials
  9. 8 Acknowledgements
  10. A Some dynamical systems theory for biological oscillators
    1. Bifurcations and excitable systems
    2. Perturbing the phase
    3. Two demonstrations of Malkin theorem
  11. B Complementary Discussions
    1. Scaling Laws
    2. Number of waves in a growth model
    3. Doppler period shift calculations
    4. Conditions for oscillations for negative feedback oscillators with delay
    5. Delayed model with noise
    6. Practical details on the ERICA model

###### List of Figures

  1. 1 Sketches of Malpighi illustrating the process of somite formation in chick embryos
  2. 2 Definition of Kymographs
  3. 1 The French Flag Model
  4. 2 Fly segmentation
  5. 3 Vertebrate segmentation
  6. 4 Sketches of Molecular Players in Somitogenesis
  7. 5 Different wave patterns in different species
  8. 1 Qualitative view of the Clock and wavefront model.
  9. 2 Two-well landscape for the clock and wavefront model
  10. 3 Mathematical formulation of the Clock and Wavefront model.
  11. 4 Simulations of the Clock and Wavefront model.
  12. 5 Landscape Dynamics in the Clock and Wavefront model
  13. 6 Primary and Secondary waves
  14. 7 Meinhardt model and its reduction
  15. 8 Simulation of the reduced Meinhardt model
  16. 9 Flow-plots of the Meinhardt-Van der Pol model
  17. 10 Nucleation of a new border in the Meinhardt model
  18. 11 Scaling of the Meinhardt-Van Der Pol model
  19. 12 Cell Cycle model
  20. 1 Clock and Unclock
  21. 2 Lewis Phase Model
  22. 3 Contributions to the Doppler Shift of Segmentation period
  23. 4 Murray model
  24. 5 Phase Amplitude model
  25. 1 Delayed Oscillator
  26. 2 Evolved clock and switch model
  27. 3 PORD model
  28. 4 Building a Geometric model
  29. 5 Phase model with SNIC
  30. 1 mPSM culture
  31. 2 $\alpha$ model
  32. 3 Entraining the segmentation clock.
  33. 4 FitzHugh-Nagumo oscillator
  34. 5 Randomization Assay to infer coupling rules
  35. 1 Four possible scenarios for primary/secondary waves in somitogenesis
  36. 2 A possible ’Switch/Clock/Switch’ scenario
  37. 3 Geometry of a Mesp-2 like secondary wave
  38. 1 Type I excitability
  39. 2 Type II excitability

##  Chapter 0 Verterbrate segmentation for theorists: why?

The French naturalist Geoffroy Saint Hillaire noticed in the XIXth century a
universal feature of the body of many common animals [1]: they are primarily
built on the repeat of metaxmeric units along their anteroposterior axis.
Canonical examples include segments in arthropods, or our vertebrae. This
organization is so fundamental that entire phylogenetic groups have been named
in reference to units of their body plan, e.g. annelids or vertebrates. Fossil
records suggest that this segmental organization is an extreme form of organ
metamerism, that possibly accompanied the Cambrian explosion 600 million years
ago [2]. As such, metamerism can be considered a major evolutionary innovation
leading to modern animal life. The segmental organization is generally assumed
to provide multiple evolutionary advantages, for instance having multiple
connected body parts allows for versatile body movements, and division between
units allows for subsequent evolutionary specializations of individual
segments [3].

Vertebrae precursors in embryos are called somites, and the process of somite
formation is called "somitogenesis". Somites first appear as pairs of
epithelial spheres on both left and right sides of the neural tube, and
sequentially form from anterior to posterior during axis elongation, Fig. 1.
Multiple tissues derive from somites so a proper understanding and control of
somite formation might potentially lead to both fundamental advances and
practical application in regenerative medicine [4]. Somitogenesis is
particularly appealing to physicists for multiple reasons. As we will describe
below, it is now established that somitogenesis is tied to the presence of a
global genetic oscillator, called the segmentation clock [5], which is
associated with multiple waves propagating in embryonic tissues [6, 7, 8, 9].
The periodicity of this process further allows for multiple observations
within one single experiment, making it an ideal system for developmental
biophysics. Examples of experimental perturbations include recovery of
oscillation following perturbations [10, 11] and entrainment [12]. Individual
cells can oscillate when dissociated [13], and it is now clear that the
segmentation clock at the tissue level is an emergent, self-organized process
[14]. Somites are in fine well-defined physical units, so somitogenesis also
presents a nice example of interaction between genetic expression, signaling,
and biomechanical process leading to morphogenesis. Lastly, it should be
pointed out that the existence of an oscillator controlling somite formation
has been predicted theoretically using advanced mathematical concepts
(catastrophe theory) [15] 21 years before its definitive experimental proof
[5]. So vertebrate segmentation is a good example of "Figure 1" scientific
endeavor [16], where theoretical predictions suggest experiments, and where a
fruitful back and forth between experimental biology and theoretical modeling
has occurred.

Important recent advances include more controlled experimental setups such as
explants [17, 18], stem cell systems [19, 20] and even synthetic developmental
biology assays [21] such as somitoids/segmentoids [22, 23]. Feynman’s famous
quote "What I cannot create, I do not understand" is often invoked (see e.g.
[24, 25]) to motivate such in vitro reconstruction of biological systems.
Indeed, great insights can be drawn by creating and manipulating minimal
experimental models. It is however important to stress that this quote, found
on Feynman’s blackboard upon his death [26], likely reflects the mindset of a
theoretical physicist, further known for his pedagogical insights 111Schwinger
even qualified Feynman diagrams as ”pedagogy, not physics” [27]. While
experiments are of course necessary, "creation” in Feynman’s mind might also
refer to the building of a predictive mathematical model, seen as the sine qua
non for understanding. This program is best described by Hopfield [28] :

> "The central idea was that the world is understandable, that you should be
> able to take anything apart, understand the relationships between its
> constituents, do experiments, and on that basis be able to develop a
> quantitative understanding of its behavior. Physics was a point of view that
> the world around us is, with effort, ingenuity, and adequate resources,
> understandable in a predictive and reasonably quantitative fashion. Being a
> physicist is a dedication to a quest for this kind of understanding."

This tutorial aims at introducing such a quantitative understanding of
somitogenesis. Excellent reviews have been recently written on a more
biological side, e.g. [29, 30, 1], or on developmental oscillations in
general, [31], see also [32] for a review of synchronization in the present
context. We hope to provide here a modern mathematical introduction to the
field of somitogenesis, allowing for conceptual discussions framed with non-
linear models, in a language amenable to physicists. We will be careful to
relate models to experimental biology as much as we can.

In the following, we first briefly summarize the main biological concepts and
molecular players. The field is still evolving and new aspects are still
discovered to this date (we write those words in 2023). This justifies a more
theoretical and conceptual discussion. We then follow approximately a
chronological approach, describing how the (theoretical) understanding of the
field has progressed with time. Importantly, classical models proposed before
the molecular biological era have been crucial to suggest experiments and
ideas, and our ambition is to describe them in detail because they are still
relevant today, at the very least to frame the theoretical discussion.

The field has then been strongly driven by the constant experimental progress
in molecular biology, genetics, imaging, and more recently synthetic biology,
allowing scientitsts to explore more complex and refined scenarios, that we
will describe. Many of the most recent ideas described in the following also
find their origin in the era of "systems biology", with a focus on the
(emergent) properties of gene regulatory networks [33]. For this reason, there
is a bias in both experimental and modeling works, towards the signaling
aspects of the system, which we would loosely define as the dynamics of gene
expression in time and space, described by non-linear models. We will discuss
the experimental reasons why such an approach makes sense in retrospect, but
will also describe works exploring other aspects (e.g. mechanics). We will
eventually connect those models to current descriptions grounded in dynamical
systems or catastrophe theory [34], with the hope to infer some general
principles and scenarios [35, 36] (see e.g. summary Figure 1). In Appendix A,
we put together a condensed discussion of classical results on non-linear
oscillators and bifurcations, with examples relevant to the present context
(phase oscillator, phase responses, and some introduction to relaxation
oscillators and excitability). Appendix B contains calculations associated
with the main text. We also include multiple Jupyter notebooks to simulate the
multiple models presented in this tutorial :
https://github.com/prfrancois/somitetutorial

Figure 1: Sketches of Malpighi illustrating the process of somite formation in
chick embryos

#### Conventions and definitions used in tutorial

In Table 1, we summarize a couple of notations used throughout this tutorial

Symbol | Definition
---|---
$\theta_{i}(t)$ |  phase of a given oscillator at time $t$ and discrete position $i$.
$\theta(x,t)$ |  phase of a given oscillator at time $t$ and position $x$ (continuous limit of $\theta_{i}$)
$\omega_{i}$ |  frequency of an oscillator at position $i$
$\omega(x)$ |  continuous limit of $\omega_{i}$
$\Omega$ | (global) frequency of the segment formation process
$T$ | period of the segment formation process
$\phi$ | phase of an oscillator in a moving frame of reference
$\psi$ |  relative phase of an oscillator with respect to a reference oscillator (usually $\phi$)
$v$ | speed of propagation of the front
$S$ | size of somites
$L$ | size of the tissue (e.g. presomitic mesoderm)
$\tau$ | delay in the differential equations and/or the coupling
$\lambda$ | wave length of the pattern
Table 1: Some notations used in this review

When discussing biology, we follow a standard convention where specific gene
names are italicized (e.g. Mesp2, Lfng, names of pathways or gene families are
kept in normal font (e.g. Notch pathway).

For many theoretical works, we represent the spatio-temporal behaviour of the
system by so-called "kymographs", which are pictures showing the spatio-
temporal values of a variable as different colour/gray level. We will follow
the convention for representing kymographs from [17] and other works: columns
of the kymographs correspond to different times (with time increasing from
left to right), and lines to different positions in the embryo (with most
anterior cells on the top and tail bud cells on the bottom), see Fig. 2 below.
For models including growth, instead of imposing some moving boundary
condition, it is common to typically extend the tail of the embryos as a
fictitious extended region in space with a homogeneous pattern of expression,
as is represented at the bottom of Fig.2.

Figure 2: Correspondance between the observed spatial-temporal pattern of
genetic expression (in an embryo, top) and a theoretical Kymograph (bottom).
When representing the behaviour of theoretical models, as a convention, we
extend the tail expression pattern to a fictitious region posterior to tail
(dotted triangle at the bottom)

##  Chapter 1 Characterizing vertebrate segmentation : clock, waves,
morphogens

###  1 Early concepts

####  1 Vertebrate segmentation

One of the first recorded observations of somite formation is due to Marcello
Malpighi, a medical doctor who pioneered the use of the microscope for
scientific observation. In Opera Omnia [37], published in 1687, Malpighi drew
several stages of chick embryonic development, Fig. 1 (reproduced from [38]).
Somites were represented as balls of cells on both sides of the neural tube.
For the first time, it was visible from these drawings that somite formation
is a dynamic process, where somites sequentially form from anterior to
posterior as the embryo is elongating.

It took a few more centuries to get a more detailed view of embryogenesis and
of somite dynamical formation. In 1850, Remak observed that future vertebrae
arise from the fusion of the posterior part of a somite with the anterior part
of the following one [39, 1], suggesting that somites are not homogeneous and
present functional anteroposterior (another biological term for this being
’rostral-caudal’) polarity. Fast forward another century, a more precise
description of somitogenesis (the "genesis" of somites) was made possible with
progress in the manipulation of chicken and amphibian embryos, and was
motivated in parts by theoretical questions. We refer to Pourquié’s recent
review [1] for a very detailed description and now proceed in describing key
theoretical proposals of those pioneering times.

####  2 Morphogens

Turing’s seminal work on "The Chemical basis of morphogenesis" [40] represents
a conceptual turning point in theoretical embryology, . Turing introduced
several key ideas that have deeply shaped the entire field of developmental
biology, up to this day. In particular :

  * •

Turing suggested that some morphogenetic events find their origin in
differences of concentrations of chemical substances. While he explicitly
discussed in the introduction the role of mechanics in morphogenesis, he was
the first to consider a model where the chemical and mechanical aspects can be
separated.

  * •

Chemical substances driving development are called "morphogens", a term now
widely used in biology. Turing postulated that morphogens interact with each
other via reaction and diffusion. This can give rise to patterns (now
generally called "Turing patterns") at the origin of biological shapes.

The typical Turing ’interaction network’ is made of two morphogens : one
’activator’ morphogen, that diffuses slowly (thus with short range activity),
self-activating and activating a ’repressor’ morphogen, that diffuses rapidly
(thus with long range activity). A simulation of 1 D Turing mechanism with
(almost) homogeneous initial condition indeed gives rise to a periodic
pattern, where islands of the activator morphogens are limited by more broadly
expressed repressors. Turing patterns thus are a natural candidate for the
formation of metameric units, similar to the ones observed in vertebrate
segmentation. Alternation of stripes in a Turing model could either correspond
to a somite/nonsomite pattern or the anterior/posterior parts of somites.

Diffusion is crucial in the establishment and maintenance of Turing patterns,
for instance if a physical barrier is put in place, the long-range repression
effect is impinged, and new activating regions can emerge. Another key feature
of Turing patterns is their intrinsic length scale, which is a function of the
parameters such as diffusion constant. This led to a direct experimental test
of a Turing-based model of somitogenesis by Waddington and Deuchar [41] (and
later Cooke [42]), who generated amphibian embryos of different sizes by
adding/removing tissue at the gastrula stage. They observed that somite size
is scaling accordingly, i.e. bigger embryos have bigger somites in all
dimensions. This excludes a process where the length scale is set by a simple
reaction-diffusion process. Another difference can be found in the dynamical
aspect of the process: as said above, the formation of somites is sequential,
from anterior to posterior, while stripes or spots in a Turing system a priori
form simultaneously.

####  3 Positional information

Further considerations of the scaling of structures in embryos of different
sizes led to many conceptual discussions on how genetically identical cells
can take different fates, which are worth mentioning to better understand the
current theoretical framework. In 1969, Lewis Wolpert introduced the notion of
positional information in development [43]. Information here should be
understood in the colloquial sense: positional information is more akin to a
zip code or an address (rather than a physics-inspired definition of
information in relation to entropy). Wolpert’s underlying idea is that cells
have ways to "know" (or to compute) their position within an embryo, and to
differentiate accordingly. Then problems such as embryonic scaling boil down
to the problem of specification of positional information (which should
actively scale with, e.g., cell number).

Figure 1: The French Flag Model (Adapted from [44]). A graded morphogen
concentration is used as an input for cells to define three domains - via two
thresholds of activity $\Theta_{B,R}$. Those three domains are depticted using
the three colours of the French Flag. (left). If the embryo size is reduced
but the morphogen gradient properly scales, this ensures a scaled pattern of
cellular fates even if the number of cells itself is not conserved (right)

The paradigmatic example of positional information in biology is Wolpert’s
famous French Flag Model [44], Fig. 1. Imagine an embryo as a line of cells
(with the position of a given cell defined by its coordinate $x$), and imagine
that there is a graded concentration of a morphogenetic protein (let us assume
it is exponential of the form $C(x)=C_{0}e^{-x/L}$ where $L$ is the size of
the tissue to pattern). Then, cells have access to local concentration $C(x)$
and can decide their fate based on this. For instance, imagine that there are
two thresholds respectively at $\Theta_{B}$ and $\Theta_{R}$, then cells
observing concentration lower than $\Theta_{R}$ can develop into a "red" fate,
cells with observing concentrations between $\Theta_{R}$ and $\Theta_{B}$
develop into a "white" fate and cells with concentrations higher than
$\Theta_{B}$ can develop into a "blue" fate, giving rise to a paradigmatic
French Flag picture Fig. 1. The French Flag paradigm provides a parsimonious
explanation of embryonic scaling. If the number of cells is changing, one can
possibly scale patterning within an embryo by scaling the morphogen gradient
itself, which is arguably a much simpler problem to solve (both for biology
and for theorists). For instance, Crick proposed in 1970 a "source-sink" model
where a gradient of a diffusing protein is maintained at concentration $C_{0}$
at one extremity of the embryo and at $0$ at the other extremity [45]. A
solution of the 1D diffusion equation with those boundary conditions clearly
is a linear, steady-state profile, which thus naturally scales with the size
of the diffusing field. To ensure scaling, one simply needs to impose boundary
conditions, which is consistent with the existence of embryonic regions such
as organizers [44]. Such ideas led to multiple discussions on the
theory/conceptual side. For instance, it is not clear if one can separate any
informational content from the processing of this information. Some of those
early debates are summarized by Cooke [46], who observed that the proportional
allocation of cells to different tissues in embryos of vastly different sizes
can not be very easily explained with simple morphogen gradients or reaction-
diffusion models . He suggested some coupling between protein production rates
and the size of tissue might rather play a role as a ’proportion sensor’. It
should be mentioned that our understanding of such scaling properties remains
incomplete to this date.

Coming back to segmentation, a natural idea within the positional information
framework would be to assume that different thresholds of one or several
morphogens would define somite locations. The problem is that many animals
(snakes, centipedes) can have many segments (more than 200 vertebrae in
snakes). In a French Flag/positional information picture, the potential number
of thresholds needed to explain somite formation appears unlikely huge [15].
Another issue is that from one animal to the other, there is some variability
in the number of somites even within the same species, which implies a degree
of versatility with respect to the overall body plan in the encoding of somite
position [15]. Other explanations are thus needed both for the process of
segmentation itself and the underlying scaling mechanisms.

###  2 Statics and dynamics of metazoan segmentation

####  1 Establishment of the (fly) segmentation paradigm

In parallel, starting in the early 1980s, molecular details of developmental
processes in general have been established and refined with increasing
progress in molecular biology, genetics, and, later on, imaging. The fruit fly
(Drosophila) model organism is the first organism for which the key principles
of segmentation and associated genes have been identified, starting in 1981
with a groundbreaking series of papers by Christiane Nüsslein-Volhard and Eric
Wieschaus [47] (who were awarded the Medecine Nobel Prize for this work in
1995.)

In a nutshell, fly segmentation appears, maybe surprisingly, largely
consistent with the "French Flag model" view, [44], Fig. 2). Multiple
morphogenetic gradients were discovered over the years: the bicoid gradient
defines identities in the anterior part of the embryo, while posterior
gradients such as nanos and caudal define identities in the posterior part of
the embryo [48, 49]. Those gradients are generally called "maternal, since
they are initially defined by localization of RNA molecules in the egg by the
mother (and subsequent cross-regulation, e.g. caudal translation is repressed
by bcd ).

In their original papers, Wieschaus and Nüsslein-Volhard identify so-called
"gap-like" phenotypes, in which mutants have parts of their body missing.
Those gap phenotypes are due to the mutation of so-called gap genes,
themselves normally expressed in the part of the body missing in the mutants.
Gap genes’s expressions are positioned and controlled by the maternal
gradients, and consistent with this, cellular identities can be shifted
anteriorly or posteriorly by changing the levels of the maternal gradients
[50].

Downstream the gap genes, we find pair-rule genes, then segmentation genes
[51, 49], Fig. 2. The pair-rule genes correspond to periodic structure every 2
segments, while segmentation genes are expressed in all segments. Those genes
are expressed in periodic stripes corresponding to future segments and their
sub-compartments. Such striped patterns naturally evoke reaction-diffusion
mechanisms to physicists, but quite astonishingly, it turns out that those
different stripes are encoded in the genetic sequence and regulated more or
less independently from one another. As an example, an Eve2 stripe genetic
module can be identified on the fly DNA, regulated by a subset of gap genes
independently from all other stripe modules [52], Fig. 2 B. Those discoveries
thus suggested a very local and feedforward view of development and positional
information, where concentrations of morphogens dictate local fates all the
way to segmentation genes. Remarkably, it has been shown since then that the
bicoid gradients and the gap genes downstream of it contain exactly the right
amount of information (in the physics sense) to encode identity with a single
cell resolutions along the entire fly embryo [53, 54, 55, 56].

Those discoveries considerably shaped the subsequent discussions on
segmentation in vertebrates as well. First, they firmly established the
morphogen gradient paradigm, where different levels define different
identities or properties. Second, they argue against models where reaction-
diffusion processes are crucial for robust patterning. The view coming from
fly segmentation is more local and modular: the definition of cellular fates
is done through a given gap gene combination [55, 56], which is specific to
the cell location, independently from all other locations within the embryo.
Consistent with this view, there is some variability in the pattern of gap
genes’ expression (and likely regulation) from one species to the other in
"long germ band" insects (forming their segments like flies) [57, 58], see
also [59] for simulations of underlying network evolution.

That said, it rapidly turned out that flies are to some extent evolutionary
exceptions. The almost paradigmatic morphogen, bicoid, does not exist outside
of Drosophila. Long germ segmentation further appears highly derived
evolutionary: it occurs in an egg of approximately fixed size, with
segmentation genes expressed more or less simultaneously, while in most other
metazoans, segmentation is sequentially coupled to embryonic growth [60] 111it
should be pointed out though that long germ segmentation still evolved many
times independently, suggesting deep evolutionary forces are at stake to move
towards such mode of segmentation  . Gap phenotypes are also not observed in
vertebrates, suggesting that segmentation is a more global, integrated process
in opposition to a more local process where identities are defined by local
morphogen concentrations. Finally, as said above, flies have a relatively
small number of segments compared to some vertebrates such as snakes.

Figure 2:  Summary of Fly segmentation (A) Schematic of the expression pattern
of some of the main genes regulating segmentation in Drosophila. Maternal
genes, control gap genes and gap genes in turn control the expression of pair-
rule genes. (B) A simplified, hierarchical feed-forward model for fly
segmentation, reproduced from [59]. The embryo is simulated as a one-
dimensional field. The left panel shows the behaviour of the model, the
maternal profiles are imposed and the network "computes" the concentration of
downstream genes from top to bottom. The right panel shows the topology of the
corresponding gene regulation network. Genes interaction are symbolized by
arrows, regular arrows correspond to activation, T-shaped arrows to
repression. For instance, one can see how individual Eve stripes are regulated
differently. We highlight Eve2, which is activated by Bcd and repressed by
both Anterior Gt and Kr in this model, and as a consequence appears at the
interface between those two genes.

####  2 Discovery and phenomenology of the segmentation clock

All animals are evolutionary related and, as a spectacular consequence, many
of the lower level controls of the animal physiology are similar even in very
different-looking animals [2]. This is especially patent for molecular
controls of embryonic development : many developmental genes are highly
conserved, and play the exact same role in many animals. A spectacular example
are Hox genes, which prescribe anterior-posterior identities of cells in
similar ways in all animals, to the point that Hox genes were proposed as a
’defining character of the kingdom Animalia’ [61, 62]. Coming back to
segmentation, given the crucial role of pair-rule genes in fly, several groups
then proceeded to identify and study their homologs in vertebrates.

It quickly appeared that vertebrate proteins closely related to the fly hairy
genes presented patterns in developing vertebrate embryos somehow reminiscent
of what happens in the fly. For instance, her1 in zebrafish 222her stands for
’hairy-E(spl)’ which is the name of the broader family of these proteins was
first described to present patterns, with broad stripes in the presomitic
mesoderm (PSM) and narrower stripes in somite primordia [63]. In 1997,
Palmeirim et al.identified a homologous of hairy in chick (called  c-hairy),
and carefully studied its behavior in a seminal work [5] redefining the entire
field.

Palmeirim et al.proceeded to study the pattern of genetic expression of
c-hairy. Comparing embryos to embryos, they confirmed that c-hairy presents
two distinct patterns of expression. In the anterior part of the embryos,
c-hairy is expressed in the posterior half of formed somites. But the pattern
of gene expression in the non-segmented pre-somitic mesoderm (i.e. posterior
to formed somites) appears much more complex. Depending on the embryo, c-hairy
is expressed broadly in the posterior, or into increasingly narrower and more
anterior stripes of genetic expression, not unlike what happens in zebrafish
for her1 [63].

The "Eureka" aspect of this work was to realize that this pattern of gene
expression in the posterior actually corresponds to snapshots of the dynamics
of a propagating (and narrowing) wave of c-hairy expression from posterior to
anterior, which appears clearly when embryos are reordered as a function of a
pseudo-time (see schematic in Fig. 3 A, c-hairy would correspond to the green
colour). To unambiguously show that such a wave originates from a posterior
oscillator, Palmeirim et al used an ingenious trick of chick embryology. They
cut the embryo into two pieces, fixed one side of the embryo, then waited
before fixing the other side. Assuming the dynamics on either side of the
embryo are independent of what happens on the other side, this allows the
capture of two time-points of the same dynamical process (essentially a two-
point kymograph), and from there to reconstruct the entire process using
multiple embryos. This technique indeed shows that the variability of the
c-hairy pattern comes from a dynamical gene expression, since in the very same
embryo one effectively sees a stripe of c-hairy gene expression move towards
the anterior, similar to what is observed for the gene expressions in Fig. 3
A). Furthermore, fixing the two halves of embryos with a time difference of 90
mins, one sees the same pattern of gene expression of c-hairy, but with one
extra somite on the right side vs the left (compare first and last time in
Fig. 3 A). This indicates that a periodic mechanism drives the waves of
genetic expression and is indeed correlated to somite formation, as expected
from the oscillatory models proposed previously (see Section Early models).
The very same technique of fixing one half of the embryo while keeping the
other alive was later used to show the existence of a segmentation oscillator
in Tribolium [64].

Figure 3:  Segmentation in vertebrate. (A) Phenomenology of segmentation,
coupled to growth. Oscillatory gene expressions in the posterior give rise to
kinematic waves in the embryo propagating from posterior to anterior. The wave
refines and stabilizes in the anterior to structure future somites into
anterior and posterior compartments. There is a differentiation zone
corresponding to the region where cellular oscillations stop and physical
boundaries form. On the right, we show schematically the oscillation phases
for a cell staying in the posterior (bottom) and for a cell (star) reaching
the differentiation zone at the end of the time window depicted. The
oscillation is slowing down as cells get more anterior, giving rise to a
period gradient within the oscillatory zone. (B) Experimental visualization of
the segmentation clock in a mouse embryo, adapted from [17]. The Middle panel
shows snapshots of a movie with a dynamical  Lfng reporter (see Section The
molecular forest), with recapitulation of the location of the last "stripe"
for each wave. The right panel shows a corresponding Kymograph . (C)
Experimental visualization of the segmentation clock in zebrafish embryos,
with a single cell resolution, adapted from [7]. A live reporter for  Her1 is
used. (D) Her-1 oscillation for single cells in a zebrafish embryo, reproduced
from [13]. The cell is moving (relatively) from tail bud until a somite where
it is integrated and where the oscillation stops. (E) Inferred period (top)
for two single-cell oscillators in a zebrafish embryo as a function of time,
reproduced from [13]. The distance measured from the posterior (bottom) allows
to identify the relative positions of the cells within the presomitic mesoderm
(PSM), and to correlate it to their respective period. The cell staying close
to the tail bud oscillates with a constant period, while the period is
increasing as a function of time for cell moving towards the anterior,
indicative of the existence of a period gradient within the PSM.

###  3 The segmentation clock paradigm

####  1 Phenomenology of the segmentation clock

It is now generally acknowledged that the work of Palmeirim et al.showed the
existence of what is now called the "segmentation clock". In this review, by
"segmentation clock", we mean the ensemble of periodic gene expressions, at
the embryo level, which controls the periodic formation of somites. Before we
focus on molecular details in the next section, we wish to point out four
high-level components and properties underlying the segmentation clock, which
will be central to the discussion in this review. The segmentation clock
paradigm is summarized in Fig. 3 A, with experimental illustrations in
subsequent panels (B-E).

Firstly, the segmentation clock emerges through cellular oscillators, clearly
visible in Fig. 3 B-C. Cells in the presomitic mesoderm PSM display
coordinated oscillation of multiple genes, thus defining a global oscillator
at the PSM level. Importantly, cellular oscillators are synchronized but not
in phase: waves of oscillations sweep the embryo from posterior to anterior,
as first evidenced in the work of Parlmeirm et al., and can now be seen using
real-time reporters Fig. 3 B-C. Those waves are related to the fact that, as
cells get more anterior, the period of their internal oscillator is increasing
(see e.g. the oscillation in the starred cell compared to posterior
oscillation in the schematic in Fig. 3 A, and see experimental measurements of
the period in single cells in Fig. 3 D-E). There are thus parallel anterior-
posterior period and phase gradients in the PSM. One of the key theoretical
question is to figure out how those gradients are related : do cells modify
their intrinsic period (slow down) so that a phase gradient builds up, or is
there a phase gradient building up (e.g. via cell-to-cell interactions)
leading up to an apparent period slowing down ?

As the waves of genetic expression move towards the anterior, and as the local
period of the oscillators increases, the wavelength decreases, before
stabilizing into a fixed pattern. Some genes, like c-hairy discussed in the
previous section, then form a stripe pattern of genetic expressions, localized
in half a somite. The formation of those stripes appears to be tightly coupled
to the formation of a somite boundary, Fig. 3 B, middle panel. Somites
eventually display an anterior-posterior (or rostral-caudal) pattern of
genetic expression, with some specific genes expressed in the anterior half of
the somite, and some others in the posterior half of the somite, see Fig. 3 A
–within the same somite, blue gene is rostral, and green gene is caudal.
Notice that this pattern is to some extent reminiscent of pair-rule patterning
in flies, compare Fig. 3 A with Eve and Ftz in Fig. 2 . The region where
cellular oscillations stop and where, subsequently, boundaries form between
future somites, is labeled as ’differentiation zone’ in Fig. 3, and is the
second important component of the segmentation process. Specific genes are
expressed in this region. Very often in the literature, this region is
designated not as a zone, but phenomenologically reduced to a single front,
often called ’wavefront’, largely because of the initial Clock and Wavefront
model that we describe in the section Early models. Notice that the slowing
down of the cellular oscillations is tied to stable patterns in somites,
following posterior to anterior waves of genes such c-hairy. So clock and
differentiation front might not be considered as independent processes. They
seem at the very least coordinated, which raises the fundamental question of
the nature of the front and its spatial extension, a central question
discussed in this review (see Fig. 1 for a synthesis).

Thirdly, segmentation is tied to embryonic growth. Schematically, as the tail
is growing, cells move anteriorly relative to the growth zone (Fig. 3 E
bottom), so that, as said above, a phase gradient is accumulating and their
period appears to increase (Fig. 3 E top). They eventually differentiate and
integrate into somites. It is well established that embryonic growth is
connected to anterior-posterior gradients of various morphogens, and thus it
is natural to think that those gradients likely regulate somite formation in
some way, especially in line with the French Flag and the Fly paradigms where
anterior to posterior gradients largely control segment position. Since
somitogenesis is a much more dynamical process, there are two additional
questions: how do gradients control cellular oscillators themselves (e.g.
their period and amplitude ?), and how do they control the location of the
differentiation zone? Again those questions are not independent and we will
comment on them in this review.

Fourthly, vertebrate segmentation is a tissue autonomous process: interruption
of continuity of the presomitic mesoderm (PSM) - the undifferentiated tissue
from which somites derive - does not impinge somite formation. Furthermore,
local inversion of fragments within the PSM leads to an "inversion" of the
progression of somite formation. This suggests that once cells exit the tail
bud, they are largely preprogrammed to oscillate and eventually differentiate
in a precise way, and as we will see below it seems that indeed dissociated
cells behave very similarly to cells within the embryo, suggesting that many
processes are largely cell autonomous. From the theoretical standpoint, it is
not clear how this large degreee of cell autonomy eventually gives rise to
weill proportioned, multi-cellular somites.

To finish this section, it is important to point out that the existence of an
oscillator (or clock) driving the formation of somites was first predicted and
studied by Cooke/Zeeman and Meinhardt in two pioneering models, that we
describe in details in section Early models. This is a nice example in biology
where theory was far ahead of experimental biology and inspired it.

####  2 The molecular forest

The phenomenology of the segmentation waves first described in [5] and
summarized in the previous section has been confirmed and generalized to other
model organisms. Furthermore, it has been established in subsequent works that
not only the phenomenon of oscillations and waves is broadly observed, but
also that a plethora of genes is oscillating, forming multiple parallel waves
of gene expression during vertebrate segmentation [65, 66]. Listing here all
phenotypes and interactions discovered would be both impossible and
potentially confusing, but to understand the principles underlying current
modeling, it is important to summarize some of the biological players, as well
as some crucial biological mechanisms they have been suggested to regulate. It
should be pointed out that a major difficulty is that interactions are not
conserved between different species [67], e.g. a gene oscillating in one
species might not oscillate in another one. This renders the study of
molecular segmentation clock very difficult, and to this date, no clear
conserved molecular mechanism controlling the segmentation oscillator has been
established, and in fact, segmentation waves likely work in slightly different
ways in different organisms (see section Difference between species). We
summarize some important results in this section, with a special focus on
mouse somitogenesis, but will also comment on some results on other animals.

Three major signaling pathways have been implicated in the segmentation waves:
Notch, Wnt, and FGF [66]. The current consensus is that the core oscillator is
related to the Notch signaling pathway, implicated in cellular communication
[32]. Notch ligands (called deltas) are produced and membrane-bound at the
surface of cells, and interact with Notch receptors at the surface of
neighboring cells, driving transcriptional response. Lunatic Fringe (Lfng), a
glycotransferase modifying Notch activity, is at the heart of the chick
segmentation clock [68]. Misexpression of Lfng disrupts somite formation and
anteroposterior compartmentalization in chick [68], and similar phenotypes are
observed in mouse [69, 70]. Lfng does not oscillate in zebrafish though, and
studies in this organism have rather focused on other components of Notch
signaling pathway. Notch ligands (delta) are implicated in many segmentation
phenotypes. Perturbation of Notch signaling results in clear somite formation
defects [11]. Mutations of delta ligands do not prevent segmentation but
impact the coherence of segmentation waves, prompting the suggestion that the
main role of Notch signaling is to synchronize cellular oscillators [71, 72].
Indeed, real-time monitoring has since then confirmed that in delta mutants,
individual cells oscillate but are desynchronized [7]. Lfng has actually been
shown to play a role in this synchronization as well in mouse by modulating
delta ligand activity and thus Notch signaling in neighboring cells [73]. The
Hes/Her transcription factors, phylogenetically related to the fly hairy gene
mentioned above, appear to play a major role in the core part of the
oscillator[74, 75, 76]. Interestingly, serum-induced oscillations of Hes1 (a
Notch effector) are observed in multiple types of cultured cells (myoblasts,
fibroblasts, neuroblastoma) with a 2-hour period consistent with somitogenesis
period in several organisms [77], suggesting that it could be part of a more
general core oscillator based on a negative feedback loop [78]. Hes5
oscillations have also been implicated in neurogenesis [79]

Another major oscillating pathway is Wnt. Axin2, a negative regulator of the
Wnt pathway oscillates in mouse, even when Notch signaling is impaired [80].
Perturbation of Wnt signaling pathway results in segmentation phenotypes, e.g.
Wnt3a is required for oscillating Notch signaling activity. Importantly, a
posterior to anterior gradient of $\beta$-catenin (a key intracellular
mediator of Wnt transcription) is also observed [6], and crucially, mutants
with constitutive (i.e. highly expressed) $\beta$-catenin display non-stopping
traveling waves of gene expression within the PSM, suggesting that Wnt plays a
crucial role in the stopping of the segmentation waves. However, Wnt does not
oscillate in zebrafish

The last major player is FGF. Many genes related to the FGF pathway oscillate
[65], but the major feature of FGF is that it appears to control the location
and the size of somite. FGF8 presents a graded expression, from posterior to
anterior [81, 10]. FGF8 overexpression disrupts segmentation by maintaining
cells in a posterior-like state (characterized by the expression of many
characteristic markers and associated posterior morphology). Dubrulle et
al.used beads soaked with FGF8 to show that local overexpression of FGF leads
to strong segmentation phenotype in chick (monitored by looking at the
expressions of the Notch ligand c-delta) [81, 10]. If the bead is initially
placed in a posterior region, as elongation proceeds and the bead gets more
anterior, major changes are observed, with several small somites anterior to
the bead and one big somite posterior to the bead. If the bead is placed
midway in the PSM, a similar phenotype is observed but only around the bead,
up to a well-defined anterior boundary, 4 somites posterior to the first
somite boundary. Grafts of FGF beads in this region yield no phenotype.

Some genes are also (in)activated following an apparent front moving from
anterior to posterior, likely controlling somite formation. For instance, in
mouse, Tbx6 is expressed only in the oscillating PSM region[82]. Furthermore,
in the most anterior section of the presomitic mesoderm, segmentation
oscillators slow down, and genetic waves of expression either stabilize or
simply disappear. In the region where the system leaves the oscillatory
regimes, new genes are expressed, such as Mesp2. Mesp2 is first expressed in a
few broad stripes, possibly slightly bigger than a somite size, before
restricting itself to the anterior part of the somite [83, 82]. Mesp2
activates Ripply2, which then turns off Tbx6.

Somites present Anterior-Posterior (or rostrocaudal) polarity. As said above,
within a somite, Mesp2 is eventually becoming anterior (A) within a somite.
Other Notch signaling pathway genes get stably expressed in the posterior part
(P) of somites, such as Dll1 or Uncx4.1. [82]. Interestingly, the boundary
formation between somites is clearly correlated to the Posterior-Anterior
boundary between Notch signaling in the posterior part of a future somite and
Mesp2 in the anterior part of the next one [84, 85].

One issue, first discussed by Meinhardt [86] is the problem of the symmetry of
AP vs PA boundary to define the somite boundary. This is visible on kymographs
such as the one in Fig. 2 focusing only on the expression of oscillator genes
: the boundaries at steady state between the green and the blue region do not
distinguish between internal or external somite boundaries. Meinhardt
suggested that there might be a third state (X) to define the such boundary.
Experiments in zebrafish possibly falsify the existence of such intermediate
state: mutants for convergence extension 333a process of cellular convergence
towards an axis, so that, because of volume conservation, tissue is thinning
perpendicular to the axis and extending in the direction of the axis  give
rise to broad, large somites with well-defined boundaries, but only 2-cell
wide in the anteroposterior directions [87]. So in such somites, there can not
be any cell corresponding to a hypothetical X state [A possible caveat is that
those cells are polarized so that there could be subcellular divisions
allowing for the existence of the X state]. Coming back to mouse, in [85], a
solution is suggested where the clock would in fact impose a rostrocaudal
gradient of Mesp2 inside the somite, imposing a natural polarity, where the PA
border between somites is "sharper" than the AP border within somites, leading
to a local "sawtooth" pattern. This exactly fits the pattern of downstream
genes implicated in cellular adhesion [88].

Figure 4: Schematic of some key molecular players in somitogenesis, using
mouse genes as examples. Anterior is on the top, posterior at the bottom. In
mouse, there is only one wave (i.e. roughly a $2\pi$ phase shift) of genetic
expression within the presomitic mesoderm.

It is worth mentioning at this stage a few other higher-order molecular
controls modulating somitogenesis formation. Retinoic Acid (RA) is well-known
to form an anteroposterior gradient opposite to FGF in metazoan embryos. RA
mutants display smaller somites [89]. So a natural question is the impact of
RA mutation on FGF gradient and the segmentation clock, [90, 91].
Surprisingly, embryos deprived of retinoic acid form asymmetrical left-right
somites. The associated phenotype is highly dynamic: for the first 7 somites,
Lfng and Hes7 waves are symmetrical, but afterward somites on the right side
of the embryo form later than on the left side, with one to three cycle delay.
The wave pattern is asymmetrical, and Mesp2 is more anterior on the right
side. This somite asymmetry is a consequence of the general left-right, Nodal
induced asymmetry (driving in particular internal organs asymmetry) [92, 90],
so that RA appears in fact to act as a buffer of this already present
asymmetry.

There are also many interesting modulations on the formation of the somite
boundaries. For instance, it is possible to induce separation between the
rostral and caudal parts of a somite by modifying cadherin and cad11 [93],
thus reavealling a length scale half of somite size in mouse. Conversely, in
zebrafish, disruption of her7 creates somites with alternating weak and strong
boundaries, suggesting the system can also generate an intrinsic length scale
twice the somite size [94].

####  3 Visualizing oscillations in embryos

Recent years have seen the development of multiple fluorescent reporters,
allowing for the real-time observations of some of the clock components. In
mouse, the current toolbox includes reporters for Notch signaling pathway,
such as a destabilized luciferase reporter for Hes1 [9], destabilized Venus
reported for Lfng (LuVeLu) [6]. An Axin2 reporter associated with the Wnt
signalling pathway is also available [18] as well as Mesp2 and FGF Erk
reporters [95]. In zebrafish, reporters for the Notch signaling pathway are
available as well, mostly based on Her1 fluorescent fusion proteins, and a
single cell resolution to visualize oscillations has been achieved [7, 96, 97,
98]. It should be pointed out that it is not necessarily easy to combine
reporters to visualize multiple components of the system in real-time, one
reason being that some of them are based on similar fluorescent proteins and
would not be easily distinguishable in the same cells [18].

Oscillations of Notch signaling pathway in single cells present a
characteristic profile, where both the average and the amplitude of the
oscillations increase as cells mature towards the anterior PSM. In zebrafish,
the last peak-to-peak time difference is approximately twice the period in the
tailbud [97], consistent with the strong slowing down first inferred from in
situs [76]. Waves of oscillations move from posterior to anterior to the very
anterior PSM, so that the most anterior cells within a somite are the last
ones to stop oscillating (as measured by the timing of the last peak of
oscillation [97]). This contrasts with the idea of a differentiation front
moving continuously from anterior to posterior: there, within a future
presumptive somite, anterior cells are expected to differentiate (and stop
their clock) before posterior cells. Such a mechanism creates an asymmetry in
the wavefront, with a $\pi$ phase shift within a future presumptive somite,
giving a "sawtooth" pattern within the presumptive somite. This could define
anterior and posterior somite compartments [97], and relate to the previous
observation that the system can generate a length scale twice the normal
somite-size [94].

It is also possible to monitor mitotic cells in embryos, providing a natural
perturbation of the segmentation oscillator. Mitosis delays oscillation in
cells, but divided cells eventually resynchronize with their neighbors after
roughly one cycle [7]. Interestingly, sibling cells are statistically more
synchronized with one another than with their neighbors, which shows that
single-cell oscillations are rather robust and only modulated by interactions.
Lastly, there is a clear interaction between the cell cycle and the
segmentation oscillator, since mitosis happens preferentially at a well-
defined phase, when Notch activity is the lowest (which possibly provides a
natural mechanism for noise robustness in presence of equal partitioning of
proteins) [7]. In Notch pathway mutants, single cells still oscillate, but in
a desynchronized way and with a longer period. The amplitude of Notch
oscillations in mutants appears bigger than in WT, with possibly a modest
increase towards the anterior, but there is no obvious increase in period
length in those mutants.

####  4 Biomechanical aspects

When treated with Noggin (an inhibitor of another signaling pathway called
BMP), non-somite mesoderm spontaneously segregates into somite-like structures
[99]. Those have sizes similar to normal somites, and when grafted instead of
normal somites, express normal somite markers. Contrary to normal somites,
they form almost simultaneously without the need for a clock, and are not
linearly organized but rather look like "a bunch of grapes". Importantly, they
do not have well-defined rostrocaudal identities: rather, cells within those
somite-like structures display patchy expressions of rostral and caudal
markers. This suggests that normal anteroposterior patterning within somites
might in fact be one of the main outputs of the clock [100].

The biomechanical program responsible for somite segregation can thus be
triggered independently of the segmentation clock. This suggests that there is
a level of biomechanical self-organization in the system, with associated
length scales, which raises the question of the multiple scaling effects at
play and of downstream self-organization within a given somite [101].
Consistent with this, it has been recently shown in normal somitogenesis that
tension forces allow for a correction of initial left-right asymmetries in
somite size [102]. This possibly suggests an overall view where slightly
imprecise signaling mechanisms (clock, wavefront, somite anteroposterior
polarity) are later canalized/corrected/adjusted by downstream biophysical
processes, such as tissue mechanics [102].

####  5 Difference between species

While the phenomenology of somitogenesis is roughly conserved between species,
it is also worth pointing out rather striking quantitative and qualitative
differences.

The segmentation period varies widely between species: around 30 mins for
zebrafish, 90 minutes for chicken, 2 hours for mice, and 5 hours for humans
[103]. More direct comparisons between mammalian cells, have been done using
stem cell cultures differentiated into PSM cells (See the section Stem-cell
systems for more details) [104, 105]. Mouse and human cells were first
compared [104], and later on, a segmentation ‘zoo’ was designed, including
marmoset, rabbit, cattle, and white rhinoceros [105]. The segmentation clock
periods in this new zoo range from 150 mins in rabbit to 390 mins in marmoset,
and are comparable to the ones in embryos. ‘Swap’ cells where e.g. human
sequences for the Hes7 gene is introduced in mouse cells show a period
increase of 20 to 30 mins, so only a fraction of the 200 mins difference of
periods between the two species. This suggests that internal cellular
biochemistry (rather than specific coding sequences) plays a role in setting
up the segmentation period.

Those scaling dependencies appear rather specific to the segmentation clock
though: the authors estimate parameters for other genetic cascades and protein
degradation rates in mice vs humans, and, while degradation rates are slower
in human cells than in mice cells, the typical differences are at most by a
few tens of percents (while the segmentation period varies by more than two-
fold), and for some important mesodermal proteins like Brachyury (also called
T) there is hardly any difference at all. All in all, those experiments
suggest that the biochemical reactions specifically implicated in the
segmentation clocks are essentially scaled in one species vs another.
Interestingly, this scaling could be rather global in the sense that the
segmentation clock period scales with embryogenesis length (defined as the
time from fertilization to the end of organogenesis). Of note, similar scaling
of embryonic developmental steps is often observed, for instance, different
fly species living under different climates (and thus different temperatures)
present scaling developmental stages [106]. See more discussions on scaling in
the Appendix, section Scaling Laws.

Beyond the time scales of the segmentation period and development, it is worth
pointing out that the wave pattern observed in the PSM widely varies between
species, Fig. 5. In Mouse and Medaka, there is only one ‘wave’ of genetic
expression within the PSM (meaning that the oscillators close to the front are
less than one cycle phase-shifted compared to the oscillators in the tail
bud). In zebrafish, there are three waves, and in snake, there are 8 to 9
waves. This suggests that the relative clock period as a function of relative
position within PSM varies widely between species. While in mice, the period
close to the front is only slightly shorter than the period in the tailbud, in
other animals such as zebrafish and snakes, the relative period in the
anterior appears to be at least 3 times longer, possibly more [107, 108, 109].
Interestingly the period profile as a function of relative position within the
PSM is highly non-linear, almost diverging towards anterior PSM, and rather
identical between zebrafish and snake, see [108] for a comparison. This could
indicate some common mechanisms ensuring the coordination of the slowing down
of individual cellular oscillators.

Figure 5: Schematic of the different wave patterns in different species.
Adapted from [109, 8, 1].

It is proposed in [108] that the extensive number of segments in snake vs
other animals is indeed due to a relatively slower overall growth rate
compared to the segmentation clock. Imagine for instance a zebrafish growing
at half the normal rate, but with a segmentation clock keeping the same pace,
then it would naturally have twice as many segments. This scenario is
supported by the following back-of-the-envelope calculation :

  * •

assuming PSM growth is completely driven by the cell cycle, period
$T_{cycle}$, the number of generation times for the PSM to fully grow is
$n_{g}=T_{tot}/T_{cycle}$ where $T_{tot}$ is again the total developmental
time

  * •

the length of a somite approximately is $S=\alpha LT$, where $T$ is the period
of the segmentation clock and $\alpha=\ln 2/T_{cycle}$ is the growth rate of
the PSM ($\ln 2$ factor converts into time via cell division)

  * •

eliminating $T_{cycle}$ one gets $n_{g}=\frac{T_{tot}}{T}\frac{S}{L\ln
2}=n_{s}\frac{S}{L\ln 2}$ where $n_{s}$ is the number of somites (assuming a
constant period of the segmentation clock).

Now $n_{s}$ is 315 in snake and 65 in mouse, but $\frac{S}{L\ln 2}$, the
rescaled ratio of somite vs PMS is also 5 times lower in snake than in mouse,
so that both effects compensate and the number of generation $n_{g}$ is the
same independent of the organism. This suggests a picture where $n_{g}$ is
constant across species for other reasons, and that inter-species variability
in the number of stripes indeed primarily comes from different values of
$T/T_{tot}$ or similarly $T/T_{cycle}$. Notice that, if the segmentation clock
period gradient within the PSM is (once rescaled) the same in all species
irrespective of PSM size, then if a cell spends relatively more time (in cell
cycle units) to go from tail bud to the front compared to other species, it
accumulates much more extensive phase gradient, which results in more waves
within the PSM, consistent with what is seen in snake (see more detailed
calculations in Appendix Number of waves in a growth model).

Going into more molecular details, it turns out that there is quite some
variability/plasticity between species in the genes oscillating [67].
Microarrays [65] identify 40 to 100 oscillating genes in the PSM, mostly
involved in signaling and transcription. In mouse, genes in Notch, Wnt and FGF
pathways oscillate, but in zebrafish it seems only Notch pathway clearly
oscillates. Phase relations between pathways also appear to vary between
species. Interestingly, only Hes1 and Hes5 orthologs appear to oscillate in
the three species considered in [67] (mouse, zebrafish, and chick), meaning
that there is likely "very limited conservation of the individual cycling
genes observed", and consistent with the hypothesis that the Hes gene family
includes the "core" oscillator. Needless to say, those differences might
matter a lot when modeling the segmentation process. There could be big
differences between segmentation processes in different species, and for this
reason, it is all the more important to discuss, contrast and compare multiple
models. Also, since individual cycling genes are likely, not conserved, this
justifies more top-down approaches, focused on higher levels, that can
eventually be related to actual gene expressions, rather than bottom-up
approaches too closely tied to the molecular implementation in a given
species.

##  Chapter 2 Early models

We now review models of vertebrate segmentation spanning more than 40 years of
theoretical work. We start with two pioneering models proposed before the
discovery of the segmentation clock : the Cooke/Zeeman clock and wavefront
model, and the reaction-diffusion Meinhardt model. Those two models frame the
conceptual discussion and still inspire experiments to this date, but they are
also useful reference points for subsequent models. We also review in this
section a cell-cycle model, proposed shortly after the discovery of the
segmentation clock, to some extent as an alternative explanation and also
providing a slightly different viewpoint (see also review in [110]).

###  1 The clock and wavefront framework

In 1976, Cooke and Zeeman [15] proposed a "clock and wavefront" model for
somite formation to recapitulate many aspects known at that time. In a
nutshell, the model argues that a simple way to build a spatially periodic
pattern (e.g. vertebrae) is to imprint a spatial record of a time-periodic
signal (i.e. a clock).

####  1 Qualitative view : wavefront

Such imprint is done with the help of a moving variable coupling positional
information to developmental time :

> "There will thus be a rate gradient or timing gradient along these columns,
> and we shall assume a fixed monotonic relation (non necessary linear)
> between RATE of an intracellular evolution of development process, and local
> positional information value experienced by a cell at the time of setting
> that rate."

It is not difficult to imagine such a variable in the context of embryonic
development since in many metazoans, growth happens in the anterior to
posterior (AP) direction, with anterior cells laid down before posterior ones.
This is represented in Fig. 1 A : here we define it as the age of the embryo
when the cell is born and positioned, counted from the beginning of embryonic
growth (anterior cells have age $0$, posterior cells have higher age), so that
the "positional information value" is linear in the position. While they were
not known at the time of the Cooke and Zeeman publication, we know now that
Hox genes [61] encode a similar discretized version of such coordinates, and
are likely controlled by a more continuous variable [111, 112]. Notice that
if, for some reason, the growth rate is twice as small, cells laid at a given
distance from the head are twice ’older’ compared to a reference embryo, so
that the positional information value grows at a doubled rate in absolute unit
in space. Thus positional information naturally scales with embryo length,
Fig. 1 A right. This naturally solves the scaling problem mentioned in Section
Early concepts.

Cooke and Zeeman propose that such positional information variable could then
be used to set the time for future developmental transitions. A simple model
would be that a developmental process is triggered after a time proportional
to the positional information value defined in Fig. 1 A. Phenomenologically,
this results in what we would call today a timer [112, 113], where the time at
which the process happens at a given position is proportional to the relative
position along the A-P axis.

In such a case, one would observe a developmental wavefront, moving along the
anterior-posterior axis. Thus in this picture developmental time (when a cell
is positioned along the AP axis) defines positional information, later setting
the stage for a kinematic wave of developmental transition moving from
anterior to posterior. [Importantly from a physics standpoint, the term wave
does not refer to any oscillation here, but rather is, to quote Zeeman, the
"movement of a frontier separating two regions" [114], see Box for the
definition of primary and secondary waves]. Again, an important aspect of such
proposal is that the kinematic wave would move at a speed scaling with the
embryo size since a temporal coordinate related to growth is properly
positioned relatively to an embryo of any size, Fig. 1 A right, consistent
with experiments where the number of cells is artificially reduced[42].

####  2 Qualitative view : clock

However, such a kinematic wave moves smoothly from anterior to posterior,
while the aim is to define discrete units (somites). To induce such change,
Cooke and Zeeman propose to introduce a periodic variable or "clock". A simple
description of the mechanism is illustrated on Fig 1 B. Imagine there is a
global oscillator in the embryo, or at the very least that there are
synchronized oscillators so that

> [pre-somite cells] "are entrained and closely phase-organized (…) because of
> intercellular communication."

Now assume that the front is moving from head to tail with a speed $v$. The
assumption is that as the front moves, it interacts with the clock to switch
the local state of a cell from undifferentiated (not somite) to differentiated
(somite). Importantly, the timing of the transition depends on the phase of
the clock when the front passes, to ensure a synchronous commitment.

To fix ideas, let us assume that a segmental boundary is formed if and only if
the clock interacts with the wavefront at phase $\phi=\phi^{*}$ (Fig 1 B).
Then starting from an initial segmental boundary where the front is present
(phase $\phi=\phi_{*}$ at $x=x_{1}$), the clock goes on ticking (period $T$)
while the front is passing. No boundary is formed until the clock reaches
again the phase $\phi=\phi_{*}+2\pi$, i.e. after waiting for the period $T$.
During that time, the front has moved from position $x=x_{1}$ to position
$x_{2}=x_{1}+vT$, where the next segmental boundary is formed. This entire
process is then:

> "converting the course of the wavefront into a step function in time, in
> terms of the spread of recruitment of cells into post-catastrophe behavior."

It is thus clear that segments of size $S=vT$ are sequentially formed.
Importantly, this process recapitulates the minimum phenomenology of somite
formation. Somites form periodically in time, and sequentially in space.
Future somite boundaries are encoded in the tissue by the kinetics of the
wavefront and the clock, so way before boundaries form. Notice that as soon as
we assume the existence of a clock with period $T$ and of a wavefront of speed
$v$, the size of the pattern to be proportional to $S=vT$ by dimensional
analysis, irrespective of the details of the model, so that if the clock
period $T$ period is fixed, the size of the segment is proportional to $v$
(which should thus scale with embryonic size) . See Appendix section Scaling
Laws for discussions of other possible scaling laws.

Figure 1: Qualitative view of the Clock and wavefront model. (A) A temporal
coordinate is imposed on the embryo, via some monotonic process (e.g. growth),
defining positional information. Different colours indicate different values
of the temporal coordinate, notice that cells along the same anterior
posterior position have the same coordinate. Also if the embryo is smaller
(right), the temporal coordinate should scale with the size of the embryo. (B)
From top to bottom, one cycle of the segmentation clock (left), as the
wavefront (vertical dashed line) progresses with speed $v$ along the temporal
coordinate defined in A (right). Phase $\phi^{*}$ of the clock defines when
the new somite boundary is formed, here at position $x_{2}$

####  3 Mathematical model : Wavefront

Cooke and Zeeman’s paper is also groundbreaking because it uses seminal
mathematical notions to describe developmental transitions. The model is
inspired by catastrophe theory, a branch of applied mathematics concerned with
a systematic classification of qualitative changes in behaviors of dynamical
systems. There the state of a cell is defined as a vector in a
multidimensional space, which generally localizes on a small number of
attractor domains (defining different cell states). The idea is that cells
move smoothly within each attractor domain, but developmental transitions
occur when cells abruptly change their attractor domain (akin to a
"catastrophe" [115], see also the work of Zeeman [116, 34]). As pointed out in
[117], there is no explicit equations provided for their model, but their
exact reasoning can easily be put into equations, which we do in the
following.

Cooke and Zeeman graphically suggest in their Fig. 4 [116] that somite
formation is induced by a bistable/cusp catastrophe, and that space and time
define the two parameters controlling the transition. Calling $t$ the time,
$p$ the positional information (which should be related to the anteroposterior
position in the embryo, higher $p$ being more posterior), and $z$ the variable
representing the state of the cell, let us then define a potential :

$F(t,p,z)=z^{4}/4-pz^{2}/2+\mu tz$ (1)

This functional form is identical to the one generated by the so-called
"Zeeman Catastrophe Machine" [34] (see also section Generalization : Zeeman’s
primary and secondary waves). A cell at the local position and time $p,t$ has
a state variable $z$, driven by the landscape defined by Eq. 1 (Fig. 2). All
cells are independent and each cell has its own landscape and state variable
$z$; it is implicitly assumed here that a positive value of $z$ corresponds to
an undifferentiated state, while negative values correspond to a
differentiated (somite) state. For simplicity, we put time and space in
different monomials in Eq. 1, which is not generic, but we will comment on
more general forms below.

Figure 2: Representation of the two-well landscape depending on time and space
defined by Eq. 1. The somite state corresponds to the left well, and the
undifferentiated state to the right well

The equilibrium points are given by the solution of the third-order polynomial
equation $\frac{\partial F}{\partial z}=z^{3}-pz+\mu t=0$. Assuming the system
is such that it rapidly stabilizes, we first see that for
$t\rightarrow-\infty$, the system is in a "positive" $z\propto t^{1/3}$ state
(so corresponding to an undifferentiated state) while for
$t\rightarrow\infty$, the system is in a "negative" $z\propto-t^{1/3}$ state
(corresponding to a somite state). Using classical algebra, it is not
difficult to show that for $p<0$, the system is monostable, i.e. $z$ can only
take one stable positive value, so can not differentiate. The interesting
behaviour occurs for $p>0$, for which there is a bistable region (i.e. $z$ can
take two stable values), delimited by $p=\left(\frac{27}{4}\right)^{1/3}(|\mu
t|)^{2/3}$. The most interesting behavior occurs along the line
$p=\left(\frac{27}{4}\right)^{1/3}(\mu t)^{2/3}$, which corresponds to the
saddle-node bifurcation where the high $z$ (i.e. undifferentiated) state
disappears (this line corresponds to what Zeeman calls a "primary"
developmental wave in [114], see section Generalization : Zeeman’s primary and
secondary waves). Inverting the expression, and assuming the system quickly
relaxes to a steady state, at time $\mu
t_{c}(p)=\left(\frac{4}{27}\right)^{1/2}p^{3/2}$, the system at position $p$
has no other choice than to suddenly jump from the positive to the negative
state (Fig 3 A-B). Notice this jump happens (much) later for higher $p$.

In this view, there would be a kinematic differentiation front, continuously
moving at higher $p$ values as a function of time, which is what Cooke and
Zeeman refer to when they say the actual differentiation wavefront involves :

> "a kinematic ‘wave’ controlled, without ongoing cellular interaction, by a
> much earlier established timing gradient."

Cooke and Zeeman point out that such variable $p$ could be easily set up by a
smooth, anteroposterior (timing) gradient.

####  4 Mathematical model : Clock

To make a somite, we shall not need a smooth wave propagation, but rather a
simultaneous differentiation for a block of cells - for a range of different
positions in the embryo $p$. To account for such "block" differentiation, one
needs to introduce a clock. There are multiple ways to put that into
equations, but to fix ideas, let us thus consider the following addition to
the cusp catastrophe model :

$\dot{z}_{p}=-\frac{\partial F}{\partial z}-k\delta_{T}(t)=-\mu
t+pz_{p}-z_{p}^{3}-k\delta_{T}(t)$ (2)

Figure 3: Mathematical formulation of the original Clock and Wavefront model.
(A-D) The blue curve indicates the possible steady-state values of the state
variable $z$ from Eq. 2 as a function of space and time. The actual dynamics
of variable $z$ are sketched with an arrow. High $z$ corresponds to
undifferentiated cells, and low $z$ to somites. When $t$ is high enough the
system goes through a saddle-node bifurcation from a bistable to a monostable
system, and $z$ suddenly jumps from high to low value. In the absence of a
clock, this transition happens at a later time for more posterior cells
(compare A and B). (C-D) The effect of the clock (red arrow) is to
periodically lower $z$, so that cells close to the bifurcation will jump from
the high to low state branch. Ensemble of cells close enough to the
bifurcation jump at the same time, thus defining discrete bloks. This is
illustrate in Panel with a yellow line Figure 4: (A) Kymograph for $z$ in the
absence of the clock, Bistable and Monostable zones are indicated for
reference (B) In presence of the clock (red arrows), modeled as periodic kicks
uniform in space, blocks of cells are simultaneously induced from high to low
$z$ state, modeling somite commitment. Notice that somite commitment happens
below the bifurcation line of (A), which is indicated by a dotted while line,
thus corresponding to the wavefront (C) Effect of a slower clock. In this
case, some cells reach the bifurcation line before the next pulse of the
clock, so that the front follows the bifurcation line with the periodic
commitment of smaller blocks (D) Changing time and space dependencies of the
control parameters changes the shape of the bifurcation line and of the front.

where we consider the time evolution of the state $z_{p}(t)$ for a cell with
positional information $p$. Here, $\delta_{T}(t)$ is a function periodically
kicking the value of all $z$ (magnitude $k$) towards a more negative value. In
such a situation, for cells close enough to the jump (saddle-node
bifurcation), the periodic kicking might induce differentiation earlier than
$t_{c}(p)$ (Fig 3 C-D). In particular, following a tick of the clock, we
expect multiple cells close to bifurcation to jump simultaneously to the
negative $z$ state, defining a somite in Cooke and Zeeman’s view. More
posterior cells with higher positional information $p$ initially stay in the
high $z$ state, but as they get closer to the bifurcation they will eventually
jump. Notice that in physics terms, the differentiation timing exactly
corresponds to the first passage time from the right well to the left well in
the time-evolving landscape of Fig. 2, under the control of the clock
periodically kicking towards the left. A 3D plot in Fig. 3 E further
summarizes the overall dynamics in the spirit of Fig. 4 of the initial Cooke
and Zeeman paper [15].

####  5 Simulated Clock and Wavefront model

Fig 4 displays actual simulations of Eq. 2 under various conditions, see also
attached Notebook. Fig. 5 also illustrates what happens within a landscape
description (see also Supplementary Movie 1). The bistable/monostable regions
are illustrated in Fig 4 A by simulating the system without the clock. Fig 4 B
shows what happens with the clock, where blocks of cells jump in a coordinated
way as desired. Notice that the new somite boundary after each pulse is always
below the bifurcation line, i.e. in the absence of the clock, cells would be
committed later compared to a situation with the clock. Interestingly, there
is a balance between the position of the bifurcation line and the
period/strength of the signal induced by the clock, a situation not studied in
[15]. For instance, if the clock is either weaker or slower enough, it can
happen that some cells will reach the bifurcation line between two cycles of
the clocks, leading to a "jagged" front, 4 C. The intuition for this result is
simpler: in the limit of no clock, the cells only transition when they go
through the bifurcation, so if the clock is both slow and weak, only cells
very close to the bifurcation would periodically transition to the
differentiated state.

Cooke and Zeeman further comment on an interesting geometrical feature of the
wavefront: as can be clearly seen from Fig 4 , the front is not a straight
line, which means that the speed of the wavefront is not constant in the
coordinate defined by the positional information $p$. Here, the saddle-node
bifurcation happens for $p\propto t^{2/3}$, so we expect the speed of the
differentiation front (in units of positional information) to be proportional
to $t^{-1/3}$ as well, i.e. going to $0$. If positional information is
directly proportional to the actual position, this means that that boundary
$i$ is located at a position scaling as $(iT)^{2/3}$, and thus the size of a
somite $i$ would then be $S_{i}\propto i^{-1/3}T^{2/3}$, so that the size of
somites would go to $0$ as well. This could explain why somites can get
smaller during development. This scaling law comes from the fact that position
and time are in separate coefficients in the polynomial of $z_{p}$ in Eq.2, a
choice we made here for simplicity. A more generic model would be to mix time
and space dependency, e.g. we can add a temporal dependency in the linear term
$z_{p}$ that modulates the front speed and shape, see e.g. Fig. 4 D : the
speed front would then go to zero and a stable boundary would form separating
the monostable and the bistable region, thus leaving a permanently
undifferentiated region.

Figure 5: Time space dependency of the cellular states in the Landscape
defined by Eq. 1, with same conventions as Fig. 2. Different lines correspond
to different positions (top is more anterior), and different columns to
different times. Green beads correspond to the time evolution of the system
without the clock, and orange bead to time evolution with the clock as defined
by Eq. 2. A cycle of the clock is completed every three columns. The
background colour is a function of the state of the cell in the Clock and
Wavefront model (light green: undifferentiated, light blue: differentiated)

Lastly, it is worth mentioning that in Cooke and Zeeman’s view, the clock is
an external pacemaker, essentially independent from the catastrophe
controlling differentiation, and could go on oscillating with minimal impact,
even in differentiated cells. Remarkably, the clock has an effect on the state
of the cells only close to the primary wave defined by the saddle-node
bifurcation. There are important experimental consequences of this
observation: for instance, if one could find an external way to manipulate the
variable $z_{p}$, one could induce somite formation without a clock, for all
cells within the bistable zone. Conversely, one should be able to largely
manipulate features of the clock (such as the period) without impacting the
potential driving the dynamics of the variable $z_{p}$. The most direct way to
test this would be to change the clock period, to see how this impacts the
speed of the regression and the size of the somites. However, there could be
new features arising in a regime where the clock is very slow, or has only a
weak influence on $z_{p}$: as illustrated in Fig 4 C, one can obtain a mixed
system with both discrete and continuous jumps for weak or slow clocks.

This example illustrates one issue in defining the wavefront: depending on the
parameters, the jump in $z_{p}(x,t)$ can be discrete within a block of cells,
continuous, or both. Thus the actual wavefront of differentiation is an
emergent feature of the interactions of the system, that might not be easily
associated with some simple observable (e.g. a given level of a morphogen).
There is an even more general lesson here: processes that are independently
regulated (here the clock on the one hand and the possible states of the cell
$z_{p}$ on the other hand) might become more coupled close to a bifurcation
(i.e. at criticality [118]), with important phenotypical consequences. For
this reason, it might in fact be desirable that both the clock and the
kinematic wave induced by the $z$ jump are in fact coordinated upstream in
some way. For instance, one could imagine models where the ‘constant’ term in
the right-hand side of Eq. 2 could also depend more explicitly on $p$ and on
the phase of the clock, or we could imagine that the strength of the clock
increases with clock period to prevent a situation like Fig. 4 C. Conversely,
a weaker clock might in fact be desirable, for instance, the jagged line in 4
C could be used to define anteroposterior polarity within one somite, so again
requiring some level of fine-tuning or coupling between the clock and the
primary wave.

####  6 Generalization : Zeeman’s primary and secondary waves

The Clock and Wavefront model is related to an earlier proposal by Zeeman
regarding the existence of "primary" and "secondary" waves for spatially
extended dynamical systems [114]. Zeeman proposes a much more general theory,
with illustrations from epidemiology, ecology, and developmental biology.

The general idea is to consider the propagation of a boundary separating two
regions with different steady states.

> "By a wave, we mean the movement of a frontier separating two regions. We
> call the wave primary if the mechanism causing the wave depends upon space
> and time."

An example offered by Zeeman in the context of embryonic development is a
field of cells, where initially cells are in a B state, but where cells can
also exist in an A state because of bistability. A primary wave can then
propagate from a region of A cells into a region of B cells if cells lose
their ability to be in the B state. This can happen for instance via a saddle-
node bifurcation, say in response to a disappearing morphogen. For this
reason, in this review, we will associate primary waves with bifurcations and
will be slightly more generic by including bifurcations associated with the
disappearance of oscillating states. Secondary waves are defined as such

> "We call the wave secondary if it depends only upon time, in other words it
> is series of local events that occur at a fixed time delay after the passage
> of the primary wave."

For instance, in a pandemic context, a primary wave would consist in the
propagation of a disease in a population, while the secondary waves would
consist of the delayed appearance of symptoms. This example illustrates in
particular how the secondary wave might reveal the existence of a hidden
primary wave. Similarly, in biology, the actual differentiation of cells might
be a secondary wave following a primary wave directing cells to go to
different fates depending on positional information depending on space, and
time.

To fix ideas and be more quantitative, let us consider a slightly more general
potential than Eq. 1, similar to the example that Zeeman uses in Fig. 5 of
[114]

$F_{\epsilon,\alpha}(t,p,z)=\epsilon(z^{4}/4-(p+\alpha t)z^{2}/2+\mu tz)$ (3)

with the associated dynamics $\dot{z}=-F^{\prime}(z)$, with various examples
displayed in Fig. 6, see also attached Notebook. Initially, all cells are in
the same state (at $t\rightarrow\infty$), and then as bifurcation occurs cells
end up in two different states, clearly visible in Fig 6. The primary wave
then coincides with the bifurcation line from bistability to monostability
separating the two regions. Notice that the wavefront in the Cooke Zeeman
model is such a primary wave and that the role of the clock is mainly to
anticipate the "catastrophic jump" associated to such primary wave.

The case $\epsilon=1,\alpha=0$ gives the same example as Fig 4 A. There, the
primary and secondary wave essentially coincides because there is a very fast
relaxation of $z$ following the jump from high to low $z$ values on the
saddle-node bifurcation line. As pointed out above, this is a bit of a
particular case because the polynomial coefficients should rather mix space
and time, so that a more general case is displayed in the middle panel of Fig.
6, where $\epsilon=1,\alpha=0.02$. In such a case, the bifurcation line does
not move completely towards the posterior, so the primary wave "invades" a
portion of the field before stabilizing, leading to the sharp and fast
definition of two regions. For slow dynamics of $z$, e.g. $\epsilon=0.001$ in
the right panel of Fig. 6, the dynamics of domain separation is not sharp and
there rather is a refinement process. The primary wave is identical to the
middle panel of Fig. 6, but because of the smallness of $\epsilon$ the
dynamics take a long time to relax to smaller values of $z$, leading to the
slow propagation of a secondary differentiation wave. Noteworthy, the final
steady state in the latter case is identical to the former one but will take a
much longer time to reach, giving the feeling that some boundary sharpens,
while it was in fact defined much earlier by the primary wave.

Figure 6: Different dynamics of the primary and secondary waves described by
Eq. 3. On the left panel and middle panel, primary and secondary wave are
essentially simultaneous ($\epsilon=1$), the right panel has same $\alpha$ as
middle panel but with a much slower $\epsilon$, giving rise to an identical
(hidden) primary wave and a much later secondary wave

###  2 Meinhardt’s model

In a series of papers in the 70s, an alternative view was defended by Gierer
and Meinhardt, who proposed that reaction-diffusion processes combining
activator and inhibitors were at the origin of segment formation in metazoans
[119]. In 1977 Meinhardt applies them to fly, proposing the following model
[120, 121] :

$\displaystyle\dot{A}$ $\displaystyle=$ $\displaystyle cA^{2}/H-\mu
A+D_{a}\Delta A+\rho_{0}$ (4) $\displaystyle\dot{H}$ $\displaystyle=$
$\displaystyle cA^{2}-\nu H+D_{h}\Delta H+\rho_{1}\ $ (5)

where $\Delta=\frac{\partial^{2}}{\partial x^{2}}$ is the one-dimensional
diffusion operator. This model is ‘Turing-like’, with an activator $A$ that
self-activates and activates a repressor $H$, both diffusing. Later, in 1982,
Meinhardt argued that the addition of a segment from a growth zone, with
subcompartmentalization, required new mechanisms to produce an alternation of
Anterior and Posterior states within one segment. In particular, it is very
natural to assume there is an oscillator generating such alternation, that can
further be coupled to an external morphogen. Meinhardt calls this the
"pendulum-escapement model" :

> "Imagine a grandfather clock. The weights are at a certain level
> (corresponding to the local morphogen concentration). They bring a pendulum
> into movement, which alternates between two extreme positions. The
> escapement mechanism allows the pointer to advance one unit after each
> change from one extreme to the other. As the clock runs down, the number of
> left-right alternations of the pendulum and hence the final position of the
> pointer is a measure of the original level of the weights (level of
> morphogen concentration)."

The "extreme" positions of the pendulum correspond to the anterior-posterior
segment states, both being generated by an oscillator and modulated by the
presence of an explicit morphogen to control the pattern (e.g. the number of
segments). So while Meinhardt proposes the existence of a clock his work
differs from the Cooke and Zeeman model in a subtle but crucial way. In the
Cooke and Zeeman model, the oscillator defines blocks of cells corresponding
to somites. In Meinhardt’s model, the oscillator defines alternating fates of
genetic expression, in modern terms corresponding to somite compartments
(anterior and posterior).

To model such alternation, Meinhardt essentially combines his fly segmentation
model reproduced above with its own negative mirror image, to include another
alternating fate. Remarkably, the addition of this fate allows for the natural
emergence of oscillations. More precisely, Meinhardt assumes that two
variables are present, called $A$ and $P$ (that correspond respectively to
anterior/posterior markers of somites). $A$ and $P$ also activate fast
diffusing variables $S_{A}$ and $S_{P}$, respectively limiting extension of
$A$ and $P$, so that the pairs $(A,S_{A})$ and $(P,S_{P})$ define two (so far
independent) Turing systems. Meinhardt then adds mutual exclusion between the
two Turing systems, via a repressor $R$ which is activated similarly to both
$A$ and $P$, Fig. 7 .

Figure 7:  Meinhardt model and its reduction to the Meinhardt-VanDerPol Model.
A and P are anterior and posterior genes within the same segment, they are
mutually exclusive via the interaction with an extra variable R. SA and SP are
diffusing genes limiting the expansion of respective genes A,P. See the main
text for detailed equations.

####  1 Mathematical formulation

We could not find an explicit mathematical description of this model from
Meinhardt, but it can be reconstructed both from Meinhardt’s other similar
models and from the BASIC code used to generate his figures, found in appendix
of [86], Fig. 7, left. Meinhardt’s model can thus be described with 5
variables :

$\displaystyle\frac{dA}{dt}$ $\displaystyle=$
$\displaystyle\rho_{0}-d_{A}A+\frac{cA^{2}}{RS_{A}}$ (6)
$\displaystyle\frac{dP}{dt}$ $\displaystyle=$
$\displaystyle\rho_{0}-d_{P}P+\frac{cP^{2}}{RS_{P}}$ (7)
$\displaystyle\frac{dR}{dt}$ $\displaystyle=$
$\displaystyle\frac{cA^{2}}{S_{A}}+\frac{cP^{2}}{S_{P}}-\beta R$ (8)
$\displaystyle\frac{dS_{A}}{dt}$ $\displaystyle=$
$\displaystyle\gamma_{A}(A-S_{A})+D_{A}\Delta S_{A}$ (9)
$\displaystyle\frac{dS_{P}}{dt}$ $\displaystyle=$
$\displaystyle\gamma_{P}(A-S_{P})+D_{P}\Delta S_{P}$ (10)

Because of the presence of $R$, in the absence of diffusion, the whole system
oscillates, while in the presence of diffusion a stabilizing wavefront
propagates, converting the temporal oscillation into a spatial one [86].

The initial Meinhardt model requires 5 variables, so is rather complicated to
analyze. But we can use its natural symmetries to simplify it and extract the
core working mechanism.

To make a better sense of what happens, let us take $d_{A}=d_{P}=d$,
$\gamma_{A}=\gamma_{P}=\gamma$, and $D_{A}=D_{P}=D$. We start with a quasi-
equilibrium assumption on $R$ so that

$\beta R=\frac{cA^{2}}{S_{A}}+\frac{cP^{2}}{S_{P}}$ (12)

This gives

$\frac{d(A+P)}{dt}=2\rho_{0}-d(A+P)+\beta$ (13)

This suggests performing a new quasi-static assumption

$A+P=\frac{\beta+2\rho_{0}}{d}=C_{0}$ (14)

Notice then that $A$ and $P$ are inversely correlated, corresponding to the
intuition that they repress one another.

Similarly, we can make a quasi-static assumption for the variable
$S_{A}+S_{P}$ so that

$S_{A}+S_{P}=\frac{\beta+2\rho_{0}}{d}=C_{0}$ (15)

(basically, we make the system fully symmetrical in $A$, $P$) This allows
using symmetries in the equations to eliminate completely either $A$ or $P$.
Keeping for instance $A$, Meinhardt’s reduced model then is:

$\displaystyle\frac{dA}{dt}$ $\displaystyle=$
$\displaystyle\rho_{0}-dA+f(A,S)$ (16) $\displaystyle\frac{dS}{dt}$
$\displaystyle=$ $\displaystyle\gamma(A-S)+D\Delta S$ (17)

wifh
$f(A,S)=\beta\left(1+\frac{(C_{0}/A-1)^{2}}{C_{0}/S-1}\right)^{-1}=\beta\frac{A^{2}}{S}\frac{C_{0}-S}{C_{0}-S+(C_{0}-A)^{2}}$.
The simplification of the model is illustrated in Fig. 7 .

Notice the similarity with the initial fly model in Eqs. 4-5 : there still is
auto-activation of $A$ and repression by $S$, in particular when $A$ and $S$
are small. But the additional modulation
$\frac{C_{0}-S}{C_{0}-S+(C_{0}-A)^{2}}$ is equal to
$\frac{S_{P}}{S_{P}+P^{2}}$. This illustrates the symmetry with respect to $P$
and suggests additional non-linear effects when both $A,S$ are close to
$C_{0}$.

A simulation of this model is shown on Fig. 8 and indeed recapitulates
properties of the full Meinhardt model, see attached Notebook. Interestingly,
in the absence of diffusion, the $A/S$ dynamics is a typical relaxation
oscillator, as can be clearly seen from Fig. 8 B (see below and in the
Appendix A for general discussions on relaxation oscillators). $A$ oscillates
between two values approximately equal to $0$ and $C_{0}$, and $S$ slowly
relaxes towards $A$, Fig. 8 B. Like standard relaxation oscillators, when $S$
passes a threshold, it induces a "jump" of $A$ towards a new value
($0\rightarrow C_{0}$ or $C_{0}\rightarrow 0$) and a symmetric part of the
cycle occurs.

One also sees the effect of the $f$ function described above on the
nullclines, (i.e. lines for which respectively $\dot{A}$ and $\dot{S}$ are 0
in absence of diffusion) in Fig. 8) B right. Close to $A\sim 0$, the $A$
nullcline (blue) diverges; this corresponds to the regime where
$f(A,S)\propto\frac{A^{2}}{S}$ so that using term Eq. 16, one gets
$A^{2}/S\propto A$, i.e. $A\propto 1/S$. Close to $A\sim C_{0}$ a new regime
occurs : in this regime $f(A,S)\sim\frac{C_{0}-S}{(C_{0}-A)^{2}}$, so that
again using 16 and definitions of $C_{0}$, one gets after some algebra that
$(C_{0}-S)\propto\frac{1}{C_{0}-A}$. This regime essentially is the
"symmetrical" regime on the posterior variable $P$ of what happens for the
anterior variable $A$.

Those two regimes provide the two branches of a relaxation oscillator driving
the AP alternation. When adding diffusion on $S$, a boundary from high to low
$A$ is stable and nucleates a moving front stabilizing the pattern.

Figure 8: Simulation of the Meinhardt model. (A) Kymographs of the variables
$A$ and $S$, respectively, obtained with parameter values
$\beta=1.5,\rho_{0}=0.012,d=1,\gamma=0.01,$ and $D=0.01$. The initial
condition is an induced boundary ($S=1$ in 3 anterior cells, $S=A=0.1$
everywhere else). (B) Example trajectories of $A$ and $S$ in a cell, and flow
diagram of the $A,S$ system in a single oscillating cell, with the limit cycle
trajectory in red. Nullclines for $A,S$ are shown, blue for $A$ and green for
$S$.

####  2 Generalization : Meinhardt - Van der Pol model

A behavior similar to the Meinhardt model can be observed with many other
(symmetrical) relaxation oscillators, which are better suited for a more
precise study of what happens. This was later rediscovered by [122] and the
associated patterning mechanism was called a "progressive, oscillatory
reaction diffusion" (PORD) model (see section Somite AP patterning: Inverse
problem approach below).

Let us for instance consider the following Meinhardt-VanderPol model, based on
the addition of a diffusive term to the slow variable of a classical Van Der
Pol/Rayleigh oscillator (see Appendix) :

$\displaystyle\epsilon\frac{dA}{dt}$ $\displaystyle=$ $\displaystyle
A-S-A^{3}/3$ (18) $\displaystyle\frac{dS}{dt}$ $\displaystyle=$
$\displaystyle\lambda(A-\mu S)+D\Delta S$ (19)

Figure 9: Flow-plots of the Meinhardt-Van der Pol model for different values
of the control parameter $h$. nullclines for $A$ is in blue, nullclines for
$S$ in green, and trajectories are in red. For $h=h_{c}\sim\lambda.0.94$ the
system underges a Hopf bifurcation.

This model has only one non-linearity, the $A^{3}$ term in Eq. 18. We can
interpret this model "biologically" with A self-activating, and repressed by
$S$ (itself activated by A). Instead of the repression by $P$, this model
introduces a cubic degradation term for $A$ which makes sure that $A$ non-
linearly goes to $0$ once $|A|$ is big enough. Notice however that both $A$
and $S$ can take either positive or negative values and that the initial
symmetry of Meinhardt’s model is in fact conserved if one flips the signs of
$A,S$. Here, only $S$ diffuses (to stay consistent with Meinhardt) and we
simulate the system on a line of (discrete) cells and the pattern is stable
(consistent with the recent observation that Turing patterns are stable with
only one diffusing variable on discrete grids [123]). We also introduce a
parameter $\lambda$ allowing to modulate the dynamics of the slow variable (in
particular the period of the oscillator).

Moving now to phase space analysis, when $D=0$, the system jumps back and
forth the sigmoidal branches of the $A$ nullclines (Fig. 9 left), like a
classical Van Der Pol oscillator (see e.g. [124] Chapter 10, and Appendix). To
understand what happens when there is diffusion, let us treat the $D\Delta S$
term as an external control parameter $h$. Phase plane analysis immediately
reveals that when $h$ passes a threshold value (approximately equal to
$h_{c}\sim 0.94\lambda$), the system $(A,S)$ system undergoes a Hopf
bifurcation, due to the fact that the $S$ nullcline moves vertically and
intersects one "bistable" branch of $A$ (the negative $A$ branch in Fig. 9
middle). Notice that since the $A,S$ system is fully symmetrical, a similar
bifurcation happens when $h<-h_{c}$, with the system stabilizing on the
positive $A$ branch. So we expect that wherever the second spatial derivative
of $S$ reaches this threshold, the system stops oscillating and, depending on
the sign of this second derivative, stabilizes in a branch of either positive
or negative A. Also notice that, interestingly, the system stays excitable
even for $h>h_{c}$ (Fig. 9 right, see Appendix for the definition of
excitability).

To understand what happens when $D>0$, it is first useful to consider the
steady state situation. We see an alternation of stripes of $A,S$, where $A$
jumps from almost constant values and $S$ presents a smoother, oscillatory
profile. In particular, for $S$ we get at steady state:

$D\Delta S(x)-\lambda\mu S(x)=-\lambda A$ (20)

A crude approximation is to consider that $A$ takes almost constant positive
and negative values ($A\simeq\pm A_{0})$, then in one stripe (centered with
$0$) we expect, solving the equation, that $S\simeq\pm\left(A_{0}/\mu-
S_{x}\cosh(\sqrt{\frac{\lambda\mu}{D}}x)\right)$ at steady state. At a stripe
boundary, $A$ switches sign, so that $S$ has to be equal to $0$ by continuity
of its derivatives. This imposes that
$A_{0}/\mu=S_{x}cosh(\sqrt{\mu/D}(x_{0}/2))$, and thus defines $S_{x}$ as a
function of $x_{0},A_{0}$ which respectively correspond to the size of the
pattern and the scale of $A$ at steady state. $S_{x}$ and $x_{0}$ can not be
defined by the steady state equation in a self-consistent way, and emerge from
the dynamics. Notice that $A$ jumps while $S$ stays continuous, so as a
consequence, the control parameter has to be spatially discontinuous at steady
state.

It is then useful to plot the dynamics of the control parameter to see how
such a discontinuity appears and how the pattern forms. We show kymographs of
$A,S$ and rescaled control parameter $|D\Delta S|/\lambda$ in Fig. 11. We see
a "checkerboard" pattern of the control parameter along the front; in
particular, at well-defined, discrete times, the control parameter quickly
moves above $h_{c}\sim 0.94$ in blocks, defining stabilized regions.

The precise dynamics explaining stabilization are rather complex, as might be
expected for a system defining its control parameter through the second
derivative of a bistable variable. To our knowledge, there is no precise
mathematical study of this process. We will thus limit ourselves to a
qualitative and intuitive description of what happens, Fig. 10. Let us focus
first on the boundary of a region that has just formed. We see that posterior
to this region, the control parameter $|h(S)|>h_{c}$, so that the
discontinuity in the pattern is established and stable, Fig. 10 A. This
induces a spatial gradient of control parameter $h$: close to this
discontinuity, the region is oscillating (like the posterior) but is close to
the bifurcation point Fig. 10 B. Such dynamics of the control parameter make
sense, since after the jump, there is a discontinuity in $A$ between the
stable and the oscillating region, and we thus expect $S$ to follow in a
"smoother" way, with an increase in its second derivative.

The absolute value of the control parameter is slowly increasing in this
region in a graded way so that oscillations stabilize in more and more cells.
Eventually, the posterior oscillation (where the control parameter still is
around $0$) jumps on the other branch, Fig. 10 C, left. This creates two
domains in $A$ (one positive, one negative), between posterior cells which
have just jumped and more anterior cells where the oscillation is close/past
the bifurcation on the other branch.

Finally, because of the relaxation-oscillator dynamics, $S$ follows $A$ with
delay. This creates a sudden increase of second derivatives of $S$ at the
interface between positive and negative $A$ , and eventually a spatial
discontinuity in both $A$ and in the control parameter Fig. 10 D-E ensues.
This both nucleates the next stable region and stabilizes this region that
never jumped, and the process iterates forming a stable alternation between
regions of low and high $A$, with $S$ following $A$ in a "smoother" way.
Notice in particular that a new block stabilizes in three steps: first, a
small stable region is nucleated close to a newly formed boundary (Fig. 10 A),
then the next stable boundary is induced Fig. (Fig. 10 C) and lastly the
interior of a newly defined block between two stable boundary stabilizes (Fig.
10 D-E).

Figure 10: [

-20pt]Time evolution of the spatial pattern in the Meinhardt model. The posterior is on the left. The range of control parameter for which the system is oscillating is indicated by green lines. In (A-B), a small region left of the boundary stopped oscillating, creating a spatial gradient in S.In (C), the jump of the oscillator in the posterior-most region nucleates a new boundary that moves towards the right. The control parameter crosses the Hopf line, stabilizing the boundary around position 80, and the oscillation stabilizes in an entire block for higher positions (D). The control parameter keeps increasing left of the new boundary (E), leading to a situation symmetrical to (A).

The scaling law of this process is of particular interest. As pointed out by
[122], the speed of the front is an emerging quantity of the diffusion of the
stabilizing zone, induced in our example by the changes in the control
parameter $h$. There are a priori at least two possibilities here. First, the
speed $v$ could emerge independently from the clock, like in the initial clock
and wavefront model, so that the size of the pattern would be $S=vT$. The
other possibility could be that the speed and the clock are coupled by
diffusion so that there is a (pattern) wavelength proportional to $\sqrt{DT}$,
where $T$ is the period of the clock and $D$ is the diffusion constant. This
would then give a wavefront speed proportional to $\sqrt{D/T}$. Going back to
simulations, our numerical studies reveal that the wavelength of the pattern
is almost exactly proportional to $T$ over more than one order of magnitude of
period change (data now shown) so that the wavefront speed does not depend on
the period, similar to the former hypothesis. This is visually illustrated in
Fig. 11 : the slope of the stable region in the kymographs does not depend
much on the control parameter of the period $\lambda$. This suggests that the
front speed is purely diffusion-driven, like many other models in physics and
biophysics, see e.g [125], while the nucleation of the new stable zone is
driven by the relaxation oscillation.

Figure 11: Behaviour and scaling of the Meinhardt-VanDerPol model for
different values of parameter $\lambda$ from Eqs. 18-19. Kymographs of the
variables $A$ and $S$ are represented. In the third column, we plot
$\frac{|D\Delta S|-h_{c}(\lambda)}{\lambda}$, showing the jump in the control
parameter at the front.

####  3 Biological Interpretation of the Meinhard’s model

As pointed out by Meinhardt himself :

> "In the model I propose, the oscillation (between A and P), the wavefront
> (separating the oscillating and the stable pattern), as well as the
> spatially stable periodic pattern (of A and P), result from one and the same
> mechanism."

This simplicity in the equations explaining multiple aspects of the process
obviously has a strong appeal for physicists, especially when reduced to two
variables. As such it provides important insights into biological mechanisms,
both by setting a modeling framework and by suggesting predictions.

First, the pattern in Meinhardt’s model is clearly stabilized by interactions
of consecutive domains where $A$ is present/absent. So spatial diffusion is
crucial to form and stabilize the boundary. This somehow contradicts the
kinematic view of somites formation associated to the robustness to various
embryonic manipulations (graft, spatial boundaries), with the caveat that
those manipulations are at the tissue scale and might not be the best to
falsify local mechanisms at the cellular scale.

Second, as explained above, there is no discrete formation of a block of cells
defining somites like in the Cooke and Zeeman model. Rather, somites are
assumed to be defined a posteriori as the concatenation of one anterior
compartment with a posterior one. Then, the alternation of a APAPAP pattern
does not define unambiguously a somite, since boundaries should be defined for
the P to A boundary but not for the A to P boundary. Meinhardt, therefore,
suggests that there might be a third oscillating variable (called X) so that
the real alternation is of the form APXAPX, unambiguously defining the somite
boundary. In fact, Meinhardt points out another potential mechanism, where the
system might rather detect the temporal succession of P to A in opposition to
A to P to trigger boundary formation :

> "Imagine a ship in a channel system with locks. A lock can be in two states.
> Either the lower gate is open and the upper gate is closed or vice versa. In
> neither state can a ship pass through. But in one state the ship can enter
> into the lock and after the switch to the other state, the ship can pass. In
> one state, the transition is prepared but blocked. In the other state, the
> block is released, the transition can take place, but no preparation of the
> next transition is possible. For the sequential activation of control genes
> I assume that, for instance, in the P-state a substance X is produced that
> activates the subsequent gene, but that its action is blocked. In the
> A-state, the block is released but X is no longer produced. Only with a P-A
> transition the activation of the subsequent gene can take place due to the
> simultaneous release of the block and the presence of the substance X. In
> contrast, activation of subsequent control genes can not occur if cells
> remain permanently in the P- or in the A-state."

In modern terms, this describes by essence a phase detector downstream of an
incoherent feedforward loop network (see e.g. [126]), where $P$ activates $X$
but represses its downstream target, while $A$ derepresses the target. $X$ is
produced only when $P$ is fading out and $A$ increasing.

Like the Clock and Wavefront model, the differentiation wavefront is emerging
from the dynamics. One can first approximate the wavefront in Meinhardt’s
model as the point where the oscillation stops (or the limit cycle
disappears), but as seen from simulations, this is not a continuous front,
rather, due to the relaxation oscillation, it jumps in a discontinuous way
from one boundary to the other, later-on stopping oscillations in-between.
This jumping process is not so different qualitatively from the pulses of the
clock inducing transitions in the Cooke and Zeeman model. The dynamics of
motion are very different though: in the Cooke and Zeeman model, the
competency zone for transition to bistability is defined by the external
positional information variable $p$, while in the Meinhard’t model, it rather
is a self-organized ’domino’ effect where one stable region nucleates the
following one with the help of the ongoing relaxation oscillator and
diffusion. This creates a difficulty for scaling/changing the size of the
pattern. In particular, in Meinahrdt’s model there is no external positional
information variable independent from the oscillation. Meinhardt anticipates
this potential difficulty by introducing a modulation to his model, adding a
spatial dependency in equation $A$ of the form :

$\frac{dA}{dt}=\rho_{0}-d_{A}A+\frac{cA^{2}}{R(S_{A}+\Theta(x,t))}$ (21)

This threshold $\Theta(x,t)$ de facto defines some external positional
information in the system, which can modulate the speed of the clock and as
such the size of the pattern. Meinhardt suggests a simple model so that
$\Theta$ is essentially monitoring the number of cycles in relation to a
morphogen gradient. By adjusting the slope of the morphogen gradient in a
size-dependent way, scaling with embryonic size can be obtained. This model
can be further adapted to account for further specialization of some segments
as a function of time (e.g. in the case of insects, some segments will give
rise to wings while other ones will give rise to halteres, due to the
expression of so-called Hox genes).

###  3 Cell Cycle model

In the early 90s, Stern and co-workers proposed that the segmentation clock
could be in fact related to the cell cycle [127]. This comes from a series of
clever experiments in chick showing very striking features [128, 129] :

  * •

One single heat shock produces several segmental anomalies, restricted to one
or two consecutive segments, but separated by 6 to 7 somites - corresponding
to roughly 9 hours of development. This suggests the existence of a long
temporal cycle implicated in segment formation, with a length corresponding to
the time required to form 6-7 somites. Then if this cycle is initially
perturbed, the perturbation would be repeated every 6 to 7 somites,
corresponding to the period of the oscillator.

  * •

The 9 hours period was later shown to correspond to the length of the cell
cycle, strongly suggesting that it is coupled to somite formation.

  * •

A single progenitor cell in the tail bud injected with dye gives rise to
several clusters of cells in the PSM and in somites, with a 6 to 7 somite
periodicity [130, 127]

This suggests the following picture: progenitors in the tail bud constantly
divide and lay down cells in the PSM in an ordered way so that cells at the
same anteroposterior position are roughly at the same phase of their cell
cycle. A 6-7 somite periodicity thus recapitulates spatially a phae gradient
of the cell cycle. Then, the cell cycle is coupled to somite formation, for
instance, there might be a special phase $\phi_{*}$ of the cell cycle for
which cells form a boundary when they reach the anterior. Now we need to
assume that one cell cycle phase (say $\phi_{S}$) is specifically sensitive to
heat shock (while other phases of the cycle would not be), which could well
happen for discrete events in the cell cycle (e.g. a transition between G1 and
S/G2/M). So when heat shock occurs, it disrupts all cells in $\phi_{S}$, not
only the older cells in the anterior but also the cells just laid in the
posterior a few cell cycles later. When those perturbed cells end up
differentiating into somites, theoretically at phase $\phi_{*}$, their
disrupted cell cycle results in segment anomalies. The cells just posterior to
this anomaly were not in phase $\phi_{S}$ at the time of the heat shock, so
are laid down normally and form somites at $\phi_{*}$. Then, one full cell
cycle later, cells that were again at $\phi_{S}$ at the time of the heat shock
would theoretically reach $\phi_{*}$ but are disrupted again. This explains
why one single heat shock disrupts several segments in a periodic way.

In [131], McInerney et al.proposed a mathematical implementation of the cell
cycle model for somitogenesis. Essentially, the goal is to understand with a
realistic biochemical model how a spatial gradient of cell cycle phases can
translate into blocks of simultaneously differentiating somites. In
particular, this model is not concerned with the formation of stripes or AP
somite polarity (contrary to Meinhardt’s model). From a modelling standpoint,
the challenge is to find how a continuous periodic process (such as the cell
cycle, with a spatial gradient of phases) can give rise to a discrete output
(spatially extended somite blocks), and as such, while details differ, this
model is in fact very close to the initial Clock and Wavefront vision. This
model is also of particular interest from a conceptual standpoint because many
subsequent models implement similar ideas with different hypotheses on the
nature of the oscillator or of the front.

The model relies on the combination of two continuously moving fronts with a
simple, two-component biochemical network, encoded into the following
equations :

$\displaystyle\frac{\partial u}{\partial t}$ $\displaystyle=$ $\displaystyle
f(u,v)$ (22) $\displaystyle\frac{\partial v}{\partial t}$ $\displaystyle=$
$\displaystyle g(u,v)+D\frac{\partial^{2}v}{\partial x^{2}}$ (23)

The $f$ and $g$ functions encode generic signaling dynamics where $u$ self-
activates, and is activated by $v$, while $v$ is repressed by $u$. After
dimensionless reduction, one gets :

$f(u,v)=\frac{(u+\mu v)^{2}}{\gamma+\kappa u^{2}}\chi_{u}-\frac{u}{\kappa}$
(25)

and

$g(u,v)=\frac{1}{\epsilon+u}\chi_{v}-v$ (26)

Two fronts moving with speed $c$ are encoded into a spatial-temporal
dependency of the activations on $u,v$ :

$\chi_{u}=H(ct-x+x_{1})\qquad\chi_{v}=H(ct-x+x_{2})$ (27)

where $H$ is the Heaviside function. A cell cycle gradient is imposed by the
fact that $x_{2}<x_{1}$ : so cells become competent to express $u$ before they
are competent to express $v$. Practically, the couple $(\chi_{u},\chi_{v})$
can only take three values $(0,0),(1,0)$ and $(1,1)$. Those three values
correspond to three spatially distinct regions of the embryo, respectively
corresponding to the posterior of the embryo (region $I$), a somite definition
zone (region $II$), and the anterior of the embryo (region $III$).

Figure 12: (A-C) Simulation of the cell cycle model for somitogenesis with
parameter values
$\mu=0.0001,\gamma=0.01,\kappa=10,\epsilon=0.001,D=60,$and$c=0.00125$. (A)
Kymograph of the variable $u$, with blocks of cells moving from low ($u=0$) to
high ($u=1$) state. (B) Kymograph of the variable $v$ showing spatially
extended transient pulses. (C) Propagation of the fronts, shown as the sum of
activations $\chi_{u}+\chi_{v}$. Three distinct regions are indicated: the
posterior ($I$), the somite definition zone ($II$), and the anterior ($III$).
(D-G) Keeping the rest of the parameters the same as in (A), we change one
parameter value in the simulation. (D)$D=30$. (A). (E) $D=120$. (F)
$\kappa=20$. (G) $\mu=0.0003$.

It is useful to first study the behavior of the $u,v$ system for constant
values of $(\chi_{u},\chi_{v})$ corresponding to different regions. The
posterior of the embryo (region $I$) is the simplest case: $\chi$s functions
are $0$ so that the only steady state is $(u,v)=(0,0)$.

In region $II$, the steady state value of $v$ still is $0$, but self-
activation of $u$ creates a new stable steady state
$u_{*}=\frac{1}{2}(1+\sqrt{1-4\gamma/\kappa})$. We will also assume for
simplicity that $\gamma<<\kappa$ so that $u_{*}\sim 1$. We thus expect region
$II$ to be a region of bistability where $u$ can go to stable values $0$ or
$\sim 1$, depending on the initial conditions of both $u$ and $v$.

In the anterior, region $III$, $v$ can no longer be strictly $0$. Assuming
initially $u\sim 0$, then we see that $v$ initially rises quickly to a (high)
value $\sim 1/\epsilon$. Then since $v$ activates $u$, the value of $u$ starts
increasing as well. The final state depends on parameters, but if $\epsilon$
is small enough, $v$ is transiently so high that it strongly activates $u$,
which can become high enough to sustain its own production. In turn, $u$ then
squashes down $v$, to get a steady state not far from the steady state value
$u_{*}\sim 1$, corresponding to a differentiated, somite state. McInerney et
al.proposed that the system is in fact monostable for those parameters,
leading to a high, sustained steady state approximately equal to
$\sim(u_{*},1/(\epsilon+u_{*})$. Notice these dynamics also create a transient
"pulse" of $v$, before going back to a close to $0$ value of $v$, so akin to
an excitable system (see Appendix).

So we see that in region $I$ we expect $u$ to be $0$, in region $II$, $u$ can
be essentially $0$ or $1$, while in region $III$, $u$ is essentially $1$.

Now the idea underlying the full model is that the Heaviside function combined
to diffusion will induce a sudden transition of $u$ from $0$ to $1$ in a block
of cells in the region $II$, via a spatially extended pulse of $v$. Let us
assume that at $t=0$, $v$ is roughly $0$ everywhere, $u\sim 1$ for $x<0$ and
$u=0$ otherwise. Then as time increases, $\chi_{v}$ is turned on close $x\sim
0$, so that there is a sudden pulse of $v$ there. If the diffusion constant is
very high, this $v$ pulse is going to diffuse very quickly towards higher $x$,
leading the pulse to be spatially extended. What happens then depends on the
balance of the parameters, but after some time (say $\tilde{t}$), this pulse
of $v$ induces a transition from $u=0$ to $u=1$ values in a region close to
$x\sim 0$. The size of this region depends on parameters. If the diffusion is
fast enough, induction occurs in the entire region where $\chi_{u}=1$, but if
this region is big enough (or diffusion too small), this will happen only in
part of the region where $\chi_{u}=1$ (see Fig. 12).

When $u$ is high enough to self-sustain, it pushes $v$ back towards $0$ in
this region. So we end up with an entire region where, after the pulse of $v$,
all cells "commit" simultaneously to a high $u$ state, corresponding to
discrete somite formation. Once this transition has occurred, $v$ is again $0$
everywhere, $u\sim 1$ for a more extended region, and $u=0$ otherwise, so that
this process can start again.

In the initial model, it was assumed that the transition happens in the entire
region where $\chi_{u}=1$ because of fast diffusion. What sets the size of the
block then is the time $\tilde{t}$ for $v$ to activate $u$ everywhere in this
region, and the size of the block of the activated bock of cells then is
$c\tilde{t}$. But if diffusion is not fast enough or the region where
$\chi_{u}=1$ is too big, the $v$ pulse will propagate from $x=0$ and activate
cells in a more localized region. The size of the pattern thus is a complex
function of all parameters, including diffusion (see different examples in Fig
12 D-G, and attached Notebook).

In summary, this model allows for the formation of somites by the generation
of periodic pulses close to the anterior PSM boundary, synchronously expressed
in a field of cells, triggering commitment to somite fate (modeled via a
bistable variable $u$). Notice that the dynamics of $u$ thus is very similar
to the variable $z$ in the Cooke and Zeeman model, Eq. 1. $v$ plays the same
role at the pulsatile clock in Eq. 2, interpreted as the cell cycle. It is
also worth comparing how primary waves (in the Zeeman sense [114]) are encoded
in both models: in the initial Clock and Wavefront model, the $(t,p$)
potential associated to the cusp catastrophe was creating an emerging
transition from a bistable to a monostable region, while here, a similar
primary wave is created by the region II to region III transition, Eq. 27,
when $v$ is activated and ensures that $u$ is monostable. In other words, the
primary wave is defined by $\chi_{v}$.

The big difference comes from the dynamics of variable $v$. First, similar to
Meinhardt’s model, diffusion of $v$ is crucial to define the pattern
(switching $u$). $u$ also shuts down $v$. This ensures coordination between
the state variable and the clock, a possibility we alluded to at the end of
the description of the clock and wavefront model. It is also noteworthy that
the oscillator is in fact not explicitly modeled in this $u,v$ model, and
rather emerges as a consequence of the sliding window $\chi_{v}$ which creates
a pulsatile window of expression of $v$ in region II. So there is no explicit
need for, say, a posterior oscillation (in the region I) like in Meinhardt’s
model. It is, in particular, not entirely clear how the differential sliding
window would practically connect to the phase of the cell cycle oscillator,
and how the initial proposal that heat shocks disrupt specific phases in the
cell cycle would be accounted for in this model.

##  Chapter 3 Phase models

On the one hand, the vast number of molecular players implicated in
somitogenesis is daunting from a theoretical standpoint, since it is not clear
how and what to model in a predictive way. On the other hand, the
phenomenology of the segmentation behavior still is relatively simple, with
waves of genetic expression sweeping from posterior to anterior, leading to
patterning. This suggests first following the spirit of classical models
described in Section Early models to focus on rather phenomenological models,
not specifically tied to actual genes. Similar issues arise for oscillators in
neuroscience and physiology and motivated the development of a "phase-based"
approach to describe more explicitly the segmentation clock dynamics, which we
briefly summarize here (see also Appendix, treatments of various complexities
can be found in [132, 133, 134, 135]). In line with our previous observation
that the clock seems to be tightly connected to Zeeman’s primary wave of
differentiation, one challenge is to tie those phase descriptions to both
clock stopping and patterning.

###  1 From chemical equations to phase

Consider a (biological) oscillator described by equations in the space of its
components, e.g. mRNA/protein concentrations :

$\frac{d\mathbf{X}}{dt}=F(\mathbf{X})$ (1)

Given some initial conditions, the system relaxes to the limit cycle, which is
a closed curve in the space of concentrations. The position on this curve can
thus be indexed by a single parameter. We define the phase $\phi(t)$ of an
oscillator by :

$\frac{d\phi}{dt}=1$ (2)

and express the phase $\phi$ modulo $T$, where $T$ is the period of the
oscillator, to account for the fact that the system is periodic (notice there
are other conventions, i.e. one can rescale time so that the period is either
$1$ or $2\pi$). In this formalism, the phase of an oscillator is nothing more
than a (rescaled) time variable on the limit cycle. For instance, if we
rescale time so that the period of the oscillator is $2\pi$, phase $\pi/2$
means that the oscillator is at the $1/4$ of its cycle, phase $\pi$ means that
the oscillator is at half its cycle, and phase $2\pi=0\quad\mathrm{mod}\quad
2\pi$ is the initial phase corresponding to the full period. Notice that
phases also correspond to positions in the space of protein concentrations,
i.e. $\phi(t)=\phi(\mathbf{X}(t))$ where $\mathbf{X}(t)$ is the value of
protein concentrations at time $t$ on the limit cycle.

There are now two important observations from the modeling standpoint :

  1. 1.

It is possible to extend the definition of phase for points outside of the
limit cycle. Imagine for instance that at a given time $t_{p}$ you first
perturb the system, e.g. by making a change
$\mathbf{X}(t)\rightarrow\mathbf{X}(t)+\Delta\mathbf{X}$, then let the
oscillator relax. Eventually, the system will go back to the limit cycle,
where you have defined a phase using Eq. 2. But then, since the phase is
nothing more than time, from this phase on the limit cycle, you can go back in
time on the trajectory you have just followed to define a phase corresponding
to the initial condition $\mathbf{X}(t)+\Delta\mathbf{X}$ at time $t_{p}$.
This way, you can define a phase for all vectors $\mathbf{X}$, even outside
the limit cycle, defining so-called "isochrons", or lines with identical
phases.

  2. 2.

for any limit cycle oscillator, the amplitude is stable (so not easily changed
by a perturbation) while the phase is neither stable nor unstable [132]. Thus,
weak perturbations of an oscillator only change its phase.

Those two properties essentially mean that, for many purposes, the behavior of
a (perturbed) limit cycle oscillator can be entirely captured by its phase
behavior, which remarkably allows us to go from complex dynamics in a high
dimensional system to only one phase variable for a given oscillator. For
instance, imagine two coupled oscillators, then if their coupling is
relatively weak, the perturbations induced by each oscillator onto one another
will stay close to the limit cycle in the initial mRNA/protein space, and one
can use the isochron theory to translate any coupling into effective phase
equations. While it is clear that this substitution is not trivial, and
computations of phase responses can be quite tricky (and has to be done
numerically for a more complex system, see Appendix A for the Adjoint method
and Malkin theorem), some generic simplifications also arise from the
periodicity of the coupling and symmetry in the equations [133] (see
Appendix). One can then use such formalism to study all kinds of effects, from
entrainment to changes of the intrinsic period. In summary, if the limit cycle