is not too perturbed and the coupling not too strong, all properties of the
oscillators under various hypotheses can then be defined in terms of the
phase, which allows for more powerful treatment, here in the context of
segmentation clock/waves.

###  2 Clock and Unclock

It is a good place to briefly mention nuances on the ‘clock’ notion for
biology, introduced by Winfree in the context of the circadian clocks [136]. A
‘simple clock’ is characterized by the fact that

> "its possible states can be arranged in a recurring sequence, in which each
> state induces the next in a repeating cycle. The "state space" is a circle
> and there is no state which is not represented in that space. The clocks of
> home and industry are of this kind; the only variable quantity in the clock
> is the angular position of its meshed gears, and to get to a new position it
> must pass forward or backward through all intermediate positions."

An example with discrete states is shown in Fig. 1 A. Phase descriptions
introduced above correspond to a continuous limit of such clock definition.
Winfree notices that neuronal oscillators represent potential biological
examples of such clocks because they can by and large be modeled by one
variable that ‘resets’ at the end of its cycle (see Appendix A).

Winfree contrasts this with multiple ‘unclocklike’ behaviors. The first
example he gives is limit cycles (possibly in a high dimensional space) Fig. 1
B. The clearest difference is in the nature of states: in a clock, states are
phases, while in a limit cycle, even extended with isochron theory [134],
there are ‘phase-less’ states, such as the center of the cycle (intersection
of green and blue nullclines in Fig. 1 B). One way to explicitly show the
‘unclock’ behavior is to induce a (strong) perturbation of the limit cycle to
reach this point. However, finding the ‘right’ perturbations towards the
phase-less point can be challenging. An indirect way to prove the existence of
such point rather is topological, to get so-called Type 0 phase resetting in
response to strong perturbations, which manifests itself by a sudden 2 $\pi$
jump in Phase Resetting Curves (see Appendix for definition) [136, 134, 137].
It sounds a posteriori surprising that limit cycles are not considered as
clocks, but as we will see in this review, there are indeed very relevant
differences between phase models and explicit ODE-based models in the
segmentation context.

Another possible ‘unclocklike‘ behavior in response to the perturbation is a
long-term memory effect in response to a perturbation. The harmonic oscillator
is not a clock for this reason: it is defined by its amplitude and its phase,
and as energy is injected, its amplitude increases and does not return to its
initial value. Similarly, any transient stimulus that would impact long-term
properties of an oscillator (e.g. the period) excludes the simple clock, see
e.g. Fig. 1 C. Surprisingly, such effects have been observed in multiple
biological systems, see e.g. fly circadian clocks [136], or for ‘period’
response in human clapping experiments [138].

In the context of somitogenesis, one should also distinguish the individual,
cellular oscillators, from the global periodic one. When talking about "the
segmentation clock", it is important to specify which level is actually
studied, because the ‘clock’ or ‘unclock’ like properties of the oscillatory
behavior at the embryonic level might not be the same as the properties at the
single cell level (even though they are of course related).

Figure 1:  Clock and Unclock (A) The standard Clock model: a system where
states (discrete or continuous) are visited periodically, in a circular way.
States are represented by orange dots. (B) Unclock behavior: a system with a
limit cycle (orange trajectory) is not a clock, in particular, because there
is a phase-less point at the center (intersection of green and blue
nullclines). (C) Unclock behavior: even the simplest harmonic oscillator is
not a clock because there are multiple possible cycles (orange lines). In
particular, the amplitude of the oscillation can be changed by a transient
stimulus.

###  3 The Lewis Phase Model (LPM)

We now have the theoretical framework to discuss phase-based models, aiming at
phenomenologically describing the dynamics of somitogenesis. The first
historically important model was in fact introduced in the appendix of the
Palmeirim et al.paper describing the segmentation clock for the first time
[5]. There, Julian Lewis briefly described a simple mathematical model
recapitulating the observations of c-hairy behavior. This model is important
for at least three reasons :

  * •

it is the first example of a "phase-model" for segmentation clock, tying
clock, pattern formation, and "wavefront"

  * •

Because of this, many subsequent models of the segmentation clock can be
related to and contrasted with this initial model, as will clearly appear in
the remainder of this tutorial

  * •

similarly, many standard observations and calculations (and Ansatz) for the
segmentation clock can be illustrated first on this model

For this reason, this model deserves a complete section and analytical study
that we perform here. In the remainder of the text, this Julian Lewis Phase
Model will be subsequently abbreviated in LPM.

The LPM assumes that each cell is a phase oscillator, with a local phase
$\theta(x,t)$ ($x$ being the position of the cell, increasing $x$
corresponding to the posterior to anterior direction, and $t$ the time). The
phase also depends on a variable called maturity $m$, taken to be, in rescaled
units, $m=x+t$. In particular, notice that cells of given, constant, maturity
$m=m_{0}$ define a moving front $x=m_{0}-t$ of speed $-1$. Then the assumption
is that the instantaneous angular velocity/frequency of the clock is a simple
function of $m$, namely that

$\dot{\theta}=r(m(x,t))$ (3)

where $r(m)$ is a smooth function so that $r(m)\rightarrow 1$ for $m<<0$
(corresponding to the posterior of the embryo) and $r(m)=0$ for $m>>0$
(corresponding to the anterior of the embryo). A kymograph of this model is
shown in Fig. 2A. In this model, each cellular oscillator is slowing down with
time, but since $m=x+t$, this slowing down is time and space-dependent, giving
rise to a pattern, with a differentiation front moving from top to bottom.

Figure 2: Simulation and calculations for the Lewis phase model. Antero-
posterior axis $x$ increases from bottom to top (A) Kymograph of
$\cos(\theta)$ from a numerical simulation with $\sigma=10$. A spatial region
from $-L$ to $0$ (initial front position) is shown. (B) Analytically
calculated spatial phase profile $\theta(x)$ at different time points. (C)
Analytically calculated phase gradient at different time points and comparison
to the approximation for $x<<0$ described in the text. (D) Function
$\psi(m)=\phi(m,t)-t$, describing the phase profile for long enough times $t$.
For the numerical solution, the phase in the moving frame of reference
$\phi(m,t)$ was taken from the diagonal of the kymograph. (E) Divergence of
the wavelength $\lambda(x)$. Analytic and numerical solutions calculated from
the phase gradient; asymptotic solution calculated with the asymptotic
$\psi(m)$.

To make this intuition more quantitative, it is helpful to see what happens
for simple forms of $r$ where all calculations can be done analytically. In
the initial Palmeirim paper, no calculation is performed and the model is only
simulated, but it turns out it is completely analytic so it is worth solving
here explicitly to explain in detail what happens.

Lewis first assumes that

$r(m(x,t))=\frac{1-\tanh(m)}{2}=\frac{1}{1+\exp{2(x+t)}}$ (4)

one gets by direct integration (assuming $\theta(x,0)=0$ for all $x$)

$\theta(x,t)=\frac{1}{2}\left(t+\ln\frac{e^{x}+e^{-x}}{e^{x+t}+e^{-(x+t)}}\right)$
(5)

We then have the following interesting limits :

  * •

for $x\rightarrow-\infty$,
$\theta(x,t)\sim\frac{1}{2}\left(t+\log{e^{-x+(x+t)}}\right)\sim t$,
indicating that the clock is oscillating with angular velocity $1$ in the
"posterior" region. Conversely, for $x\rightarrow+\infty$ we get
$\theta(x,t)\sim\frac{1}{2}\left(t+\log{e^{x-(x+t)}}\right)\sim 0$, so that
the phase is stationary in the "anterior" region.

  * •

for $t\rightarrow+\infty$ and $x$ negative enough ,
$\theta(x,t)\sim\frac{1}{2}\left(t+\log{e^{-x-(x+t)}}\right)\sim|x|$,
indicating a stationary phase linear in $|x|$ for long times (Fig. 2B). This
indicates that after a long enough time, a periodic stationary pattern of size
$L=2\pi$ is reached.

A very interesting feature of this model is noticed in [5]:

> "Note that the model correctly mimics the curious way in which successive
> waves of c-hairy1 expression, initially broader than one somite, appear to
> move rapidly forward in the presomitic mesodern, squeezing up together as
> they go."

Such behavior is not predicted by any of the classical models in Early models
and thus is accounted for and explained for the first time by the LPM.

The local wavelength at position $x$ is given by the inverse of the phase
gradient i.e.

$2\pi\lambda(x,t)=\left|\frac{\partial\theta(x,t)}{\partial
x}\right|^{-1}=2\left|\tanh(x+t)-\tanh(x)\right|^{-1}$ (6)

Again, from this equation it is clear that for $x$ negative enough and
$t\rightarrow+\infty$, we have
$2\pi\lambda(x,t)=2|\tanh(x+t)-\tanh(x)|^{-1}\simeq 2/(1-(-1))=1$, so that at
stationarity the pattern is periodic with size $2\pi$. Conversely, fixing $t$
but in the limit of $x\rightarrow-\infty$ the wavelength clearly diverges with
both $\tanh$ going to $-1$, indicating the (infinitely) broad wave pattern
observed corresponding to all oscillators oscillating in phase in the
posterior-most region. The transitory region between those two domains
corresponds to the behavior described in the quote above. To further study
this transition, in a region with $x$ very negative (i.e. a region initially
in the oscillating phase), we have $\tanh(x)=-1$ which gives :

$2\pi\lambda(x,t)=2\left|\tanh(x+t)+1\right|^{-1}=1+e^{-2(x+t)}$ (7)

We clearly see from this equation how the wavelength at position $x$ goes from
$+\infty$ to $1/2\pi$ as a function of time, and conversely, we see that for a
given time, along the $x$ axis, the wavelength diverges exponentially in space
with a characteristic length of $1/2$. This length corresponds to the size of
the transitory region mentioned above, between the oscillating region and the
static region. This characteristic length can be further adjusted by
modulating the function $r(m)$, e.g we can generalize by taking

$r(m)=\frac{1-\tanh(m/\sigma)}{2},$ (8)

which gives

$2\pi\lambda(x,t)=2\left|\tanh\frac{x+t}{\sigma}+1\right|^{-1}=1+e^{-2(x+t)/\sigma}$
(9)

so that as $\sigma$ increases (corresponding to a shallower gradient), the
region of divergence of wavelength increases.

In fact, we get in both cases that

$2\pi\lambda(x,t)=(1-r(m))^{-1}$ (10)

so that (the inverse of) the wavelength is a direct measurement of the
frequency gradient. This is a very general result that can be derived directly
by making a change of variable in equation 3 to express the phase in the
moving frame of reference $\phi(m,t)$ to get

$\frac{\partial\phi}{\partial t}+\frac{\partial\phi}{\partial m}=r(m)$ (11)

This change of variable simply adds an advection term
$\frac{\partial\phi}{\partial m}$. Since we expect that $\phi(t,m=-\infty)$ is
proportional to $t$ by construction (to ensure spatially uniform oscillations
in this region), after some transitory time, we can hope to find a solution
satisfying the following Ansatz :

$\phi(m,t)=t+\psi(m)$ (12)

which nicely separates $t$ and $m$. Such Ansatz is not the full solution of
the differential equation with the imposed initial condition: rather, it is an
asymptotic behavior that we expect to reach after a transitory time, that will
be observed experimentally if we wait long enough.

This Ansatz allows us to explicitly solve the full phase profile for any form
of $r(m)$. Coming back to Eq. 12, we then have $\partial\phi/\partial t=1$ and

$\displaystyle\frac{d\psi}{dm}$ $\displaystyle=$ $\displaystyle r(m)-1$ (13)

which can be explicitly integrated to get the value of $\psi(m)$, entirely
solving the problem. $\psi$ represents the stationary phase profile in the
moving frame of reference. In particular $\psi$ captures all information on
the spatial dependency on the system. So any snapshot property measured
experimentally at a given time (e.g. phase gradient in the PSM) are related to
$\psi$. For instance, Eq. 10 just expresses the fact that wavelength is given
by

$2\pi\lambda(x,t)=\left|\frac{\partial\psi}{\partial m}^{-1}\right|$ (14)

Second, $r(m)-1$ goes from $0$ to $-1$ as $m$ varies from $-\infty$ to
$+\infty$, so that, as we integrate it, $\psi(m)\sim 0$ when $m=-\infty$ to
$\psi(m)\sim-m$ when $m=+\infty$. $\psi(m)=0$ is the region where $\phi=t$,
i.e. the oscillating region. $\psi(m)=-m$ gives $\phi(x,t)=t-m=-x$: this is
the region where the oscillator freezes to give the pattern. Notice those
behaviors are very general and do not assume anything on the shape of $r$
except it asymptotic limits.

Third, since $\psi$ depends on $m$, this means that the stationary phase
profile defined by $\psi$ simply moves with constant speed in the original
frame of reference defined by $x$. So $\psi$ captures the transition zone,
moving from the region with stationary phase toward the oscillatory region.
Fig. 2D represents the shape of $\psi$ for $r(m)$ given by Eq. 8.

Specific forms of $r$ only change the phase profile in the transition zone, as
long as the asymptotic values $0$ for $m\rightarrow+\infty$ and $1$ for
$m\rightarrow-\infty$ are fixed. So waves are observed, emerging from the
oscillating zone, and traveling from posterior to anterior before stabilizing
(notice that those waves are purely kinematic, they are not corresponding to
the propagation of any signal, rather they purely come from the change of
frequencies induced by $r$). Waves are locally moving with speed
$v=\frac{\partial\phi}{\partial t}/\frac{\partial\phi}{\partial
x}=\exp{-2(x+t)/\sigma}$ so that $v=\frac{r}{1-r}$. In particular, the wave
has infinite speed in the posterior (corresponding to the infinite wavelength)
and stabilizes as they move toward the anterior.

It is worth computing how many waves are traveling, which is a simple
biological observable (we reproduce with adaptations here a calculation done
in [139], for which the wave pattern exactly matches the Lewis Phase Model as
will be discussed in section Wavefront as a phase shock emerging from
coupling).

One can easily compute the total phase gradient between the oscillating region
and the static region, and the number of traveling waves simply is the
gradient divided by $2\pi$. We express the total phase difference between the
tail bud and position $m_{0}$ :

$\Delta\psi_{0}=-\int_{-\infty}^{m_{0}}\frac{\partial\psi}{\partial
m}\mathrm{d}m=\int_{-\infty}^{m_{0}}\frac{\mathrm{d}m}{1+\exp{-\frac{2m}{\sigma}}}$
(15)

The problem is that this integral diverges for $m_{0}\rightarrow\infty$.

But of course, no embryo is infinite, even though it helps simplifying
calculations. So the realistic way to perform this calculation is to put a
cut-off on positive $m$, since the PSM is of finite size $L$. A reasonable
assumption is to assume that $m=0$ (where $r=1/2$) corresponds to mid-PSM,
meaning that the "real" PSM starts at $m=-L/2$ and ends at $m=L/2$. One then
gets :

$\Delta\psi_{0}=\int_{-L/2}^{L/2}\frac{\mathrm{d}m}{1+\exp{-2m/\sigma}}=\frac{\sigma}{2}\ln\frac{e^{L/\sigma}+1}{e^{-L/\sigma}+1}\simeq
L/2$ (16)

assuming $L/\sigma$ large enough (meaning that the PSM length is bigger than
the typical length scale for frequency change). The number of waves is
$\Delta\psi_{0}/2\pi$. Notice that $2\pi$ is the length of the pattern in
rescaled units, which means that going back to actual physical units, the
number of propagating waves is :

$N\simeq\frac{L}{2S}$ (17)

So from this model the number of waves simply is the PSM length expressed in
units of somite size, divided by 2. This is a remarkably simple prediction
that works well for many organisms. For instance, as pointed out by [139]
using data from [108], for zebrafish, the PSM length is $48$ cellular
diameters, a somite length is $7$ to $8$ diameters so that the number of
traveling waves is around $3$, which visually fits (dynamical) data. For snake
[108], $\frac{L}{2S}\simeq 12$, which approximately fits the number of
traveling waves observed.

Intuitively this scaling proportional to $L$ comes from the fact that if
$L>>\sigma$, the integrated quantity $\frac{\mathrm{d}m}{1+\exp{-2m/\sigma}}$
is essentially equal to $1$ for $m$ higher than a few $\sigma$, and so the
phase difference with respect to the oscillator basically accumulates almost
linearly in time (and space) when the frequency of cells is very small,
between positions $m>\sigma$ and $m=\frac{L}{2}$. In other words, the $L/2S$
scaling merely reflects the fact that cells get rather slow over a significant
portion of the PSM (here half of it), so that, maybe not surprisingly, the
phase difference with respect to the oscillatory region essentially is
proportional to the length of the region of significantly slower oscillation.
Another way to frame this result is that, if the clock frequency was suddenly
dropping as cells enter the PSM, we would expect to have at most $L/S$ waves,
because this simply is the number of somites that should come out of a PSM of
length $L$. In reality, the slowing down takes some time to be significant, of
the order of half the PSM, so the effective phase shift will be of the order
of $L/2S$. So we expect this result to be in fact quite generic; as an
illustration, a different calculation is also performed in Supplementary Box 2
of [108], accounting for the exponential growth of the PSM but in the end
leads to the exact same result (see Appendix Section Number of waves in a
growth model).

Those calculations show that despite its simplicity, the LPM includes
nontrivial and observable features not present in the initial clock and
wavefront model: both a spatial and temporal dependencies of the clock, which,
as cleverly noticed by Palmeirim et al., explain in a simple way the observed
wave pattern. In fact, one could almost say this feature is contradictory with
the initial clock and wavefront model, which postulates that the clock and the
wavefront are two independent variables. In the LPM model, for finite
$\sigma$, there clearly is nothing such as a discrete wavefront, rather the
clock continuously stops, giving rise to traveling waves in the transitory
regime. Also, the oscillators eventually stop and form a spatial pattern. Both
those aspects are more reminiscent of Meinhardt’s model, and indeed suggest
that the segmentation process should be considered in its continuity from
oscillation to stabilization to be fully understood. In other words, it is
important to consider segmentation waves and not only segmentation
oscillators.

Two differences with Meinhardt’s model are nevertheless worth pointing out:

  * •

The LPM is purely cell autonomous, in the sense that the behavior of
individual cells is entirely prescribed by the dynamics of variable $m$, which
is externally controlled. In Meinhardt’s model, the slowing down of the clock
depends on the interactions of different oscillators and as such is an
integral part of the model. One can not exclude that the maturity $m$ itself
is coming from the interactions between cells or within cells, which are not
described in the LPM model. We will see an explicit example of this in Section
Wavefront as a phase shock emerging from coupling

  * •

The final state of the system $\dot{\phi}=0$ is rather strange from a
dynamical systems standpoint: this literally is a frozen oscillator. As such
the model does not converge to a well-defined attractor (contrary to both
Meinhardts’ model and the catastrophe-centric view from the initial Clock and
wavefront framework). One could for instance imagine in this model that even
after stopping, the cellular oscillators could be restarted in various ways by
taking control of the maturity $m$, or simply that local phases could be
shifted after stopping by some other processes.

###  4 Flow-based phase models

The LPM describes the PSM as a field of passive cells, swept by a frequency
gradient. But this picture does not fit the reality of growing embryos, where
progenitors divide in the tail bud to generate a flow of PSM cells moving from
posterior to anterior relative to the tail bud.

In 2001, Jaeger and Goodwin [140] accounted for this to propose a model based
on cellular oscillators recapitulating those developmental features. In this
model, cells exit the progenitor zone to enter the PSM, and their frequency
$\frac{2\pi}{T(\alpha)}$ depend on their ‘age’, $\alpha$, defined as the time
from their exit from the progenitor zone. The connection to the LPM model can
be made explicit by noticing that, in the moving frame of reference
(coordinate $m$), one has the following frequency profiles for $m>0$

$r(m)=\frac{2\pi}{T(m/v)}$ (18)

since $\alpha=m/v$ is the age of the cell at relative position $m$, assuming
growth occurs at speed $v$, and rescaling time units so that $r(0)=1$
corresponding to the taibud. A very similar model has been proposed by Kaern
et al.[141] within a more general framework of "Flow-Distributed Oscillator"
(FDO). This framework explicitly models a flow of oscillating cells injected
in a growth zone, so that the phase of cells depends in a simple way on their
age from the time of injection. One of the merits of those models is to more
explicitly connect the clock, the flow, and growth, thus providing some
explanations on observed phenotypes to explain some experimental data.

For instance, if one assumes that the period of a local oscillator depends on
its age, one can simply explain the observation that, if part of the PSM is
(artificially) reversed, it will pattern in a kinematic way, with an inverted
pattern relative to the normal one. Kaern et al.also use their flow model to
offer an explanation of the periodicity of heat shock phenotypes (that led to
the cell-cyle model). In the initial cell cycle model, to observe a periodic
pattern of disruption of heat shock, one has to assume that there is a
gradient of cell cycle phase within the PSM covering several cycles, so that
there are several clusters of cells at the heat-shock sensitive phase
$\phi_{S}$ separated by one cell cycle length (corresponding to 6-7 somites in
the chick). Kaern et al.observe that this explanation is unlikely because one
would then need a very long PSM to have up to 4 periodic anomalies due to a
single heat shock as observed experimentally. They rather assume that cells in
the progenitor zone have a uniform distribution of cell cycle phases and that
the heat shock delays mitosis. So at the moment of the heat shock, only a
fraction of the progenitors are about to go through mitosis, at phase
$\phi_{M}$ : heat shock delays them and as a consequence fewer cells divide,
and thus the PSM growth rate is impacted. Right after the heat shock,
progenitor cells slowly recover so that the growth rate transiently increases,
before coming back to the normal growth rate when all impacted cells have gone
through mitosis and other cells normally divide. In the absence of any
coupling/resetting, and assuming the cell cycle progresses uniformly in all
cells, the progenitor cells delayed by the first heat shock remain delayed, so
one would see again a decrease in the growth rate at the next cycle,
explaining how the defect can be repeated. In reality, one would expect some
coupling between cells and noise so that eventually the initial uniform
distribution of cell cycle phases in the progenitor is re-established, but if
this takes some time, one expects to observe an oscillatory growth rate of the
PSM post heat shock, with a period similar to the period of the cell cycle,
thus explaining the periodicity of the anomaly. One issue with this model is
that it does not account for the very first anomaly, which is clearly due to
cells already in PSM (and not in progenitor state): Kaern et al.suggest that
the initial heat shock also impacts the frequency of those cells, thus
implying a direct coupling between segmentation clock and cell cycle.

###  5 Delayed coupled models

We now consider elaborations of the LPM accounting for coupling delays between
oscillators (e.g. due to Notch/Delta signaling), presented in [142, 143, 144].
A line of oscillators (index $i$) is considered and the general time evolution
of the phase $\theta$ is given by

$\frac{d\theta_{i}}{dt}=\omega_{i}(t)+\epsilon_{i}(t)/2\sum_{k}\sin[\theta_{k}(t-\tau_{i}(t))-\theta_{i}(t)]+\zeta_{i}(t)$
(19)

$\omega_{i}(t)$ is a time/space dependent frequency, playing the exact same
role as $r$ in Eq/ 3. Morelli et al.[142, 143] consider the form

$\omega_{i}(t)=\omega_{\infty}(1-e^{-(i-vt)/\sigma})$ (20)

where $v$ is speed of the moving frequency gradient, so that $i-vt$ plays the
exact same role as the maturity $m$ in Eq. 3, except the spatial axis is
flipped so $i\rightarrow+\infty$ corresponds to posterior here and very
negative $i$ to anterior. For simplicity here, we first count length in
discrete units, but Morelli et al.introduce a unit conversion factor $a$.
$\sigma$ quantifies the spatial scale of changes of the frequency, similar to
Eq. 9. Notice that as $i\rightarrow-\infty$, $\omega_{i}$ would become very
negative, so one should rather assume a cut-off for $\omega_{i}$, for
instance, one can assume that $\omega_{i}$ is $0$ for $i<vt+f$, with $f\geq 0$
which gives a frequency at the front

$\omega_{f}=\omega_{\infty}(1-e^{-f/\sigma})$ (21)

Similarly, Morelli et al.define a cut-off maximum frequency $\omega_{N}$ in
the tailbud, defined at $i=N+vt$ so that :

$\omega_{N}=\omega_{\infty}(1-e^{-N/\sigma})$ (22)

So the PSM is entirely contained in the region $vt<i\leq vt+N$, and its length
is $N$ discrete units, or using a scaling factor $a$, a length $L=Na$ in
physical units.

The major difference with the LPM relies in the addition of coupling between
cells. The sum over $k$ accounts for nearest neighbors interactions, with
Kuramoto-like coupling [133]. Importantly, there is a delay in the coupling
($\tau_{i}(t)$), which depends on the local oscillator $i$. A noise
$\zeta_{i}(t)$ is also added for generality, but for now, let us put it to
$0$.

The first step is, similar to the LPM, to put oneself into the moving frame of
reference, i.e. calling $j=i-vt$, there is a fixed frequency gradient relative
to the moving frame of reference
$\omega_{j}=\omega_{\infty}(1-e^{-j/\sigma})$. We now define the phase
$\phi_{j}=\theta_{i}(t)$ in the moving frame of reference, which adds a
(discretized) drift term proportional to $v$ in the equation to get (assuming
now constant coupling and delays)

$\frac{d\phi_{j}}{dt}=\omega_{j}+v(\phi_{j+1}-\phi_{j})+\epsilon/2\left[\sin[\phi_{j+p+1}(t-\tau)-\phi_{j}(t)]+\sin[\phi_{j+p-1}(t-\tau)-\phi_{j}(t)]\right]$
(23)

One subtlety here is that $\phi_{j}$ represents the phase at position $j$ with
respect to the moving frame of reference, so does not represent the phase of
an actual, physical oscillator in this model. Rather it is a fictitious
oscillator representing the phase of the oscillators $\theta_{i}$ that we meet
at position $i=j+vt$ as we follow the frame of reference moving with speed
$v$. One of the major differences with the LPM is that one has to introduce an
integer $p=v\tau$, a length scale combining the delay and the speed of the
front. This correction accounts for the fact that during the delay $\tau$,
$\phi_{j}$ has physically moved by a distance $v\tau$ in the fixed frame of
reference, and thus is coupled to its physical neighbors at $t-\tau$, which
gives this non-local coupling for variable $\phi$, considerably complexifying
the analysis. Let us insist again that each physical oscillator $\theta_{i}$
does not change neighbors, it is only because of the motion of the frequency
gradient that the fictitious oscillator $\phi_{j}$ changes neighbors.

From there, the equation is not solvable in general but one can use classical
Ansatz similar to Eqs. 12,27, assuming the system globally oscillates and that
associated to it one observed a relative phase profile $\psi$ such that

$\psi_{j}=\phi_{j}-\Omega t$ (24)

$\Omega$ is the (yet unknown) global (constant) frequency of the process,
corresponding to the tail bud frequency. One gets by substitution into 23

$\frac{d\psi_{j}}{dt}=\Omega-\omega_{j}-v(\psi_{j+1}-\psi_{j})-\epsilon/2\left[\sin[\psi_{j+p+1}-\psi_{j}-\Omega\tau]+\sin[\psi_{j+p-1}-\psi_{j}-\Omega\tau]\right]$
(25)

The big advantage of this equation is that all time dependencies have now been
absorbed into the global frequency $\Omega$, so that the right-hand side of
Eq. 25 is a pure function of space (via the index $j$). We expect that at
steady state $\frac{d\psi_{j}}{dt}=0$ which gives :

$\Omega=\omega_{j}+v(\psi_{j+1}-\psi_{j})+\epsilon/2\left[\sin[\psi_{j+p+1}-\psi_{j}-\Omega\tau]+\sin[\psi_{j+p-1}-\psi_{j}-\Omega\tau]\right]$
(26)

This equation couples the global frequency of the process to the phase profile
$\psi_{j}$. Now the last step here, very similar to the LPM , is to consider
what happens in the tail bud, corresponding to the limit $j\rightarrow+\infty$
(in the notation of [143]). There, we assume that the tail bud is reached and
that there will be a homogeneous phase so that all $\psi_{j}$ are identical
with frequency $\omega_{N}$. This gives the self-consistency relation

$\Omega=\omega_{N}-\epsilon\sin(\Omega\tau)$ (27)

Solving the transcendental equation 27. Graphically, the root(s) of the
equation are found at the intersection of functions
$y=\omega_{N}-\epsilon\sin(\Omega\tau)$ and $y=\Omega$. One way to numerically
solve the equation is an iterative procedure: starting from an initial
condition $\Omega^{0}$, calculate
$\Omega^{k+1}=\omega_{N}-\epsilon\sin(\Omega^{k}\tau)$ until a desired
precision is reached. If there are multiple solutions, the procedure has to be
repeated with different initial conditions (i.e. uniformly spaced on the
relevant interval). The roots found are indicated with dots; $\epsilon=0.25$.
Eq. 27 defines self-consistently a collective frequency of somitogenesis,
which depends in a non-trivial way on the intrinsic tail bud frequency
$\omega_{N}$, on the coupling strength $\epsilon$ and on the delay $\tau$. It
is clear that if the delay $\tau$ or the coupling strength $\epsilon$ are very
small, $\Omega\sim\omega_{N}$ corresponds exactly to the LPM. So significant
deviations are expected for longer delays and stronger couplings. Equation 27
is a transcendental equation with typically several solutions (as illustrated
in Fig. 5). Assuming $\epsilon$ still is relatively small, we see that for
short delays $\tau$, the global frequency is lowered compared to $\omega_{N}$
(meaning that delays slow down the intrinsic oscillators). Increasing the
delay flips the sign of the sine (for delays roughly bigger than
$\frac{\pi}{\omega_{N}}$), there, the global frequency is increased, which
means that the delays now speed up the intrinsic tail bud oscillators. From
this theory, one expects a relative maximum change of period (compared to the
intrinsic period $\omega_{N}$) of about $20\%$ in both directions (depending
on the delay).

A continuous limit can also be obtained for this case, assuming a typical
distance $a$ between oscillators, calling $x=ia/L$ and defining
$\Psi(x)=\psi_{i}$ while taking the limit $a\rightarrow 0,N\rightarrow\infty$,
keeping $L=Na$ constant. While the discretized drift term
$v(\psi_{j+1}-\psi_{j})$ simply becomes $v\Psi^{\prime}(x)$, one must be
careful with the delay terms. The reason is that there is an intrinsic length
scale in the system $v\tau$ combining the speed and the delay, as can be seen
with the $p$ terms in Eq. 23. This distance should certainly not go to $0$ as
we take the limit $a\rightarrow 0$. More explicitly we have :

$\displaystyle\phi_{j+p\pm 1}-\phi_{j}$ $\displaystyle=$
$\displaystyle\Psi(x+v\tau+a)-\Psi(x)$ (28) $\displaystyle=$
$\displaystyle\Psi(x+v\tau\pm a)-\Psi(x+v\tau)+\Psi(x+v\tau)-\Psi(x)$ (30)
$\displaystyle\simeq$ $\displaystyle\Psi(x+v\tau)-\Psi(x)$ $\displaystyle\pm
a\Psi^{\prime}(x+v\tau)+\frac{a^{2}}{2}\Psi^{\prime\prime}(x+v\tau)$

From this, we can Taylor expand both $\sin$ terms in Eq. 26 close to
$\Psi(x+v\tau)-\Psi(x)-\Omega\tau$ to get the continuous self-consistent
equation for $\Omega$:

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\omega(x)+v\Psi^{\prime}(x)+\epsilon\sin[\Psi(x+v\tau)-\Psi(x)-\Omega\tau]$
(31) $\displaystyle-$ $\displaystyle\frac{\epsilon
a^{2}}{2}\sin[\Psi(x+v\tau)-\Psi(x)-\Omega\tau][\Psi^{\prime}(x+v\tau)]^{2}$
$\displaystyle+$ $\displaystyle\frac{\epsilon
a^{2}}{2}\cos[\Psi(x+v\tau)-\Psi(x)-\Omega\tau]\Psi^{\prime\prime}(x+v\tau)$

Notice the linear terms have canceled out so that the discretization has a
very small influence on $\Omega$. This formula has been derived in [145] for a
more general coupling function $h$ than the $\sin$ coupling assumed here (the
reasoning is completely equivalent). Again Eq.31 is highly nonlocal, coupling
phases at position $x$ and $x+v\tau$, so specific assumptions must be done to
solve this equation. However, like the discrete case, one can simply assume
that for $x\geq 1$ the phase profile is flat to account for the tailbud
dynamics so that all derivatives are $0$ and we get the exact same relation as
before

$\Omega=\omega(x=1)-\epsilon\sin(\Omega\tau)$ (32)

Also, just like one considered the limit $m\rightarrow\infty$ in the Lewis
Phase Model, one can consider the equivalent limit $x\rightarrow\infty$ in
this model where $\omega=0$. In the limit $a\rightarrow 0$ we get

$\Omega=v\Psi^{\prime}(x)+\epsilon\sin[\Psi(x+v\tau)-\Psi(x)-\Omega\tau]$ (33)

for which one obvious solution is $\Psi(x)=\frac{x\Omega}{v}$, indicating the
formation of a pattern of size $S=vT$, as expected from dimensional analysis
[145]. Numerical studies of the delayed models (continuous and discrete) can
be found in [143, 145].

Taking a step back, the phenomenology of this model, with traveling kinematic
waves, is qualitatively similar to the Lewis Phase Model, as can be seen from
the analogous mathematical treatments of both models. The Ansatz used are very
similar in spirit, allowing to derive a differential equation connecting the
stationary spatial phase profile to the moving frequency profile: Eq. 13 for
the Palmeirim model, and the more complicated Eq. 31 for the delayed model. To
be more quantitative in this comparison, it is useful to combine Eqs. 27 and
31 to write :

$v\Psi^{\prime}(x)=\omega(x=1)-\omega(x)+\epsilon\mathcal{H}\\{\Psi\\}$ (34)

where $\mathcal{H}$ is a complicated functional of the delay, $a$, $\Omega$
and of $\Psi$. Comparing with Eq. 13 , we see that $\omega(x=1)-\omega(x)$ is
completely equivalent to $1-r(m)$ in rescaled units (a minus sign comes from
the different directions of propagation of the fronts) and in particular, the
dominating term does not depend on the modified frequency $\Omega$. So all
influences of delays, couplings, etc… captured in the $\mathcal{H}$ term gives
only a (small) perturbation of order $\epsilon$ for the phase profile compared
to the much simpler Lewis Phase Model. Practically, it might be rather
difficult to see a difference between the theory with and without
coupling/delays by focusing only on the shape of the phase profile.

Conversely, the clear new prediction of the delayed model is the influence of
delay and coupling strength on the global frequency of the system, which is
amenable to experimental verification. Indeed Herrgen et al.[144] estimated
parameters of the delayed model from zebrafish data. Mutants of Notch pathway
were considered, as well as DAPT (a Notch inhibitor) treated embryos, and it
appears that for those mutants, clock period and segment length are increased
compared to WT, with a maximum change of around $10-20\%$. This is consistent
with the delayed theory in the sense that modifications of Notch pathway are
expected to change $\epsilon$, and the magnitude of the change is consistent
with the analytical theory. From there, saturating DAPT concentration was
assumed to put $\epsilon$ to $0$, allowing us to estimate that the uncoupled
period (corresponding to $\omega(x=1)$ in Eq. 27) was $18\%$ higher than the
observed segmentation period ($28$ min at $28$ C for zebrafish). This gives a
coupling strength $\epsilon=0.07\pm 0.04\mathrm{min}^{-1}$ and a minimal
coupling delay $\tau=21\pm 2\mathrm{min}$. Notice that the delay is very close
to the actual period of the oscillator. It is remarkable that the effects
predicted by delayed coupling theory are consistent in magnitude with
experimental data.

###  6 Doppler period shift

The phase-based models presented so far assume the existence of a moving
external frequency gradient, moving at speed $v$, which gives rise to a steady
state moving phase profile (defined by Eqs.13 for the LPM, and self-consistent
relation 31 for the delayed coupling theory).

It is, however, well-known that the size of the PSM slowly varies as a
function of time. This should have a direct observable consequence of speeding
up the segmentation period. To see this intuitively, one can first notice that
irrespective of the model considered, if the system at the embryo level is
truly periodic, the period of formation of new segments is exactly the same as
the period of the oscillation in the tailbud. But let us now assume that the
PSM starts shrinking, in the sense that anterior cells stop oscillating
earlier than expected for a fully periodic system. Then, since clock stopping
means segment formation, the next segment is expected to form earlier than in
a normal situation, thus the period of segmentation should decrease.
Alternatively, from the front standpoint, in the periodic case, the speed of
the moving frequency gradient (e.g. like in the LPM) is constant, so that an
observer at the front will meet a peak of the wave with the same frequency as
the tail bud oscillation. But if the front speeds up, an observer at the front
meets a peak of the wave more frequently. This means that the frequency of
segment formation is increasing, or that the segments will form with a shorter
period than the tail bud oscillation. This is reminiscent of the Doppler
effect in physics and has been studied experimentally [96]. As we will see,
biology is more complicated because it turns out that this Doppler effect is
partially compensated. To understand what happens we start with a dedicated
mathematical study detailed and expanded from [146, 147] (see also Appendix
B2) before we discuss experimental data and their implication.

####  1 Contributions to Doppler shift

We start with the expression phase profile in the moving frame of reference
with constant speed $v$ (equivalent to the LPM model, Eq. 11)

$\frac{\partial\phi}{\partial t}+v\frac{\partial\phi}{\partial
x}=\omega(x,L=\bar{x(t)})$ (35)

Here, we have defined the PSM length $L=\bar{x}$, which can change as a
function of time. We added a $L$ dependency in $\omega(x,L)$, which now
defines the frequency of the oscillator at position $x$ when the PSM has size
$L$. This will account for possible changes in the frequency profile with PSM
size. For instance, in [96, 147] it is assumed that $\omega(x,L)$ is a
function of the ratio $x/L$, i.e. the frequency gradient $\omega$ scales with
the PSM size. However, other functional forms might be possible so we will
stay generic for now. The tail bud is at $x=0$. For simplicity, we also assume
that for $x<0$ all cells are effectively synchronized, in particular,
$\left.\frac{\partial\phi}{\partial x}\right|_{x=0}=0$, and that
$\omega(0,L)=\omega_{0}$ does not depend on $L$, defining a constant reference
angular velocity/frequency. Notice that here, $x$ is the coordinate in the
moving frame of reference (corresponding to variable $m$ in Eq. 11).

The observed angular velocity at the front is by definition :

$\Omega_{A}=\frac{d}{dt}\phi(\bar{x}(t),t)=\left.\left(\frac{\partial\phi}{\partial
t}+\frac{d\bar{x}(t)}{dt}\frac{\partial\phi}{\partial
x}\right)\right|_{\bar{x}(t)}$ (36)

and the period of segmentation is $2\pi/\Omega_{A}$. If $\bar{x}$ does not
depend on time, we know from our previous Ansatz that
$\left.\frac{\partial\phi}{\partial t}\right|_{\bar{x}}=\omega_{0}$
corresponding to the frequency in the tail bud (See e.g. Eqs. 12 and 24). But
if $\bar{x}$ is changing, we can already see that there is one added
contribution coming from $\frac{d\bar{x}(t)}{dt}$, which is mathematically
similar to the Doppler effect in physics (see the end of this section for a
discussion of the common points and differences with physics). We will see
that other contributions arise, in particular, there will be a correction to
$\left.\frac{\partial\phi}{\partial t}\right|_{\bar{x}(t)}$.

Since the PSM size is changing, we can no longer use an Ansatz assuming there
is a moving phase profile with constant shape. Fortunately, it is possible to
explicitly integrate the equations, by considering a fixed position $z=x-vt$
in the static frame of reference. A cell at this position enters the PSM
($x=0$) at $t_{0}=-z/v$. So coming back to the absolute position $z$ we get
from Eq.35:

$\frac{\partial\phi(z,t)}{\partial t}=\omega(z+vt,\bar{x(t)})$ (37)

We directly integrate this equation keeping $z$ constant to get

$\displaystyle\phi(z,t)$ $\displaystyle=$
$\displaystyle\tilde{\varphi}(z)+\int_{t_{0}}^{t}\omega(z+vt^{\prime},\bar{x}(t^{\prime}))dt^{\prime}$
(38) $\displaystyle=$
$\displaystyle\varphi(t-\frac{x}{v})+\frac{1}{v}\int_{0}^{x}\omega\left(x^{\prime},\bar{x}\left(t-\frac{x-x^{\prime}}{v}\right)\right)dx^{\prime}$

where $\tilde{\varphi}(z)$ is the initial phase of the cell as it enters the
PSM and $\varphi$ the same function expressed as a function of
$(t-\frac{x}{v})$. It is convenient to simplify this term by immediately using
the boundary condition $\left.\frac{\partial\phi}{\partial x}\right|_{x=0}=0$,
meaning a flat phase gradient in $0$. We have from 38

$0=\left.\frac{\partial\phi}{\partial
x}\right|_{x=0}=-\frac{1}{v}\dot{\varphi}+\frac{1}{v}\omega(0,\bar{x(t)})$
(39)

so that

$\dot{\varphi}=\omega(0,\bar{x(t)})=\omega_{0}$ (40)

i.e. as expected the tail bud simply oscillates with frequency $\omega_{0}$,
and we thus have (taking as an arbitrary initial condition $\varphi(0)=0$)

$\displaystyle\phi(x,t)$ $\displaystyle=$
$\displaystyle\omega_{0}(t-\frac{x}{v})+\frac{1}{v}\int_{0}^{x}\omega\left(x^{\prime},\bar{x}\left(t-\frac{x-x^{\prime}}{v}\right)\right)dx^{\prime}$
(41) $\displaystyle=$
$\displaystyle\omega_{0}t+\frac{1}{v}\int_{0}^{x}\omega\left(x^{\prime},\bar{x}\left(t-\frac{x-x^{\prime}}{v}\right)\right)-\omega_{0})dx^{\prime}$
(42) $\displaystyle=$ $\displaystyle\omega_{0}t+\psi(x,t)$ (43)

defining

$\displaystyle\psi(x,t)=\frac{1}{v}\int_{0}^{x}\omega\left(x^{\prime},\bar{x}\left(t-\frac{x-x^{\prime}}{v}\right)\right)-\omega_{0})dx^{\prime}$
(44)

Eq.42-43 are equivalent to the now familiar Ansatz defined e.g. in Eqs. 12,
24, since we see from Eq. 43 that $\psi(x,t)$ is the phase profile relative to
the tail bud phase. In previous derivations, this phase profile did not depend
on time, e.g. integrating Eq. 13 of the LPM with proper units we had in a
similar way :

$\displaystyle\psi(x)=\frac{1}{v}\int_{0}^{x}(\omega(x^{\prime})-\omega_{0})dx$
(45)

The big difference is that $\omega$ now depends on time in a very complex way:
we see from Eq. 44 that the relative phase profile $\psi$ depends on all past
PSM lengths $\bar{x}\left(t-\frac{x-x^{\prime}}{v}\right)$. So this problem is
highly non-local; below we will see some simplifications that can arise for
special forms of $\omega$. It is useful to define the time derivative of this
phase profile

$\Lambda(x,t)=-\frac{\partial\psi(x,t)}{\partial t}=-\frac{\partial}{\partial
t}\frac{1}{v}\int_{0}^{x}\omega\left(x^{\prime},\bar{x}\left(t-\frac{x-x^{\prime}}{v}\right)\right)dx^{\prime}$
(46)

Notice that as expected, $\Lambda$ is $0$ for constant $\bar{x}$

This allows us to simply compute from Eq. 43

$\frac{\partial\phi}{\partial t}=\omega_{0}-\Lambda(x,t)$ (47)

and

$v\frac{\partial\phi}{\partial
x}=\omega(x,\bar{x}(t))-\omega_{0}+\Lambda(x,t)$ (48)

[Here we used the fact that the integrated $\bar{x}$ term in Eq. 44 depends on
$\left(t-\frac{x-x^{\prime}}{v}\right)$ so that for this term
$\partial_{t}\bar{x}=-1/v\partial_{x}\bar{x}$, which allows to express
$\frac{\partial\Psi}{\partial x}$ as a function of $\Lambda$ ]

Now injecting this into Eq. 36 we get :

$\displaystyle\Omega_{A}$ $\displaystyle=$
$\displaystyle\omega_{0}-\Lambda(\bar{x},t)+\frac{\dot{\bar{x}}}{v}\left[\omega(\bar{x},\bar{x})-\omega_{0}+\Lambda(\bar{x},t)\right]$
(49) $\displaystyle=$ $\displaystyle\Omega_{P}+\Omega_{W}+\Omega_{D}$ (50)

defining

$\Omega_{P}=\omega_{0}$ (51)
$\Omega_{W}=-\Lambda(\bar{x},t)=\left.\frac{\partial\psi(x,t)}{\partial
t}\right|_{x=\bar{x}}$ (52)
$\Omega_{D}=\frac{\dot{\bar{x}}}{v}\left[\omega(\bar{x},\bar{x})-\omega_{0}+\Lambda(\bar{x},t)\right]=\dot{\bar{x}}\left.\frac{\partial\psi(x,t)}{\partial
x}\right|_{x=\bar{x}}$ (53)

There are three more or less intuitive contributions in this sum. The first
contribution is $\Omega_{P}=\omega_{0}$ : this is the usual angular frequency
in the tail bud (posterior), which is the "standard" term when the PSM size is
constant.

$\Omega_{D}$ is the so-called "Doppler" effect. It is proportional to the
shrinkage speed of the PSM, and composed of two terms. First,
$\omega(\bar{x},\bar{x})-\omega_{0}$ quantifies the difference in angular
frequency between front and tail bud. The contribution to the Doppler effect
is rather intuitive: since cells at $\bar{x}$ oscillate with angular frequency
$\omega(\bar{x},\bar{x})$, within a small time interval $dt$, there is a
(negative) phase difference accumulation relative to the tail bud equal to
$(\omega(\bar{x},\bar{x})-\omega_{0})dt$. Now when the PSM shrinks by a small
quantity $d\bar{x}$, this defines a small time $d\tau=d\bar{x}/v$ for which
those cells do not accumulate this phase difference, and thus, per unit of
time a relative increase of phase
$\frac{\dot{\bar{x}}}{v}(\omega(\bar{x},\bar{x})-\omega_{0})$, which
corresponds to an increase of angular velocity as expected. So this term
exactly is the intuitive Doppler contribution described in the preamble of
this section.

The less intuitive contributions come from
$\Lambda(x,t)=-\frac{\partial\psi(x,t)}{\partial t}$ both independently in
$\Omega_{W}$ and as a term in the Doppler term. By definition, $\Lambda(x,t)$
is the change of the phase profile due to the shrinking of the PSM. Since the
phase profile is directly related to the wavelength of the pattern, it is
called the "dynamical" wavelength in [146, 147]. $\Lambda$ is particularly
difficult to compute (analytically) since it contains information about the
entire history of PSM shrinkage, so extra assumptions are needed to compute
it. In the Appendix B2, we compute $\Lambda$ for two cases: when the frequency
gradient is instantaneously scaling with PSM size (frequency scaling), which
is the hypothesis made in [96, 146, 147], and when the frequency gradient is
left unchanged and a front of $0$-frequency moves towards the anterior without
changing it (frequency cropping).

Experimentally, Soroldoni et al.[96] visualized the segmentation clock
dynamics in a zebrafish embryo. Using moving Regions of Interest (ROIs), they
monitored oscillations both in the posterior and anterior end of the PSM for
almost 20 cycles and noticed that for 9 oscillations in the posterior, one
observes 10 oscillations in the anterior, suggesting a non-stationary process.
The period of segment formation is however the same as the period of the
oscillator in the anterior, suggesting they are indeed the same process.

The difference between anterior and posterior periods is accompanied by a
shrinkage of the PSM by about 60 % while 13 segments are formed, allowing to
quantify the Doppler shift. A phase map is then experimentally derived to
visualize the profile $\psi(x,t)$. It is very clear that these profiles change
with time: for instance, over 500 mins (corresponding to 12 oscillations in
the posterior), the phase difference $|\psi(\bar{x},t)-\psi(\bar{x},0)|$
decreases from more than $5\pi$ to roughly $3\pi$ (Fig S5 in [96]). It is
quite noticeable though that there are different regimes for $\psi$, for
instance by visualizing $\partial\psi/\partial t$, one sees that initially,
$\partial\psi/\partial t$ (corresponding to the term $\Lambda$) is non-zero in
a broad anterior region of the embryo, before reaching $0$ later on after
$300$ mins. Quantitatively, the posterior frequency $\Omega_{P}$ can be
directly measured from the moving ROI in the posterior, and is around
$0.15min^{-1}$ for most of the time, slowly decreasing to roughly
$0.12min^{-1}$. To estimate other contributions, one has to experimentally
measure $\partial\psi/\partial t$ and $\partial\psi/\partial x$, which can be
done from kymographs of the oscillations. The Doppler contribution
$\Omega_{D}=\dot{\bar{x}}\left.\frac{\partial\psi(x,t)}{\partial
x}\right|_{x=\bar{x}}$ varies with time between $0.025min^{-1}$ and
$0.05min^{-1}$, so quite significantly compared to $\Omega_{P}$. The Dynamic
wavelength $\Omega_{W}=\left.\frac{\partial\psi(x,t)}{\partial
t}\right|_{x=\bar{x}}$ is of the order of $-0.03min^{-1}$ and gets to $0$
after 300 mins. Those values are summarized on Fig. 3.

So the Doppler contribution and the Dynamic wavelength are rather high in
magnitude (roughly $1/3$ of the posterior frequency) but almost compensate so
that the segmentation period is only roughly 10 % shorter than the posterior
period. To understand how this compensation happens, going back to Eq. 50, one
can make the following assumptions from the data [147]:

  * •

a constant shrinking PSM rate
$\beta=\frac{\dot{\bar{x}}}{v}=\frac{\bar{v}}{v}$

  * •

a ‘scaling’ frequency gradient with PSM size, i.e.
$\omega(x)=\omega_{0}U(x/L)$ where $L$ is the current size of the PSM

One then gets a more compact expression for $\Omega_{A}$ :

$\Omega_{A}=(1+\beta)(1-\Delta)\omega_{0}$ (54)

introducing $\Delta=\beta\int_{0}^{1}U(x)(1+\beta\xi)^{-2}d\xi$ (see Appendix
B)

We see better what happens from this expression. There are two effects. The
$1+\beta$ corresponds to a ‘traditional‘ Doppler effect. In the Appendix we
show that this contribution would exactly be the Doppler effect corresponding
to a shrinking PSM with a $0$ frequency in the anterior, thus shortening the
period. However, there is an added factor $(1-\Delta)$, which clearly goes in
the other direction, i.e. lengthening the period. This is due to the dynamical
wavelength effect [96], originating in this model from the frequency scaling
[147]. The origin of this effect is a bit easier to understand from the
frequency standpoint: if the frequency profile scales with the PSM size, in
means that as the front moves, the frequency observed right at the front is
actively lowered compared to a situation where the frequency gradient is held
constant (see e.g. frequency cropping calculation in Appendix). As noticed in
[147], a similar effect can be obtained for a fixed observer in classical
waves equation with a refractory index increasing with time, which gives local
frequency going to $0$ in the limit $t\rightarrow\infty$.

Experimentally, $\beta$ approximately is $1.15$, meaning that the PSM shrinks
as fast as the tail bud grows. That means that without frequency gradient
scaling, we would (initially) expect $\Omega_{A}$ of the order of
$2\omega_{0}$, much faster than what we observe. But, fitting the frequency
gradient with a functional form
$U(x)=\sigma+(1-\sigma)\frac{1-e^{k(x-1)}}{1-e^{-k}}$ for $0\leq x\leq 1$, one
finds experimentally $\sigma=0.34$ and $k=2.07$, which gives $\Delta\sim
0.45$, so that the $(1-\Delta)$ factor almost exactly compensates $(1+\beta)$
111It is also interesting to see that $\Delta$ depends only very weakly on the
shape of the gradient, varying $\sigma$ gives values of $\Delta$ between $0.4$
and $0.5$.

In physics, compensation of big numbers to give relatively smaller ones are
often called "unnatural" because they suggest some arbitrary parameter fine-
tuning. Here, the observed compensation might in fact suggest an active
scaling mechanism so that indeed $\omega(x,L)$ is a pure function of $x/L$,
the relative position of a PSM cell, as hypothesized in [147] to fit the data.
Such a hypothesis further accounts for some nontrivial features of
segmentation, for instance, the segment length reaches a maximum for some
intermediate segment numbers [147], before shrinking down at a later time.

####  2 Doppler effect or Doppler period shift?

In physics, the Doppler effect is associated with a wave propagating towards a
moving observer, or, similarly, a wave source moving relative to a static
observer. This is the familiar experience of an ER vehicle moving rapidly in
the street: when the vehicle moves towards an observer, the pitch of the sound
is higher than when it moves away from the observer. Another example in
physics is the famous "red shift" perceived in light emitted from galaxies
further away, indicating that they are moving away from us (and demonstrating
universe expansion). More quantitatively: assuming a standard wave equation in
one-dimension with constant speed $c$ given by:

$\frac{\partial^{2}u}{\partial t^{2}}=c^{2}\frac{\partial^{2}u}{\partial
x^{2}}$ (55)

A propagating solution from left to right with angular velocity $\omega$ takes
the form

$u(x,t)=u_{0}\sin(\omega(t-x/c)+\phi_{0})$ (56)

defining a phase

$\phi(x,t)=\omega(t-x/c)+\phi_{0}$ (57)

Now assuming an observer moves as a function of time, position $\bar{x}(t)$,
the observer sees the phase $\bar{\phi}(t)=\phi(\bar{x}(t),t)$ and the
corresponding observed frequency $\Omega$ is :

$\displaystyle\Omega=\frac{d\bar{\phi}}{dt}$ $\displaystyle=$
$\displaystyle\left.\left(\frac{\partial\phi}{\partial
t}+\frac{d\bar{x}(t)}{dt}\frac{\partial\phi}{\partial
x}\right)\right|_{\bar{x}(t)}$ (58) $\displaystyle=$
$\displaystyle\omega(1-\frac{1}{c}\frac{d\bar{x}(t)}{dt})$ (59)

We indeed see that if $\frac{d\bar{x}(t)}{dt}<0$, meaning that the observer is
moving towards the wave, $\Omega>\omega$ (corresponding for instance to the
higher ER vehicle pitch), and conversely is lower if the observer is moving
further from the wave. Dynamical wavelength effects can also be observed in
physics if $\omega$ in Eq. 55 is time-dependent (e.g. a time-dependent
refraction index) as discussed in [147].

So the mathematical formalism is clearly similar, however, there is an
important difference in terms of the actual meaning of the equations,
especially in the biological context. In the classical examples of Doppler
effects mentioned above (ER vehicle moving, redshift), there is a local
observer moving relative to the source of the signal and ’sensing’ the
frequency $\Omega$. In a biological context, there is a priori no such local
observer moving with respect to the source of the oscillation (e.g. the
tailbud): rather the front is a macroscopic variable, defined by a moving
location $\bar{x}$, that might only be defined on average, see e.g. the
complicated motion of the bifurcation in the Meinhardt model. Remarkably
though, there is an associated macroscopic observable at the embryo level: the
frequency of segment formation, which is measured by the experimentalist
(again with multiple caveats, since e.g. boundary formation is a discrete
process). So maybe one could rather refer to this phenomenon as a Doppler
segmentation period shift, indicating this is an emergent effect at the embryo
level to be contrasted with a Doppler effect usually tied to an observer
moving relative to the wave.

Figure 3: Contributions to the Doppler Shift of Segmentation period in
zebrafish

###  7 Wavefront as a phase shock emerging from coupling

In all phase models described so far, the oscillators are controlled by an
imposed frequency gradient and associated front, similar to the original clock
and wavefront model. Biologically, this could correspond to a morphogen
gradient (such as FGF) controlling clock frequency and imposing a threshold
for clock stopping. However, recent models have suggested different views
where both the observed frequency gradient and the clock stopping are much
more connected to the clock itself.

A remarkably simple model has been proposed by Murray et al.[139, 148] with
spatially coupled phases :

$\frac{\partial\theta}{\partial t}=\omega+A\frac{\partial^{2}\theta}{\partial
x^{2}}-B\left(\frac{\partial\theta}{\partial x}\right)^{2}$ (60)

This model can be derived as a Taylor expansion of oscillators spatially
coupled with a Kuramoto-Sakaguchi coupling of the form
$H(\Delta\theta)=a\sin(\Delta\theta)+b(\cos\Delta\theta-1)$ : the coupling
term at position $x$ would then be (Taylor expanding at order $2$ in $dx$)

$\displaystyle H(\theta(x+dx)-\theta(x))+H(\theta(x-dx)-\theta(x))$
$\displaystyle=$ $\displaystyle
a\left(\sin(\theta(x+dx)-\theta(x))+\sin(\theta(x-dx)-\theta(x))\right)$
$\displaystyle+$ $\displaystyle
b\left(\cos(\theta(x+dx)-\theta(x))-1+\cos(\theta(x-dx)-\theta(x))-1\right)$
(63) $\displaystyle\sim$ $\displaystyle
a(\theta(x+dx)+\theta(x-dx)-2\theta(x))$
$\displaystyle-\frac{b}{2}((\theta(x+dx)-\theta(x))^{2}+(\theta(x-dx)-\theta(x))^{2}$
$\displaystyle\sim$
$\displaystyle\frac{a}{2}dx^{2}\frac{\partial^{2}\theta}{\partial
x^{2}}+bdx^{2}\left(\frac{\partial\theta}{\partial x}\right)^{2}$ (64)

so that we see that $A=\frac{a}{2}dx^{2}$ and $B=bdx^{2}$.

$\omega$ represents the frequency in the posterior part of the PSM. The $A$
term is a classical phase diffusion term. The novelty arises from the
quadratic $B$ term. If $B>0$, we see that any local phase gradient tends to
slow down the local oscillator. So imagine that we have a region with a sharp
phase gradient next to a region with a flat phase profile: we would then
expect that the oscillators at the boundary between the two regions would get
very slow, thus expanding the region of high phase gradient (see Fig. 4). In
other words, regions with high phase gradients would tend to propagate,
suggesting a plausible mechanism for wavefront formation.

Figure 4: Simulations of the Murray model : kymographs on the left and phase
profiles for different times on the right.

To be more quantitative, Eq. 60 clearly has a steady state solution, that we
interpret as the pattern in the anterior part of the embryo. Assuming that the
phase is stationary and linear spatially as $x\rightarrow-\infty$ (so
corresponding to anterior, like in the delayed phase model from [143]), we get

$\left.\frac{\partial\theta}{\partial
x}\right|_{x\rightarrow-\infty}=\sqrt{\frac{\omega}{B}}$ (65)

In the posterior part of the embryo, we assume synchronous oscillations, so
that

$\left.\frac{\partial\theta}{\partial x}\right|_{x\rightarrow+\infty}=0$ (66)

Defining as a new variable the gradient of the phase
$\Psi=\frac{\partial\theta}{\partial x}$ we see that $\Psi$ has to interpolate
between $0$ and a finite value. In fact, as intuitively explained above, one
will observe a propagation of the region with non-zero $\Psi$, towards the
region of $0$ $\Psi$, corresponding to the stabilization of the phase gradient
giving rise to a periodic pattern! More precisely we get for $\Psi$ :

$\frac{\partial\Psi}{\partial t}+2B\Psi\frac{\partial\Psi}{\partial
x}=A\frac{\partial^{2}\Psi}{\partial x^{2}}$ (67)

This equation is well-known in physics and is called Burgers’ equation [149].
The left-hand side of this equation looks similar to standard advection
equations such as Eq. 11 or Eq. 35 but with a speed proportional to the value
of the propagating field $\Psi$. The Burgers’ equation has been applied to
various contexts, e.g. traffic flow modeling, or boundary propagation and
there are standard methods so solve it analytically [150]. A stochastic
version of Burgers’ equation is the famous Kardar-Parisi-Zhang model of
stochastic growth of surface [151], which is also related to stochastic models
of particle motions in 1D such as the asymmetric exclusion process [152].
Burgers’ equation is clearly non-linear, but can still be solved using the
same stationary Ansatz as before :

$\Psi(x,t)=w(x-vt)=w(m)$ (68)

assuming a propagation of the pattern from left to right with (yet unknown)
speed $v$. We then get

$-vw^{\prime}+2Bww^{\prime}=Aw^{\prime\prime}$ (69)

which directly integrates into an order one differential equation

$-vw+Bw^{2}=Aw^{\prime}+C$ (70)

Now we use the boundary condition that $w=w^{\prime}=0$ at $x=+\infty$ and
$w^{\prime}=0,w=\sqrt{\frac{\omega}{B}}$ at $x=-\infty$. This directly gives
$C=0$ and the value of the speed $v=\sqrt{\omega B}$. With the value of $v$ it
is straightforward to integrate Equation 70 and one gets immediately (properly
shifting the initial value for $m$)

$w(m)=\frac{\sqrt{\omega/B}}{1+e^{\sqrt{\omega B}m/A}}$ (71)

Now $\Psi=\partial\theta/\partial x=w(m)$ with $m=x-vt$ so that we recover the
simple equation

$\frac{\partial\theta}{\partial m}=\frac{1}{v}\frac{\omega}{1+e^{vm/A}}$ (72)

where we made the dependency on $v=\sqrt{\omega B}$ explicit. Rescaling time
and space to that $\omega=v=1$, we then get

$\frac{\partial\theta}{\partial m}=\frac{1}{1+e^{m/A}}$ (73)

which is completely equivalent to Eq. 13 with

$r(m)=\frac{1}{1+e^{-m/A}}=\frac{1}{1+e^{(-x+t)/A}}$ (74)

This exactly is the frequency profile assumed in the LPM Eq. 4 ! (modulo a
flipping of the direction of the $x$ axis). Thus this model is
indistinguishable from the LPM by only considering stationary propagation of
the front.

To really contrast this model with the LPM, one should test the mechanism
underlying front propagation related to the $B$ term, in particular its
influence on the front speed. The first interesting prediction relates to the
scaling of the pattern with the clock period [139, 148]. Since $v=\sqrt{\omega
B}$, and the somite size $S=2\pi v/\omega$, we get :

$S=2\pi v\sqrt{B/\omega}=\sqrt{2\pi B}\sqrt{T}$ (75)

where $T$ is the period of the tail bud oscillator. So if we were to modify
the period, one would get a scaling law of power $1/2$ for the somite size
with respect to the period, while if speed and period vary independently we
expect $S$ and $T$ to be proportional, so rather a power $1$.

Another way to test the model is to directly change $B$. In particular, if we
locally decrease $B$, we expect a locally smaller speed, and thus a pattern
with a smaller wavelength. Conversely, if $B$ is locally increased, the front
speed is higher and the pattern should have a bigger wavelength. As said
above, we expect coupling to be related to Notch signaling, so that those
predictions on the $B$ term could be tested by performing transplants of cells
with different Notch levels. Such transplants have been done in [153], where
cells expressing high levels of Delta ligands induce locally smaller somites,
so consistent with a smaller $B$ term. However, we notice that those
experiments are difficult to interpret in the present context without more
direct connections to the $A$ and $B$ terms. In particular, since there are
two coupling terms $A$ and $B$, with $A$ acting to keep oscillators in
synchrony while $B$ tending to increase phase gradient, it is unclear what
would be the net effect of the modification of Notch signaling pathway. But
those experiments combined with this modeling have immense merit to point out
the potential role of Notch lateral inhibition in destabilizing the phase
gradient to generate a wavefront and stabilize the pattern.

Experimentally, one can estimate relatively simply the different parameters
$\omega,A,B$ by connecting them to the physical observables we already
discussed, respectively the tail bud frequency, the size of the transition
zone from oscillation to dynamic pattern, and the speed of the front/size of
the pattern. Based on data from [76], $\omega$ obviously is the angular
velocity in the tail bud, corresponding to $0.21{min}^{-1}$ in zebrafish [139,
148]. Since $B$ is related to the speed $v$ of the front and $\omega$, it can
be estimated from somite size $S=vT$ to get $B=S^{2}/{2\pi T}$: with a somite
length of roughly 6 cell diameters ($cd$), one finds that
$B=0.19cd^{2}.min^{-1}$. The most difficult parameter to estimate is $A$,
which relates to the typical length scale of the phase gradient $L_{exp}$.
Using a linear approximation, one finds $L_{exp}=4A/\sqrt{\omega B}$, and
experimentally one finds that $L_{exp}\sim 48cd$, which gives $A\sim
2.4cd^{2}min^{-1}$. Murray et al.[139, 148] also inferred parameter values for
different species using published data [108], in particular, rescaling the
period to generate the dimensionless equation

$\frac{\partial\theta}{\partial t}=\gamma\frac{\partial^{2}\theta}{\partial
t^{2}}-16\gamma^{2}\left(\frac{\partial\theta}{\partial t}\right)^{2}+1$ (76)

with $\gamma=S/4L_{exp}=\pi B/8A$. The rescaled speed is $4\gamma$, and the
number of waves in the PSM can be computed from Eq.16 and is $1/8\gamma$. One
then finds that mouse and chick have relatively close $\gamma$ of around
$0.06$, giving two waves, for zebrafish $\gamma=0.04$, giving slightly more
than 3 waves and for snake $\gamma=0.02$ giving around 7 waves.

###  8 Phase-amplitude coupling and excitability for oscillation arrest

In [97], we observed that, within one future presumptive somite, the apparent
wave of clock arrest does not move smoothly from anterior to posterior part,
as expected in classical clock and wavefront or phase models (such as the LPM
or the Murray model from the previous section). Rather, the Notch wave sweeps
from posterior to anterior, so that it eventually stops in the anterior part
of each somite. This means that, as measured in terms of Notch signaling
maximum activity, the oscillation cycle dies out first in the posterior part
of a given somite and not in its anterior part where the front should first
pass. We named this process "clockwave stopping", to stress out the motion of
that last wave. So there is some interesting paradox between scales here,
suggesting that phase dynamics alone is not enough to account for the
difference between the local (somite, posterior to anterior) scale and global
(tissue, anterior to posterior) dynamics of clock stopping.

A possible explanation is that the clock stopping is coupled to some other
observable moving in the same direction as the oscillation wave, from
posterior to anterior. For instance, it is well known that the amplitude of
the cellular oscillations increases from posterior to anterior, as described
in [97], Fig. 4. To explain the local vs global paradox in the direction of
the clock stopping, it is then natural to assume that the overall oscillatory
signal (combining phase and amplitude) controls some bifurcation leading to
clock arrest.

This can be captured by the following simple phenomenological model, that was
used in [97], Fig. 4. Adding a variable $A$ for amplitude and defining the
overall clock signal $s(A,\theta)=A(1+\cos\theta)$, we write in the PSM frame
of reference:

$\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle\omega(m)$ (77)
$\displaystyle\dot{A}$ $\displaystyle=$ $\displaystyle\Theta(m)\lambda_{0}A$
(78)

As usual, $\omega(m)$ accounts for the slowing down of the clock with variable
$m=x+t$ in rescaled units. The precise functional form does not matter as long
as $\omega$ goes to $0$ for big $m$. $\Theta(m)$ is the Heaviside function
equal to $1$ if $m>0$ and $0$ otherwise. So for tail bud cells ($m<0$),
$\lambda(m)=\Theta(m)\lambda_{0}=0$ and $A$ stays constant (and non zero).
Then when cells are injected in the PSM (i.e. $m>0$), $\lambda=\lambda_{0}$,
the amplitude increases in PSM (as the clock is slowing down). So those two
equations simply capture both the clock slowing down and the amplitude
increasing. To account for the clockwave stopping, we impose the following
stopping conditions :

  * •

as soon as $s(A,\theta)>A^{*}$, the bifurcation is crossed

  * •

once the bifurcation is crossed, the oscillation completes its cycle up to
phase $\theta=\pi$

So the primary wave here is defined by the condition $s(A,\theta)>A^{*}$. The
second condition phenomenologically reproduces the situation after a
bifurcation from an oscillating to an excitable system, such as the system on
the right-most panel of Fig. 9 (see Appendix for discussion on excitable
systems). The idea is that a new fixed point appears on (or very close) to the
cycle, because of relaxation-like dynamics the overall oscillations need to be
fully completed before the system stabilizes.

Fig. 5 shows a simulation of this simple model, which strikingly recapitulates
many features of the zebrafish clock arrest, in particular a "sawtooth"
pattern at the front, and the maximum intensity in the anteriormost part of
each somite, corresponding to the location where the clock stops last at each
wave. We represent dynamics in three cells in Fig. 5 D : green is a posterior
somite cell of a given somite, while blue and orange respectively are
posterior and anterior cells of the next somite. Notice that the blue and
orange cells go through an additional cycle compared to the green cell and
that the orange cell, despite being very close to the green one, cycles once
more with higher intensity. This possibly solves the Meinhardt paradox between
AP and PA boundaries: from the orange cell standpoint (anterior A), there is
an additional cycle in the blue cell vs the green one, both being committed to
posterior fates. So because of this asymmetric wave dynamics, one would simply
need a mechanism to coordinate the epithelialization of cells locally going
through the exact same number of oscillations to define a proper somite. More
generally, this suggests that excitable dynamics close to the front might
indeed play a role, not only to explain the wave pattern but also to control
the bifurcation itself. In fact, as will be described in section Landscape
geometry of segment formation, such bifurcation could also be more directly
coupled to the slowing down dynamics itself. We will discuss further evidence
for such excitability of the segmentation oscillators in section Hacking the
segmentation clock.

Figure 5: Simulation of Eqs. 78-77, showing a sawtooth pattern at the front
and recapitulating the posterior to anterior stopping dynamics in a given
somite. (A) Simulated Phase (B) Simulated Amplitude (C) Total Signal. Notice
the strong asymmetry and the "sawtooth" pattern. Green, orange, and blue lines
correspond to positions displayed on panel B, the green cell is assumed to be
one somite earlier than the blue and orange ones (D) Simulated oscillations in
three cells indicated in (C).

##  Chapter 4 From systems biology to geometric models

In our discussion so far, we have considered relatively abstract models, e.g.
focused on the dynamics of the phase in the previous section. But many "low-
level" components (genes, coupling effects, possibly mechanical effects) are
driving the segmentation process.

With the advances in systems biology (both theoretical and experimental), it
is possible to detail more explicitly processes implicated in segmentation. It
would be a tremendous achievement to follow a bottom-up approach, identifying
precisely all components of the system, modeling them, and from there
inferring the observed dynamics. We are still not there, and the historical
approach has been much more top-down, with at least two modeling steps :

  * •

relate phase dynamics to possible biochemical mechanisms

  * •

relate possible biochemical mechanisms to actual genes

We make here an explicit distinction between biochemical mechanisms and genes
because the same biochemical process (defined in a loose sense, say for
instance a negative feedback oscillator, or a bistable system) could possibly
be implemented with different genes in different species, explaining
interspecies variability mentioned earlier [67]. This can lead to confusion
(common to many systems biology modeling work) on what is actually modeled and
predicted. For instance, we have already seen in the previous section that
very similar phase dynamics can be modeled with different biochemical
processes (e.g. via spatial coupling terms), so a detailed gene-based model
reproducing only the correct wave pattern might not necessarily be insightful
without specific experimental tests. This raises a common problem for all
models based on explicit gene dynamics: they can be falsified on at least
three levels (the genes themselves, the biochemical process they implement,
and the emerging phase dynamics). For this reason, in this section, we focus
on models exploring specific biochemical mechanisms inspired by actual genetic
interactions, eventually circling back onto "geometric models" not tied to
specific genes.

###  1 Simple delayed oscillators

Many models of segmentation clocks are based on so-called negative feedback
oscillators. There is a large (experimental and theoretical) literature on
such oscillators in a broader systems biology context, notably for circadian
clocks modeling [154, 155, 156] and cell cycle [157], see also [29] for a
recent review on developmental timing.

Many molecular networks implicated in somitogenesis could plausibly implement
such negative feedback. For instance, her1 and her7 in zebrafish appear to
both oscillate and negatively regulate their own expression [74, 75] as well
as the ortholog of her1 in mouse, Hes7 [158, 159]. Hes1 in mouse is also
clearly implicated in a negative feedback loop [77].

When modeling negative feedback loop models to get oscillations, it is
tempting to write a system of RNA/protein regulating its own production in the
following way (using notations from [78, 160]):

$\displaystyle\frac{dm}{dt}$ $\displaystyle=$ $\displaystyle f(p)-cm=u(m,p)$
(1) $\displaystyle\frac{dp}{dt}$ $\displaystyle=$ $\displaystyle am-bp=v(m,p)$
(2)

with $f^{\prime}(p)<0$ to model negative self-interaction. However, it can be
proved that such a system can not oscillate: assume there is a closed orbit
$S$ in 2D, then the flux
$\oint_{S}(u(m,p)\mathbf{m}+v(m,p)\mathbf{p})\cdot\mathbf{n}\mathrm{d}s$ is
$0$ since $\mathbf{n}$ is always normal to the tangent of the orbit. But from
the divergence theorem, this integral is also equal to
$\iint\limits_{R}(f^{\prime}(p)-b)\mathrm{d}m\mathrm{d}p$ which is strictly
negative inside the domain $R$ limited by $S$, so there is a contradiction,
and thus there can not be any closed orbit. This simple reasoning is called
Poincaré-Bendixson criterion (See discussion in [161]).

There are several ways to fix this. For instance, negative feedback loops in
gene regulatory networks are not necessarily direct and could be via extra
steps (e.g. a phosphorylation cascade, a spatial segregation step, or an extra
transcriptional step where a gene activates its own repressor instead of
repressing itself directly). Mathematically, this increases the dimensionality
of the system which is no longer limited by the Poincaré-Bendixson criterion
but also becomes more complicated to study. Examples of such models include
the Goodwin model [162], or the MFL model [163], both of which are still
simple enough to be studied analytically. Alternatively, the addition of
positive feedbacks allows us to circumvent the Bendixson criterion to get
oscillations with a 2D model (see e.g. [164]). Coming back to the segmentation
oscillator, direct binding of Hes7 protein to its own promoter has been shown
[159], excluding the presence of a clear intermediate step and there is no
evidence either for positive feedback in this part of the network. A
parsimonious modeling choice rather is to assume the presence of explicit
delays in the negative feedback loop [154, 160, 78, 165, 166], i.e. modifying
Eq. 2 such that :

$\frac{dm}{dt}=f(p(t-\tau))-cm$ (3)

Such delays could be due to the combination of various processes, including
transcription/translation kinetic delays. An even simpler model consists in
making a quasi-static approximation on either $m$ or $p$ to get an effectively
reduced system of the form

$\epsilon\frac{db}{dt}=f(b(t-1))-b$ (4)

This equation is identical to the famous Glass-Mackey oscillator [167, 168],
proposed in a broader physiological context to model phenomena ranging from
dynamical respiratory to hematopoietic diseases. In this equation, the time is
expressed in units of the delay $\tau$. $\epsilon$ quantifies the ratio
between the delay and the time scale of accumulation of $b$. For simplicity,
we will study first this simplest version of the delayed model, but we refer
to [169] and below for more explicit versions of this model in the context of
somitogenesis.

####  1 Conditions for oscillations

We will be especially interested in the limit of small $\epsilon$, where the
delay is much longer than typical degradation kinetics. Simple intuitive
results can be immediately derived. In particular, to understand why there is
an oscillation, it is insightful to take the limit $\epsilon\rightarrow 0$
which gives [78]

$b=f(b(t-1))$ (5)

If $f$ is a (well-chosen) decreasing function of $b$ , one can define two
values $b_{1}$ and $b_{2}$ such that $b_{1}=f(b_{2})$ et $b_{2}=f(b_{1})$. In
such a situation, starting say in $b_{1}$, the system stays in $b_{1}$ for a
time of order $1$, then, since $b=f(b(t-1))$ switches to $b_{2}$, stays there
for a time of order $1$, before switching back to $b_{1}$. This implements a
simple "toggle" oscillator of period $2$, or in the original units, twice the
delay of the oscillator [78]. This yields an immediate prediction: if one can
modify the delay, then one should be able to change the clock period. In fact,
one can easily derive the conditions on the parameters to get oscillations
[170, 168] to confirm the relevance of this limit. In a nutshell, the
repression function $f$ needs to be "strong" enough (e.g. in terms of
cooperativity of the negative feedback), correspondingly $\epsilon$ needs to
be small enough, and the delay $\tau$ big enough to have oscillations. In the
Appendix B3, we derive the full conditions to have oscillations.

####  2 Cycle description

Figure 1: A) Numerical trajectories of the delayed oscillator 4 with step-like
repression $n\rightarrow\infty$ and $\epsilon=0.5$. (B) Reconstruction of the
limit cycle in 2D: $b(t)$ as a function of $b(t-\tau)$ for the purple
trajectory in (A), production rate $r=10^{6}$. The arrows correspond to a
fixed time interval $\Delta t$, indicating a slowing down of the trajectory as
it approaches the origin. (C) Period of the simulated oscillations $T(r)$
compared to the approximate solution 15. For numerical calculation of the
period, the time points where the trajectory derivative changes sign (the
oscillator switches between "on" and "off" and vice versa) were detected.

In the very strong repression limit, taking $f$ as a step function ($f(b)=r$
if $b<1$, $0$ otherwise), the limit cycle can be simply described with piece-
wise linear equations with two different phases: a pure degradation phase when
transcription is off, and a pure production phase when transcription is on.
One can then analytically describe the limit cycle [171]. We get

$\epsilon\dot{b}^{-}=-b^{-}$ (6)

in the degradation phase (as long as $b(t-1)>1$) and

$\epsilon\dot{b}^{+}=r-b^{+}$ (7)

in the production phase (as long as $b(t-1)<1$) where we use $\pm$ superscript
to indicate which phase (degradation vs transcription) the system is.

Let us start at a point in the cycle where $b^{-}$ is crossing $1$ (which
means that $b>1$ for $t<0$). We then get for $t\in[0,1]$

$b^{-}(t)=e^{-t/\epsilon}$ (8)

At $t=1$, $b^{-}(1-1)<1$ so that $f$ flips to $r$. We then get

$b^{+}(t)=r+(b_{1}-r)e^{-(t-1)/\epsilon}$ (9)

with $b_{1}=e^{-1/\epsilon}$. This remains valid until $b^{+}(t_{1}-1)=1$. One
can obtain $t_{1}$ by Taylor expanding Eq. 9 :

$t_{1}=2+\epsilon\ln\left(\frac{r-e^{-1/\epsilon}}{r-1}\right)\sim
2+\epsilon\ln\left(\frac{r}{r-1}\right)$ (10)

at lowest order in $\epsilon$. Then $b$ decreases until $1$, reached at
$t=t_{2}$, closing the cycle. The equation for $b$ between $t_{1}$ et $t_{2}$
is :

$b^{-}(t)=b_{max}e^{-(t-t_{1})/\epsilon}$ (11)

with $b_{max}=b^{+}(t_{1})$. By a Taylor expansion, one finds the value

$t_{2}=t_{1}+\epsilon\ln b_{max}\simeq
2+\epsilon\ln\left(\frac{r^{2}}{r-1}\right)$ (12)

at lowest order in $\epsilon$ (notice that $b_{max}\sim r$ up to an
exponential correction in $e^{-1/\epsilon}$).

It is also useful to extract the duration of the "on" and "off" phases, in
particular

$\displaystyle t_{on}$ $\displaystyle=$
$\displaystyle\tau(t_{1}-1)=1+\epsilon\ln\left(\frac{r}{r-1}\right)$ (13)
$\displaystyle t_{off}$ $\displaystyle=$
$\displaystyle\tau(t_{2}-t_{on})=1+\epsilon\ln r$ (14)

There is a simple intuition here: $t_{off}-1$ essentially is the time to
exponentially decrease due to degradation from $b_{max}\simeq r$ to $1$, while
$t_{on}-1$ is the time to increase from $0$ to $1$ in the production phase.
Notice in particular that $t_{off}$ increases with $r$ (it takes more time to
decrease from a higher value) while $t_{on}$ decreases with $r$ (it takes less
time to reach $1$).

Getting back to the initial time units, one gets for the overall period:

$T=2\tau+\epsilon\tau\ln\left(\frac{r^{2}}{r-1}\right)$ (15)

This equation nicely combines the two-time scales of the system: the
(dominating) delay $\tau$, fixing the period, and $\epsilon\tau$ corresponding
to the transient phases of the system with magnitudes dominated by production
and time scale by degradation.

In [172] a generalized version of this calculation is derived to compute the
perturbations due to multiplicative noise, i.e. an equation of the form

$\epsilon\frac{db}{dt}=f(b(t-1))-b+b\eta$ (16)

where $\eta$ is a standard Gaussian noise
$<\eta(t)\eta(t^{\prime})>=\sigma_{0}^{2}\delta(t-t^{\prime})$, and $f$ the
same step function as above, i.e. $f(b)=r$ if $b>1$, $0$ otherwise. The main
motivation for such calculation comes from experiments, where it has become
recently possible to monitor in real-time concentrations of fluorescent
proteins, which allows for direct comparison between theory and data. It is
argued in [172] that a multiplicative noise is best suited to account for
data. We summarize the steps and results of this calculation in Appendix B4.
It is possible then to go back to experiments and to fit those distributions,
with very good agreement. In particular, in presence of noise $b_{max}$ is the
addition of a noisy production term with $b_{min}$, so we expect it to have
not only a higher average value but also higher variance, which is clearly
seen experimentally. One can also use similar self-consistent approximation to
get distributions of $t_{on},t_{off}$ and numerically estimate the expectation
value and variance of the times for each phase. Compared to the deterministic
case described above, there is a correction in the average due to the noise
autocorrelation.

####  3 Infinite period bifurcation for increased production rate

In the limit of big production rate $r$, we have seen above that $t_{on}\simeq
1$ while $t_{off}\simeq 1+\epsilon\ln(r)$. In particular if the production
rate $r$ is becoming very big, we might get a strong asymmetry between the two
phases of the cycles, with a much longer off phase. It then makes sense to
assume $r$ is very big, to define a corresponding small parameter
$\epsilon_{r}=1/r$, and to introduce the new dimensionless variable
$c=\epsilon_{r}b$. Again there is a degradation and a production phase such
that :

$\epsilon\dot{c}^{-}=-c^{-}$ (17)

in the degradation phase (as long as $c(t-1)>\epsilon_{r}$) and

$\epsilon\dot{c}^{+}=1-c^{+}$ (18)

in the activation phase (as long as $c(t-1)<\epsilon_{r}$). Notice that $c$ is
necessarily bounded between in the interval $[0,1]$.

Assuming the limit cycle starts at $c_{max}\simeq 1$, the off phase clearly
takes a time of order $-\ln\epsilon_{r}$ to reach the threshold value
$\epsilon_{r}$. The on phase has a duration roughly equivalent to the delay,
and if $\epsilon$ is small enough, $c$ has ample time to get back to a value
of order $1$. So we get to a situation where there is a (relatively) short
transcription phase where the system resets from $0$ to a high value
$c_{max}$, followed by a slow decay from $c_{max}$ to $0$. This behavior and
the logarithmic period divergence are very reminiscent of the classical "fire
and integrate" neuronal model in the homoclinic bifurcation regime (see [135]
and Appendix A). Interestingly, in the limit $\epsilon_{r}\rightarrow 0$, we
thus get an infinite period bifurcation with excitability, where it takes an
infinite time to reach $0$ from a positive value, but where an infinitesimal
fluctuation pushes the system to an excursion. See Fig. 1 for a simulation of
this system and a comparison between the simulated period and the analytical
one from Eq. 15

####  4 "Realistic" delayed models and experimental validations

In [160, 169], Monk and co-workers study delayed models in a broader systems
biology context (see also [166]). It is proposed that delayed negative
feedback loops are at the core of several systems, not only Hes1 oscillator in
somitogenesis, but also the oscillations of the oncogene p53 [173] and the
immune NK$\kappa$-B system [174]. A common feature of those systems are
"ultradian" oscillations (i.e. shorter then the 24-hour circadian period),
where one gene represses its own expression via multiple slow steps. As
pointed out in [160], the idea that a gene can repress itself to generate an
oscillation is first due to Goodwin [175], but the so-called Goodwin
oscillator requires very strong, almost step-wise non-linearity to oscillate,
which appears unrealistic. In comparison, a delayed model appears more
plausible experimentally, especially since there are natural biological
delays, e.g. transcription and translation.

It is however important to check if the experimentally measured parameters are
in fact consistent with the regime of oscillations observed and the period
measured. In [160, 169] a two-step model is studied by combining a delayed
equation for mRNA production (i.e. Eq. 3) with a rate equation for production
(Eq. 2). In mouse, half-lives of both hes1 mRNA and proteins are slightly
higher than $20$ minutes, and a reasonable estimate for transcript elongations
(accounting for delays) is another 20 mins. Numerically, this gives an
oscillatory period of roughly 2 hours, with a 20 mins shift between RNA and
protein peak, in excellent agreement with experiments! This is a similar order
of magnitude as the calculation in Eq. 15 with an $\epsilon$ of order $1$
(which assumes fast dynamics on either the protein or the mRNA);
interestingly, it is also found that the oscillation period does not depend
much on the repression threshold but depends much more on the half-life of
either proteins or RNA (period ranging from 100 to 140 mins with either half-
life ranging from 10 to 30 mins; see also [165]).

This leads to a direct experimental prediction: changing the mRNA half-life
should impact the segmentation period. To do this, in [176], mice with $30$
mins half-life for Hes7 are generated. It is then observed that segmentation
is in fact impaired after 3 to 4 cycles. Importantly, Hes7 loss of functions
presents no segment at all, which suggests that mutants with longer half-
period are still functional in the beginning of the segmentation. The authors
interpret this phenotype as a "dampening" of the oscillation: from the
theoretical analysis (e.g. Eq 27), a too-long half-life would indeed lead to
oscillation death via a Hopf bifurcation. Another route to test those models
is to modify delays. In [177], it is observed that introns in the Hes7
sequence could lead to overall delays (introns are part of the RNA that have
to be removed prior to translation). A reporter gene is constructed, with and
without introns, and it is found that while the gene with intron has similar
kinetics to the endogenous Hes7, expression of the gene without introns is
advanced by roughly 20 mins. In fact, it turns out that decreasing the $\sim$
40 mins delay from [176] to 20 mins indeed abolish the oscillations.

While those changes are consistent with a delayed model, they are not true
validation since they both destroy oscillations. A better route to validate
the model is to modify the period without breaking segmentation. This is done
in [178], where the number of introns is simply reduced. There are three Hes7
introns in mice. If only the third intron is conserved, one or two extra
somites form in the cervical region compared to WT (for more caudal regions,
segmentation is impaired and all vertebrae are fused). With the help of a
reporter similar to [177], it is found that the one intron construct turns on
Hes7 5 mins earlier than WT, suggesting shorter delay and subsequent shorter
period (leading to more somites). The oscillation itself can be monitored with
luciferase, and the one-intron clock indeed oscillates with a period of 115
mins vs 126 mins for the wild type. Consistent with this shorter period, the
somite size are reduced by roughly 10 %.

####  5 Experiments: Delay or phase shift ?

The modifications detailed in the previous section are consistent with delayed
models, but it is not clear that other models might not explain experiments
equally well. A recurring issue in the overall discussion is an ambiguity on
the meaning of "delay". Very often, the delays described experimentally
correspond to delays between expression peaks. In other words, they correspond
to phase shifts between different oscillating components. But the existence of
such phase shifts does not imply that the core mechanism is driven by a
delayed differential equation similar to Eq. 3, nor that a delay in expression
corresponds to an explicit delay in a process corresponding to parameter
$\tau$.

Such questions are especially crucial when relating theory to experiments to
understand the origin of the segmentation period [179, 105]. For instance, the
delayed transcription model was tested in zebrafish by considering mutants of
her1 and her7 [180]. Their respective gene sequences have very different
lengths, so if delays come from the time to transcribe genes, one should
observe very different segmentation clock periods in mutants. But their period
is the same, indeed suggesting the "delays" (or rather phase shifts) come from
other processes.

Another issue is that delayed models can "easily" oscillate, so are possibly
less realistic and might miss important properties. For instance, in the
circadian clock context, in [154] it has been shown that the delayed model
bypasses the need for a positive feedback loop, which might in fact be
biologically relevant [181]. Also notice that when a positive feedback loop is
present, the divergence of the flow is not necessarily negative so 2D models
can oscillate. We will see below that segmentation oscillators present
features reminiscent of excitability and relaxation dynamics, which indeed
typically require such positive feedback loops, and can thus give different
bifurcation patterns and oscillatory behaviours. From a practical standpoint,
delayed models are also often more difficult to study mathematically (see e.g.
[163, 182]), because they explicitly introduce tight couplings between
variables expressed at very different (delayed) times.

###  2 Molecular (clock and) Wavefront models

The models in the previous section mostly focused on the modeling of the
oscillation itself. But as clearly discussed in the section Phase models,
segmentation is not only about the clock: it also is about the transition from
an oscillatory system to a stable pattern of gene expression accompanied by
somite formation. So a full model of segmentation should not only focus on the
oscillatory mechanism but also account for such transition, like the original
clock and wavefront or Meinhardt models. In particular, the experimental
realization that FGF drives the dynamics of the wavefront led to the proposal
of different models more explicit from the molecular standpoint on FGF
influence on the oscillation.

####  1 Cycle model with FGF gradient

Following papers showing the clear influence of FGF on the front determination
[81, 10], Baker et al.updated the cell cycle model from [131] to include a
control for the determination front, assumed to be done by an FGF morphogen
gradient (while also being more agnostic on the clock nature). In brief, a new
variable $w$ is added to account for FGF, with the following dynamics :

$\frac{\partial w}{\partial t}=\chi_{w}-\eta
w+D_{w}\frac{\partial^{2}w}{\partial x^{2}}$ (19)

This is a classical production/degradation/diffusion system, with

$\chi_{w}=H(x-x_{n}-c_{n}t)$ (20)

where $x_{n},c_{n}$ are constant and $H$ the Heaviside function similar to the
cell cycle model [131] (Eq.27). We thus see that the production region of FGF
moves as a function of time with speed $c_{n}$. In a stationary regime, we
thus expect to observe a gradient of $w$ moving in the direction of increasing
$x$, which accounts for the experimental dynamics observed in [10].

Then the exact same model as [131] for oscillations coupled to "pulse"
generations is used (Eqs.22-23, with identical forms for $f,g$), but the model
is modified to redefine the activation of genes $u,v$ in a region of lower FGF
(i.e. past the front), namely Eq.27 is replaced by :

$\chi_{u}=H(w^{*}-w)\qquad\chi_{v}=H(t-t_{w}(w^{*},x)-t_{s})$ (21)

The $\chi_{u}$ term is relatively straightforward to interpret: $u$ can get
activated only when $w$ is lower than a constant $w^{*}$, defining the
determination front. The $\chi_{v}$ looks very different from the term in Eq.
27: $t_{w}(w^{*},x)$ is defined at the time when $w$ reach $w_{*}$ at position
$x$, and $t_{s}$ accounts for the period of the segmentation clock. So it is
assumed that the system gets activated at the front, but the $v$ pulse is
delayed for a (segmentation) clock period to activate.

It turns out both those terms are in fact mathematically equivalent to the
terms in Eq. 27, so the model is not fundamentally different but simply ties
parameters of the initial model [131] to the dynamics of $w$. To see this, let
us define $x^{*}(w^{*},t)$ the position at which $w$ reaches $w^{*}$ at time
$t$. Since the gradient of $w$ is moving with speed $c_{n}$ we thus have

$x^{*}(w^{*},t)=x^{*}(w^{*},0)+c_{n}t$ (22)

With this, we can redefine the $\chi_{u},\chi_{v}$ production terms in terms
of $x$ to get

$H(w^{*}-w)=H(x^{*}(w^{*},t)-x)=H(c_{n}t+x^{*}(w^{*},0)-x)$ (23)

so that $\chi_{u}$ is completely equivalent to the initial cell cycle model by
taking $x^{*}(w^{*},0)=x_{1}$ in Eq. 27.

Noticing that from Eq. 22 we have

$c_{n}t_{w}(w^{*},x)=x-x^{*}(w^{*},0)$ (24)

we get

$\displaystyle\chi_{v}$ $\displaystyle=$ $\displaystyle
H(t-t_{w}(w^{*},x)-t_{s})$ (25) $\displaystyle=$ $\displaystyle
H(c_{n}(t-t_{w}(w^{*},x)-t_{s}))$ $\displaystyle=$ $\displaystyle
H(c_{n}t-x+x^{*}(w^{*},0)-c_{n}t_{s})$

we see that we get again a similar equation to Eq. 27 by choosing
$x_{2}=x_{1}-c_{n}t_{s}$. This also makes clear that the $x_{2}$ position
where the pulse of $v$ is activated simply is the position reached by the
front after a period of the clock $t_{s}$.

One interest of the model though is that it offers clear tests of what FGF
concentration changes should entail on the front dynamics since they directly
affect variable $w$. For instance, to model perturbations by a bead soaked
with SU5402 similar to [81], one can assume that $w$ is transiently going down
(due to a sink term for instance). This directly affects the terms $\chi_{u}$
and $\chi_{v}$ which are activated later than they should be, so that one
transient big somite is formed in this model. However, notice that the initial
experiments from [81] are not quantitative and that there likely is a rich
dynamics of FGF inhibition (e.g. diffusion), which could be important to
explain other experimental aspects, such as the formation of several smaller
somites anterior to the bead.

####  2 Somite determination: RA-FGF bistability downstream of coupled
Notch/Wnt oscillator

The model discussed in the previous section is really much more focused on the
front dynamics than on the clock. For instance, the clock is modeled with
pulses of $v$, themselves activated only past the front so that there is no
real account for the multiple oscillatory pathways in the PSM way before the
front.

In a series of papers, Goldbeter and co-workers rather suggested detailed
models of the entire process from oscillation to somite formation [183, 184].
Following descriptions of oscillatory pathways, a first model was built to
couple three oscillators on three different pathways [184]. Three" modules"
are explicitly considered and modeled: a Notch module relying on a negative
feedback loop where Lfng prevents its own activation, a negative feedback loop
in the Wnt pathway implicating Axin2 via Protein-protein interactions, and a
negative feedback loop in the FGF pathway. Interestingly, each of those
modules is proposed to oscillate by themselves, with close but different
periods, but it is hypothesized that they are coupled in the embryo to form a
global oscillator.

Downstream of this global oscillator, Goldbeter and Pourquié [183] proposed
that the bistability of the traditional clock and wavefront model (Eq. 1 and
Fig. 2) is implemented by mutual repression between FGF and Retinoic Acid
(RA). Both proteins are known to form an anteroposterior gradient with
inverted polarity (Fig. 4). Mutual repression between genes indeed is a
classical motif displaying bistability [185, 186]. Just like the classical
Clock and wavefront model, the clock needs to deliver a synchronous signal in
the future presumptive somite to trigger segmentation, and it is proposed that
Mesp2 (which is not explicitly modeled) could also play a role in this
synchronous activation. Remarkably, this model provides a simple explanation
of the time dependency defining the bistability region in Eq. 1 : it simply
arises here from the dynamics of the FGF mRNA and proteins, which degrade as
cells get more anterior [10]. Consistent with the bistability hypothesis, RA
response elements (RARE) sharply turn on at the determination front [90].
However, as said earlier the precise role of RA is unclear since mutants of RA
still form somites [90, 91], but other proteins could play a similar role.

###  3 Somite AP patterning: Inverse problem approach

As said in the introduction of this section, there are multiple levels to
consider to describe the process, and it could be useful to first focus on
possible biochemical processes/networks rather than on specific genes. A
popular route in systems biology is to frame this as an "inverse" problem
question, i.e figuring out (in a relatively unbiased way) what kind of
networks :

  * •

gives rise to oscillations

  * •

stabilizes into patterns of genetic expressions, e.g. similar to c-hairy in
chick or pair-rule genes in insects.

The Meinhardt model offers an example of such a model but presents issues, in
particular, because it heavily relies on diffusion/cell couplings to stabilize
a pattern. So one can wonder if other solutions exist. A popular way to
address such problems is to computationally explore the space of models
performing specific functions [187, 188, 189, 190]

####  1 Clock and switch model

We used simulated evolution of models of gene networks[191, 192], an approach
that led to great success to predict general features in many contexts from
biochemical adaptation [193] to immune recognition [194]. Evolution is
simulated as a standard mutation-selection process on a population of
networks, modeled mathematically using ordinary differential equations. Random
mutations include changes in parameters and changes in network connectivity
(addition/removal of genes in the network and interactions between them), see
[195] for a review of the method and [192] for a recent Python-based
implementation. We used this method to evolve networks giving rise to spatial
patterns under the control of either static or dynamic morphogens [196], in
the absence of cellular communications or diffusion. Importantly, the only
evolutionary pressure (or fitness) was to create patterns of gene expression
of one specific gene $E$ along a one-dimensional axis, so we neither imposed
constraints on the type of genetic architecture we were looking for, nor on
the type of patterns (e.g. on the size or nature of stripes).

Evolution is done under the control of a moving input parameter of the form
$G=H(x-vt)$, in the initial work $G$ was effectively modeling a sliding
gradient, but in the following, we will show results for $H$ being a step
(Heavyside) function for simplicity. A simple network architecture leading to
many domains of expression of $E$ spontaneously evolves Fig. 2. The simulated
evolutionary pathway was highly reproducible from one simulation to the other,
Fig. 2 A, in the sense that successful simulations systematically follow
similar evolutionary trajectories, that we describe below.

The first step in evolution is for a gene $E$ to become bistable via a
positive feedback loop. In terms of equations, this simple network can be
described with :

$\frac{dE}{dt}=\rho\mathrm{max}\left(G,\frac{E^{n}}{E^{n}+E_{0}^{n}}\right)-\delta_{E}E$
(26)

The $\mathrm{max}$ term encodes a transcriptional dynamics akin to a logical
OR, meaning that gene $E$ can be activated either by $G$ or by itself (via a
Hill function). Such logic is necessary to ensure that gene $E$ can be first
activated, then is able to self-sustain when the $G$ morphogen has
disappeared, Fig. 2 A and B Step 1, leading to bistability. In systems
biology, such mechanisms are usually called "switch" and generally rely on
self-activation like here, or on mutual repressions [185, 186], (see also
[197] for a detailed discussion on how to model bistable systems).

Then, the second step is the evolution of a stripe module, where a gene $R$
limits the extension of $E$ Fig. 2 A and B Step 2, giving equations of the
form :

$\displaystyle\frac{dR}{dt}$ $\displaystyle=$
$\displaystyle\rho_{R}G-\delta_{R}R$ (27) $\displaystyle\frac{dE}{dt}$
$\displaystyle=$
$\displaystyle\rho\mathrm{max}\left(G,\frac{E^{n}}{E^{n}+E_{0}^{n}}\right)\frac{1}{1+(R/R_{0})^{m}}-\delta_{E}E$
(28)

This stripe module relies on what is called an incoherent feedforward loop
[198, 199], because the influence of G is positive on both R and E, but then R
"incoherently" represses E. Concretely, the E switch is turned on by G, then
off by R, ensuring a localized expression of E. Such a module is heavily
studied in multiple (evolutionary) contexts, see e.g. [200].

From there, new regulations evolve upstream of $R$, e.g. new repressions.
Remarkably, eventually, evolution finds that a very simple way to generate a
complex pattern of expression of $E$ is to close a negative feedback loop with
delay on $R$, Fig. 2 A and B Step 3, to get a system of Equations of the form
:

$\displaystyle\frac{dR}{dt}$ $\displaystyle=$
$\displaystyle\rho_{R}\frac{G}{1+(R/R_{*}(t-\tau))^{p}}-\delta_{R}R$ (29)
$\displaystyle\frac{dE}{dt}$ $\displaystyle=$
$\displaystyle\rho\mathrm{max}\left(G,\frac{E^{n}}{E^{n}+E_{0}^{n}}\right)\frac{1}{1+(R/R_{0})^{m}}-\delta_{E}E$
(30)

Looking in detail, Eq. 29 is very similar to Eq. 4, and with proper parameters
(selected by evolution) gives rise to oscillations of $R$, as long as $G$ is
high enough. Then $R$ influences $E$ from Eq. 30: when $R$ is low, the $E$
switch is activated by $G$ but when $R$ is high, $E$ is abolished and the
switch is turned off. This means that the system alternates between regions
where $E$ is active or inactive depending on the value of $R$. Lastly, when
$G$ eventually disappears, the final state of the system depends on the
remaining activity of $E$ : if $E$ is high enough, it sustains its own
expression and gives rise to a local gene expression, Fig. 2 C. Because the
time of disappearance of $G$ depends in a linear way on the position in the
embryo, at steady state, one obtains an alternation between regions of
sustained $E$ and regions where $E$ could not maintain itself due to $R$. Thus
the network effectively turns the $R$ temporal oscillations into a striped
pattern of $E$. Interestingly, when those simulations are repeated, one sees
the same kind of evolutionary pathway emerging again and again leading to an
oscillator combined with a bistable system, with only changes in the precise
implementation of the different parts of the network. For instance, in some
cases, we saw the emergence of a Repressilator [201] as the controlling
oscillator.

So the evolved network implements phenomenological dynamics closer to
Meinhardt’s model than the traditional clock and wavefront model: in
particular, just like the Meinhardt model, the clock dies out, and a bistable
component discretizes the phase of the clock into a pattern, that we interpret
as providing AP polarity to somites. But there are key differences. In
Meinhardt’s model, the exact interactions between A and P underlie both
oscillations and patterning, and the front emerges from non-homogeneous
initial conditions and cell-to-cell interactions. Conversely, here this
evolved "clock and switch" model does not require any cell coupling to produce
a pattern. It is a purely local model. Also, the striped pattern emerges from
the interplay of two, independent, genetic modules (first the switch due to
the bistability of E, second the clock due to the negative feedback of R),
under the control of a moving morphogen $G$.

The nature of the wavefront (or primary waves) in such a clock and switch
model is worth discussing. In phase models, there is an unambiguous position
where the frequency of the oscillators goes to $0$, but we already saw that
there is some ambiguity on the front definition for more explicit models (like
Meinhardt’s), which we face again here. In the clock and switch model, there
is a spatially extended transition zone, from the oscillatory behavior to the
bistable behavior, controlled by $G$. In fact, there are at least two primary
waves. It is because the generic way to go from an oscillatory to a bistable
expression is through two bifurcations, as illustrated in Fig 2 A, and those
will typically not happen for the same values of $G$ without symmetries in the
equations.

As $G$ decreases, the system first goes through a Hopf bifurcation, killing
the oscillation, then through a saddle-node bifurcation, creating a bistable
system. In between those two bifurcations, there is a region of $G$ where
there is only one single fixed point for $E$. If the system spends too much
time in this region, any level of molecular noise would jiggle $E$ around so
that $E$ would not be able to remember the phase of the clock (position on the
cycle). So, for the bistable system to keep information about the phase of the
oscillation, there are only two possibilities: either noise should be strongly
suppressed in this region (e.g. with the help of some extra variable
reinforcing the phase encoding), or the transition should happen rapidly (e.g.
$G$ presents a step-like behavior, effectively collapsing the two primary
waves into one, as simulated here). The first scenario appears more consistent
a priori with data since there is a clear gradient in the PSM of various
observables (from morphogens to frequencies), suggestive of a rather gradual
process. Notice the two bifurcations correspond to the two modules that are
evolved: the clock needs to die for the bistable system to turn on, and the
morphogen should adjust those transitions.

Figure 2: (A) Evolved network, corresponding to Eqs. 29-30. $R$ represses $E$
and self-represses with delay. $E$ self-activates. Circle numbers indicate the
order of appearance of each interaction during evolution (see corresponding
kymographs in B). We show below the network a scheme of the associated
bifurcation diagram: when activator $G$ decreases, the system goes through a
Hopf bifurcation (H) and then a saddle-node bifurcation, i.e. transitioning
from a clock to a bistable system. Notice the modular structure of the network
recapitulating both the bifurcation diagram and the evolutionary steps. (B)
Kymograph for $E$ under control of a moving, step-wise $G$ dynamics, for the 3
evolutionary steps corresponding to the addition of three interactions
indicated in A (C) Two single-cell trajectories with different dynamics of $G$
for the final evolved network. Depending on the timing of $G$ disappearance,
the system ends up either in a high or low $E$ state

####  2 Inferred PORD and hybrid models

For small enough networks, another way to figure out network topologies
performing a specific function simply is a direct enumeration, a method
pioneered in [189]. The advantage of this method is that it provides a more
extensive sense of simple possible mechanisms, with the drawback that
solutions might not be "evolvable". Also notice that by design such an
approach can not explore networks with more than a couple of nodes because of
combinatorial explosion.

When this approach is applied to stripe formation in a somitogenesis-like
context, Cotterell et al.[122] first recover the clock and switch mechanism
described in the previous section (that they call "Clock and Gradient Model")
but they also exhibit another mechanism that they call "progressive,
oscillatory reaction diffusion" (PORD) model. Equations for this mechanism are

$\displaystyle\frac{dA}{dt}$ $\displaystyle=$ $\displaystyle
ReLu_{s}(k_{1}A-k_{2}R+F+\beta)-\mu A$ $\displaystyle\frac{dR}{dt}$
$\displaystyle=$ $\displaystyle\frac{k_{3}A}{1+k_{3}A}-\mu R+D\Delta R$

where $ReLu_{s}(x)=ReLu(\frac{x}{1+x})$ and $ReLu(x)=max(x,0)$ is the
rectified linear function (commonly used in machine learning). $F(x,t)$
depends on both time and space and is assumed to be a sliding morphogen
gradient, encoding FGF. $A$ and $R$ are respective activator and repressor.
$R$ is clearly activated by $A$ only. The $RELU_{s}$ function for $A$ thus
encodes a saturating activating function, positively influenced by $A$ (and
$F$) and negatively influenced by $R$. So the network dynamics with $A$ self-
activating and $R$ diffusing and repressing $A$ is very similar to the reduced
version of Meinhardt model Eqs. 16-17, see phase space and kymographs for a
system with constant F in Fig. 3, to be compared with Figs. 9, 10. It should
however be mentioned that this PORD model is more realistic than the Meinhardt
model in the sense that it relies on more explicit transcriptional
biochemistry, in particular with fewer symmetries in the equations.

Otherwise, the phenomenology of this model is largely identical to the
Meinhardt model: in absence of diffusion, the system displays relaxation
oscillations, with A switching rapidly between $0$ and its maximum value, and
$R$ more smoothly oscillating. When diffusion is included, a pattern emerges
from short-range interactions, where bistable stripes of $A$ form, and
activate $R$ which diffuses and stabilizes the pattern. This requires no
external morphogen gradient, even though a modulation by morphogen gradient
$F$ is included (an external control could also be added to the Meinhardt
model to regulate other aspects such as stripe size, see Eq. 21). Again, as we
discussed in the section on the Meinhardt model, similar behaviors can be
obtained with various relaxation oscillators.

Figure 3:  PORD model with constant $F$. On the left, the phase portrait is
represented, with $A$ nullclines in blue and $R$ nullclines in green, and the
limit cycle is in red. On the right, we show kymographs, temporal behavior in
one posterior cell, and final pattern. The behavior of this system is
qualitatively similar to the Meinhardt-VanderPol model see Figs. 9, 10

There are several possible experimental evidence for the FGF influence. For
instance, in this model $F$ can lead to local self-sustained activation of
$A$. This appears consistent with the effect of local cutting in the embryo,
which induces a FGF response and a stripe of Lfng. Because of the way the
pattern is formed (similar to the "domino" effect mentioned for the Meinhardt
model), local changes of FGF might "propagate" and influence pattern formation
very far from the initial perturbation. This prediction is validated by the
fact that a cut in the embryo gives a diagonal, striped pattern far from the
cut.

It is, however, not entirely clear how the imposed motion of the initial
gradient of $F$ in [122] couples to the wavefront, since both those speeds are
governed by different, uncoupled parameters ($v$ is imposed externally, and
the wavefront progression is a purely local variable, as detailed in [122]).
Pantoja-Hernandez et al.[202] include a more explicit mechanism for the
coupling between a moving morphogen and the relaxation-diffusion dynamics.
They consider a model of the form

$\displaystyle\frac{dA}{dt}$ $\displaystyle=$ $\displaystyle
P_{A}(A,R,\beta)-\mu_{A}A+D_{A}\Delta A$ (31) $\displaystyle\frac{dR}{dt}$
$\displaystyle=$ $\displaystyle P_{R}(A)-\mu_{R}R+D_{R}\Delta R$ (32)

where

$\displaystyle P_{A}(A,R,\beta)$ $\displaystyle=$
$\displaystyle\frac{\beta+(A/K_{1})^{n_{1}}}{1+(A/K_{1})^{n_{1}}+(R/K_{2})^{n_{2}}}$
(33) $\displaystyle P_{R}(A)$ $\displaystyle=$
$\displaystyle\frac{(A/K_{3})^{n_{3}}}{1+(A/K_{3})^{n_{3}}}$ (34)

which again gives rise to relaxation oscillations. They show that for large
values of $\beta$ compatible with oscillations, alternating spatial patterns
do not appear, while they can form spontaneously for smaller values of
$\beta$. They thus suggest that a moving wavefront of $\beta$ might control
somite formation: for higher $\beta$, no pattern forms, but when $\beta$ is
graded spatially and reaches smaller values, the Meinhardt-PORD mechanism
kicks in and a pattern stabilizes. So one gets here an external control of the
transition from a pure oscillatory system to pattern formation, reconciling
gradient-induced transitions with pattern formation, into a "hybrid model". It
is suggested that this provides a more robust way to control pattern
formation.

####  3 Experimental evidence for a switch

Such systems-level approach yields important predictions that can be tied to
specific gene dynamics. AP identities within somites are plausible candidates
for the two steady states found in the Clock and Switch model, and should be
associated with specific genes. Notice that anterior and posterior fates
within somites have mutually exclusive markers, so that the evolved positive
feedback loop on $E$ (Eq. 30) could be a mutually repressing system [186,
185]. Experimental evidence is consistent with the idea that anterior-
posterior somite fates are bistable identities discretizing the clock phase.
Genes of the Notch signaling pathway are indeed implicated both in the clock
and in anteroposterior definition, e.g. in zebrafish, deltaC oscillates and
becomes restricted to half a somite, consistent with the behavior of gene $E$
[107]. Dll1 in mouse appears to be necessary and sufficient for the definition
of posterior somite fate [83], and is maintained after somite boundary
formation. It should nevertheless be pointed out that in the clock and switch
model, $E$ does not necessarily have to be expressed in the entire PSM: it is
only required that it oscillates right before the clock stops to encode the
pattern. Mesp2 precisely displays such behavior: it is controlled by the
clock, expressed right at the front, and localizes in the anterior part of the
somite [82] (see also [85] for a study combining theory and experience showing
how a spatial Mesp2 gradient can be established within one future somite in
response to the clock). Interestingly, anteroposterior markers are expressed
stably before boundary formation, and there is multiple evidence that they
control genes implicated in segmental border formation [203]. If somite cells
are indeed bistable, one expects that those fates could be stably induced
independently of the clock, which has been realized in an in vitro system for
Mesp2 and Uncx4.1 [23].

Lastly, evidence from some mutants is also very consistent with the clock and
switch model. Convergence-extension mutants in zebrafish have two-cell wide
somites, and, quite spectacularly, display single-cell wide rows, expressing
anterior or posterior markers in alternation [87]. Furthermore, those
extremely narrow somites nevertheless segregate. This excludes any predefined
length scale in the process (e.g. due to diffusion) that would define somite
length. This is rather consistent with a cell-autonomous, local mechanism,
where the final pattern is more dependent on the time course of a morphogen
like $G$ which, here, would lay down alternating fates. It is also quite
spectacular that only two fates are expressed in such mutants, showing that no
stably expressed "third" state is necessary for somite formation, and that the
clock information is clearly discretized, again completely consistent with the
clock and switch model.

###  4 Landscape geometry of segment formation

Figure 4: (A) Building a geometric model: a system continuously translates
from a dynamical, oscillatory state, towards a multistable one, under the
control of a sliding morphogen $G$. (B) If the transition is non-linear, the
system goes through a Hopf then a saddle-node/pitchfork bifurcation.
Corresponding flows in the phase plane are represented for lower values of $G$
from right to left. The limit cycle is shown in light gray in the 2D plane,
stable fixed points in green, and unstable fixed points in red. A kymograph
for $x$ is shown on the right, with the "blurred" transition zone highlighted.
(C) If the transition is more linear, the oscillation disappears through a
SNIC bifurcation, where new fixed points appear on the cycle (giving rise to
excitability). Notice the sharper kymograph, we also highlight the asymmetry
in the wave.

More and more data are accumulated on embryonic development (e.g. via single
cell RNAseq [204]), and low dimensional representations of data using
unsupervised machine learning techniques (such as UMAP [205]) are becoming
increasingly popular. Mathematical constraints on low-dimensional systems
aiming at modeling development are under intense scrutiny. Models are
generally of the so-called ’Morse-Smale’ type and can be fully characterized
by the way attractors (corresponding to cellular states) appear or disappear
[34]. One can also characterize the way morphogens can modulate cell fates,
and derive experimental predictions from small, "gene-free" models. Recent
examples of this approach include predictions on C. Elegans vulva formation,
and neural progenitor differentiation [206, 207].

This is very reminiscent of the initial clock and wavefront model, based on a
landscape description (e.g. Fig. 2). Assuming that in each cell there is
indeed a well-defined oscillator in the posterior and that oscillations
eventually stabilize into different fates, we expect from dynamical systems
theory that there is only a limited number of possible effective models
describing the cells’ behaviors, and associated bifurcations (in the dynamical
system sense).

Concretely, how can such bifurcations occur? Developmental systems are
regulated by different enhancers at a different time of development. Many
dynamical features of the development are correlated to changes in enhancer
controls, coordinating the expression of genetic ’modules’ [208, 209]. This
has been clearly shown experimentally for gap genes control in fly (where
different enhancers correspond to different initial and final positions of gap
gene Kni [210]), and led to the suggestions that in other insects, there could
be both "dynamic" and "static" enhancers, the former responsible for the
initiation of the patterns and the later for its stabilization. A model of
Tribolium segmentation based on those ideas recapitulates well both the normal
developmental transition and mutants [211].

Elaborating on those ideas, we studied the transition of a system from a set
of enhancers controlling an oscillation toward another set of enhancers
controlling a bistable system [212], Fig. 4. This can be done with explicit,
gene-network models, or with purely geometric descriptions of a 2D flow.
Interestingly, the geometric description captures well all the interesting
features of more explicit gene network models, as we now describe.

In detail, let us consider two 2D flows, one corresponding to an oscillator
$D(X)$ and the other corresponding to a multistable system $S(X)$, where $X$
is a 2D vector. We assume that the dynamics of the system are given by the
following equation :

$\dot{X}=\Theta_{d}(G)D(X)+\Theta_{s}(G)S(X)$ (35)

where the weights $\Theta_{d}$ and $\Theta_{s}$ are functions of an external
parameter $G$. For instance, if we take $\Theta_{d}(1)=\Theta_{s}(0)=1$ and
$\Theta_{d}(0)=\Theta_{s}(1)=0$, then system $X$ moves from the oscillator $D$
to the multistable system $S$ as $G$ decreases from $1$ to $0$. Corresponding
attractors are illustrated in Fig. 4. Practically, $G$ will be again assumed
to be a morphogen traveling over the embryo, so will take the form
$H(x,t)=H(x-t)$ where $H$ is a decreasing function from $1$ to $0$. Other
enhancers can be of course added in a similar way.

We used various polynomials for $\Theta$s and standard models for both the
oscillator and the multistable system. In particular, for the oscillator
$D(X)$ we can use a Poincaré model best written in polar coordinates (Fig. 4 B
and D, left-most flow plot):

$\displaystyle\dot{r}$ $\displaystyle=$ $\displaystyle r(1-r)$ (36)
$\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle 1$ (37)

while possible equations for the bistable system is (Fig. 4 B and D, right-
most flow plot)

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle x(1-x^{2})$ (38)
$\displaystyle\dot{y}$ $\displaystyle=$ $\displaystyle-y$ (39)

ensuring that $(-1,0)$ and $(1,0)$ are stable fixed points.

Let us immediately point out that flows here are very symmetrical, but the
symmetry can be broken without losing the general properties observed [212]
(see also discussion in section From geometric back to phase models ?). We
also checked that the bifurcation diagrams observed and the general properties
of the system are largely preserved when using different attractors (e.g.
different types of non-linear oscillators).

Such a system recapitulates the sequence of two primary waves (Hopf+saddle-
nodes) first obtained through computational evolution, by taking non-linear
functions $\Theta$ and further stabilizing the origin (Fig. 4 B). However, we
also obtained a much more generic behavior, when the transition from the
oscillatory attractor to the bistable one is linear (e.g. with
$\Theta_{d}(G)=1-\Theta_{s}(G)=\theta_{G}$)

$\displaystyle\dot{x}$ $\displaystyle=$
$\displaystyle\theta_{G}(1-r)x-\theta_{G}y+(1-\theta_{G})x(1-x^{2})$ (40)
$\displaystyle\dot{y}$ $\displaystyle=$
$\displaystyle\theta_{g}(1-r)y+\theta_{G}x-(1-\theta_{G})y$ (41)

with $r=\sqrt{x^{2}+y^{2}}$.

It is useful to start with an intuitive explanation of what happens then. If
the $\Theta$s vary smoothly with $G$, e.g. $\Theta_{G}=G$, the dynamics of the
system around a critical value $G_{c}~{}1/2$ is a "perfect" mix of an
oscillator and a multistable system. Right below $G_{c}$, the system does not
oscillate, but by continuity, half of the dynamical flow (corresponding to the
limit cycle) still is present in the dynamics, meaning that for some initial
conditions, the system might still follow part of the cycle. This is
reminiscent of what is observed in excitable systems which have a single fixed
point/steady state, but when perturbed suitably, come back to it through long
excursions very similar to cycles (see Appendix). Such systems are very common
and have been well-studied in theoretical neuroscience [135].

With this geometric modeling, we in fact obtained excitable so-called type I
oscillators [135], which arise through a bifurcation called a SNIC (Saddle-
Node on Invariant Cycle, see Appendix). SNIC bifurcations exactly correspond
to the intuitive picture of a system being at the same time oscillating and
stable: they arise when a fixed point appears on a limit cycle or, conversely,
when a pair of stable/unstable fixed points cancel out on a closed trajectory
to give rise to oscillations (explaining the name of the bifurcation). Right
after a SNIC, the system is normally monostable, but in our case, if both $D$
and $S$ are symmetrical enough, a saddle-node bifurcation will generate a
second fixed point, giving rise to bistability encoding fates, Fig. 4,
Geometric Model 2. So similar to the Clock and Switch model there still are
two primary waves, but with a SNIC bifurcation instead of a Hopf one to
destroy the oscillator.

By varying the functions $\Theta(G)$, one can easily compare different
bifurcation scenarios with the same "boundary conditions", i.e. dynamic and
static flows $S,D$. Kymographs are shown in the right column of Fig. 4B-C for
two scenarios. One can then compare a Hopf+saddle-node scenario to a
SNIC+saddle-node one. In short, the SNIC bifurcation scenario is in general
more robust to all kinds of perturbations (noise, shape of traveling morphogen
$G$) than the Hopf scenario (as quantified in [212]), suggesting that it might
be evolutionarily adaptive. The main reason is that at the Hopf bifurcation,
the dynamics tend to "collapse" at the fixed point so that the phase
information from the cycle can be easily lost through perturbations (as is
visible in the highlighted blurry phase in the kymograph for the Hopf model in
Fig. 4 B). Conversely, in the SNIC scenario, since fixed points appear on or
close to the cycle, information about the phase is more directly encoded and
less easy to perturb (as can be grasped from the overall sharper kymograph in
Fig. 4 C for the SNIC scenario).

Strikingly, the behavior of the system close to the SNIC bifurcation is
qualitatively reminiscent of what is observed in actual cells. As we get
closer and closer to the SNIC bifurcation in this model, the period of the
cycle diverges to infinity (explaining why it is sometimes called Saddle Node
with Infinite Period bifurcation (SNIPER) [213]). So SNIC/SNIPER bifurcations
induced by the transition to bistability offer a very natural explanation for
the considerable period slowing down and the associated kinematic waves
observed in embryos (in particular zebrafish and snake). Another interesting
aspect is that close to the bifurcation, the SNIC system behaves like a
relaxation oscillator, which gives rise to rather asymmetrical wave profiles,
looking more like a sawtooth profile (see e.g. the highlighted "‘teeth" in the
yellow phase of the model 2 kymograph in Fig. 4 C). Such asymmetrical waves
are observed experimentally [97] and motivated the phenomenological coupled
phase-amplitude model (which assumed a form of excitability) described in
Section Phase-amplitude coupling and excitability for oscillation arrest. This
further provides a possible answer to the somite polarity conundrum, like
explained at the end of section Phase models : because of those asymmetrical