# Sequential pattern formation governed by

signaling gradients

###### Abstract

Rhythmic and sequential segmentation of the embryonic body plan is a vital developmental patterning process in all vertebrate species. However, a theoretical framework capturing the emergence of dynamic patterns of gene expression from the interplay of cell oscillations with tissue elongation and shortening and with signaling gradients, is still missing. Here we show that a set of coupled genetic oscillators in an elongating tissue that is regulated by diffusing and advected signaling molecules can account for segmentation as a self-organized patterning process. This system can form a finite number of segments and the dynamics of segmentation and the total number of segments formed depend strongly on kinetic parameters describing tissue elongation and signaling molecules. The model accounts for existing experimental perturbations to signaling gradients, and makes testable predictions about novel perturbations. The variety of different patterns formed in our model can account for the variability of segmentation between different animal species.

###### pacs:

82.39.Rt, 87.17.Pq, 87.18.Hf, 89.75.KdThis is a preprint of the corresponding article published as

Physical Biology 13, 05LT03 (2016).

Morphogenesis, the formation of shapes and patterns in the developing embryo, relies on the tight coordination of cellular actions [1]. During embryonic development, spatial profiles of signaling activity in tissues control the behavior of cells such as proliferation, migration, and differentiation [2, 3, 4]. In vertebrates, a vital morphogenetic process is the segmentation of the elongating body axis during which the precursors of the vertebrae are formed [5]. Segments form rhythmically and sequentially from an unsegmented progenitor tissue, the presomitic mesoderm (PSM), see Fig. 1A. During segmentation, the body axis elongates while the PSM continuously changes its length and eventually shortens. After a species-dependent number of segments is formed, the process terminates when the PSM becomes very small [6, 7]. The total segment number can vary from about ten in frog to several hundreds in snake [8]. The temporal progress of segmentation is controlled by oscillations of the cellular concentration levels of functional proteins in the PSM [9, 5]. These cellular oscillations are achieved through autoregulation of so-called ‘cyclic genes’ [10, 11, 12]. On the tissue level, these oscillations give rise to nonlinear waves that propagate through the PSM [13, 14, 15, 16, 17, 18, 19, 20], see Fig. 1B. A new segment is formed with each completed oscillation at the anterior end of the PSM, corresponding to an arriving wave [20]. Hence, in contrast to pattern formation via instabilities of homogeneous states, segmentation is characterized by a spatially inhomogeneous system with the PSM driving the patterning process. In addition, pattern formation takes place in a dynamic medium: the body axis continuously elongates during segmentation while the PSM at the tail of the body axis shortens until segmentation terminates. A theoretical model that integrates pattern formation, tissue elongation and shortening, and termination of segmentation by a self-organized mechanism is still missing.

To yield robust morphological results, this complex patterning process requires tight integration of spatial and temporal cues. What regulates the integration of tissue growth and patterns of oscillating gene expression in vivo? The elongating body axis exhibits spatial concentration profiles of several signaling molecules. In particular, the PSM displays posterior protein concentration profiles of FGF (fibroblast growth factor) and Wnt and an opposing anterior profile of Retinoic Acid [21, 8, 22, 23, 24, 25, 26, 27, 28, 29], see Fig. 1C. These signaling molecules are thought to be involved in body axis elongation [30] and regulating oscillations [31] and cell fate during segmentation [32], i.e., maintaining cells in an oscillatory state within the PSM and triggering segment formation upon arrival at its anterior end. Moreover, Retinoic Acid and FGF have been reported to display mutual inhibition and to antagonistically regulate genes that are involved in specifying segmented and unsegmented tissue [22, 23]. Basic principles of segmentation were captured by simplified models in which the PSM length is kept constant [16, 17, 33]. These studies use discrete or continuous phase descriptions of the cellular oscillators to describe the effects of cellular interactions on tissue level. Furthermore, the role of signaling activity within the PSM was studied using different approaches including reaction-diffusion models of signaling gradients without coupling to oscillators or coupled to static clocks [34, 35, 36], models of signaling dynamics that is explicitly time-dependent [37], models based on diffusible cyclic gene products [38] and detailed cell-based models [39, 40, 41]. However, whether the interplay of signaling activity and genetic oscillations alone can lead to self-organized gene expression waves and robust segmentation of the body axis in a dynamic tissue is not understood.

In this paper, we present a theoretical description of vertebrate segmentation, in which spatial concentration profiles of signaling activity regulate both growth and the frequency of cellular oscillators that control segmentation. We show that local interaction rules can lead to robust self-organization of genetic oscillations, tissue elongation, and PSM shortening. The resulting patterning system yields the correct spatial morphology and shortening of the segmenting tissue. We introduce a one-dimensional continuum description of the dynamic tissue based on phase oscillators and two signaling activities and varying in space and time, which represent Wnt/FGF and Retinoic Acid activity, respectively. The signaling activities are effective tissue-level representations of the opposing and antagonizing signaling gradients found in vivo, see Fig. 1C. The interactions of the signaling system lead to termination of segmentation after a finite number of segments. We study how the key features of segmentation depend on the kinetic parameters of the signaling system. In a second step, using a simplified scenario with a single signaling gradient and time-periodic wave patterns, we derive analytical relations that explicitly show how the characteristic length scales of segments, waves, and tissue extension arise from our model.

We introduce a curved coordinate axis along the embryo, in which corresponds to the posterior tip of the PSM, see Fig. 1A. We consider a set of cellular oscillators in the PSM, described by their phase in the oscillation cycle. Growth of the tissue and frequency of the oscillators is regulated by a signal that moves with the cell flow and is degraded, see Fig. 2A (black and gray arrows). A second signal that emanates from the formed segments and quickly diffuses triggers additional degradation of . The spatio-temporal distributions of the signals are described by the one-dimensional activity fields and . Furthermore, we introduce a phase field that represents the local state of the cellular oscillators in the PSM [16, 42]. The dynamic equations for the phase field and the concentrations and are given by

(1) | |||||

(2) | |||||

(3) |

where is a velocity field accounting for cell flow, is the intrinsic frequency of the cellular oscillators, is the coupling strength. and are the diffusion constants of and , respectively, and are their basal production rates, is the production rate of in the formed segments, and are decay rates. The rates and indicate the degree of mutual degradation of and , motivated by the antagonistic action of opposing gradients in vivo. The dependence of the production rate of on leads to additional localized production in the center of the formed segments. For simplicity, we here consider a constant length of the source region of , . In our model, the local elongation rate as well as the frequency and coupling strength profiles and are controlled by through the relations

(4) | |||||

(5) | |||||

(6) |

where , , and are a characteristic elongation rate, oscillation frequency, and coupling strength, respectively. The position of the anterior end of the PSM is defined as the point where the level of reaches the threshold level ,

(7) |

see Fig. 2B. We consider open boundary conditions for the phase field, . For the signaling activities, we consider no-flux boundary conditions at the posterior tip, .

The model proposed here shows that the dynamics of segmentation and length decrease of the PSM can arise from local cellular interactions as a self-organized process. To illustrate this, we start the system with a steady-state initial condition in which and form opposing gradients in the absence of phase dynamics, see Fig. 3A. The profile of generates a spatial profile of intrinsic frequencies through Eq. (5). In the source region, the profile of is typically flat due to the balance of production, decay, and growth. The frequency profile leads to a wave pattern, i.e., different parts of the PSM are out of phase [14, 16, 42]. The pattern is flat in the posterior and displays a characteristic wavelength in the anterior, in accordance with experiments [20]. This corresponds to no segments and a PSM of finite length with a flat phase profile. Figs. 3A,B show snapshots of a numerical solution of Eqs. (1–3) with (4–6) for different time points (see also Supplemental Movie 1). As soon as waves leaves the PSM, i.e., as they enter the region where , they become sources of according to Eq. (3), see Fig. 3A. The thus elevated levels of diffuse into the PSM and lead to an increased degradation of in the vicinity of the anterior end (at ) and the profile of shortens towards the posterior, see Figs. 2C and 3A. Consequently, the PSM shortens, see Fig. 3C. Finally, the segmentation process terminates with the PSM reaching zero size, .

To discuss the dynamic features of this model, we define the time-dependent number of waves , the segment length , and the number of formed segments [42],

(8) | |||

(9) | |||

(10) |

The number of waves is the total phase difference between posterior tip and anterior end , and the segment length is given by the local wavelength of the pattern at the anterior end. Fig. 3C shows these functions together with the time evolution of the PSM length. The time dependence of these quantities capture key features of the segmentation process as seen in experiments, such as the decrease in PSM length, the formation time of segments, and the decrease of segment length over time [20, 42]. The non-monotonic time dependence of the number of waves is not observed in experiments [20] and is related to the constant source length, which we have chosen here for simplicity.

We now show that we can obtain a variety of different oscillation patterns and morphologies by changing the dynamics of signaling. As an example, we vary two key parameters: the diffusion constant of the signaling component , and the elongation rate , which sets the scale of cell flow velocity. We define the total time of segmentation as the time at which the PSM has shortened to zero. The average number of waves in the PSM is denoted by , the minimal segment length is , and the total segment number is . Fig. 4 shows these observables as a function of and . We find that the total time of segmentation diverges at the boundary shown in Fig. 4A. Within a parameter range (gray), segments are generated in a time-periodic manner without end (‘infinite snake’) corresponding to the simplified theory Eqs. (11,12). Biologically relevant parameters are found in the green region, where the total time of segmentation is finite. This occurs if advection described by is small enough that can diffuse sufficiently far into the PSM to degrade . In this parameter region, the number of waves decreases with stronger advection, while the segment length increases, see Figs. 4B,C. The latter trend can be understood through the simplified model, see Eq. (18) and below. Interestingly, the total segment number displays a non-monotonic behavior as a function of the elongation rate with the smallest number of segments formed for intermediate values of , see Fig. 4D.

The model with two signaling gradients described by Eqs. (1–6) can generate self-organized segmentation. However, additional insights in the length scales of segments, tissue, and wave pattern arise from the model dynamics can be obtained from the study of a simplified model. In this reduced version of our model, only the posterior signaling gradient is present, see Fig. 2A (black arrows only). This reduced model cannot account for PSM shortening but captures the self-organization of genetic oscillations in wave patterns and tissue elongation. The dynamic equations for and are then given by

(11) | |||||

(12) |

The system of Eqs. (11) and (12) together with Eqs. (4–6) has a solution with a stationary activity profile and a time-periodic wave pattern , where is the collective frequency with which the wave pattern corresponding to the phase profile repeats. Note that the dependence of on makes Eq. (12) a nonlinear equation. The number of waves in the PSM and the segment length at formation are now time-independent and given by

(13) | |||||

(14) |

Interestingly, general properties can be discussed independent of knowledge of the full solution. We here consider the case in which coupling only provides a minor correction to the phase pattern. As an approximation, we set in Eq. (11). For the phase profile , we thus obtain [42]

(15) |

The condition determines the collective frequency as

(16) |

where . Consequently, the segment length is given by , where . The velocity at the anterior end and the posterior concentration can be obtained from the stationary profile , which, according to Eq. (12), satisfies

(17) |

Eliminating in Eq. (17) by using Eq. (4), we obtain an equation for the elongation rate . Integrating from to yields for the case , where . The posterior concentration is determined by Eqs. (17) and (4) at the posterior boundary as . Hence, the segment length is given by

(18) |

where . Since and , the segment length increases with increasing elongation rate. Note that is proportional to the length of the source region, which is the only length scale in the system.

The full solution for the steady state of is given by

(19) |

where and is the principal branch of the Lambert function, defined by the relation [43]. The growth field generated by through Eq. (4) corresponds to a velocity field, where cells move anteriorly and reach their maximum speed at the anterior end of the tissue. From Eq. (19), the PSM length can be obtained using the definition Eq. (7),

(20) |

The number of waves can be obtained using Eqs. (15) and (5) and approximating the velocity field by the velocity at the anterior end, ,

(21) |

Dynamic solutions of the minimal model with constant PSM length are shown in Supplemental Movie 2. The stationary profiles for and describe the periodic wave-like pattern of gene expression that is determined by a signaling gradient. These solutions are similar to the stationary ‘infinite snake’ patterns found in the full model, see Fig. 4A. The waves propagate in a dynamic tissue that expands with a velocity profile determined by Eq. (4).

In this paper, we have shown how segmentation of the vertebrate body plan can arise as a self-organized patterning process controlled by signaling activity coupled to genetic oscillators. In our full model, local rules representing cellular interactions lead to (i) self-organized wave patterns, (ii) tissue elongation, (iii) PSM shortening, and (iv) self-organized termination of segmentation with a finite number of segments. To complement our study, we used a reduced model with one signaling gradient to derive explicit relations between the time and length scales of the wave pattern and segmentation and the biochemical properties of signaling. Our work shows that by varying biochemical parameters of the signaling gradients, a variety of different patterns and morphologies can be generated. Such variations could correspond to the differences in the segmentation process and the resulting morphologies between different species, e.g., in fish, mouse, chick, snake, and frog. Among these species, the observed number of segments ranges from about ten to several hundred, the number of waves ranges from one to five [8]. In our model, the length of the signaling source region sets the length scales of wave patterns and segment length. This correspondence could naturally account for ‘scaling’ of segment size with the size of the organism if this source region occupies a characteristic proportion of the segmenting tissue. Note that the principle of self-organization proposed here does not require diffusion of the posterior signaling molecule in the tissue. Indeed, our work shows that frequency and growth profiles could emerge either from effective diffusion or from advection or from a combination of both [44, 45]. It is worth noting that our model also exhibits two recently discovered wave effects in embryonic segmentation: a Doppler effect that arises from the anterior end of the PSM moving into the waves and a ‘dynamic wavelength effect’ that denotes the decrease of the wavelength at a fixed position over time [20] (both effects can be inspected in Supplementary Movie 1).

There are several possibilities to test whether the mechanism proposed here actually underlies the segmentation process in vivo. One of these possibilities is transient up- or downregulation of the signaling molecules present in the PSM. Experimentally controlled transient reduction of Wnt signaling, for instance, leads to a shorter average PSM length, faster PSM shortening, and longer segments [29], all observations consistent with our model, see Fig. 5A (red curves). Furthermore, our model could be tested by experimentally induced transient up- or downregulation of Retinoic Acid in a quantitative and dynamic assay system, see Fig. 5B. The details of cross-regulation between signaling molecules are yet largely unknown and how, for instance, signaling activity regulates the frequency of the cellular oscillators and cell fate on a molecular level remains an open question. Our model, which is based on available experimental knowledge, yields testable predictions on tissue level that could shed light on the nature of these molecular interactions and motivate further experimental and theoretical research.

## Acknowledgments

We thank Rachna Narayanan and Luis Morelli for fruitful discussions.
This work was supported by the Francis Crick Institute which receives its core funding from Cancer Research UK, the UK Medical Research Council, and the Wellcome Trust. In addition, AO was supported by a Senior Research Fellowship from the Wellcome Trust (WT098025MA) and the Medical Research Council (MC_UP_1202/3).

## References

- [1] Wolpert L, Smith J, Jessell T, Lawrence P, Robertson E and Meyerowitz E 2006 Principles of Development 3rd ed (Oxford University Press)
- [2] Turing A M 1952 Phil. Trans. R. Soc. B 237 37–72
- [3] Rogers K W and Schier A F 2011 Annu. Rev. Cell Dev. Biol. 27 377–407
- [4] Werner S, Stückemann T, Beirán Amigo M, Rink J C, Jülicher F and Friedrich B M 2015 Phys. Rev. Lett. 114 138101
- [5] Oates A C, Morelli L G and Ares S 2012 Development 139 625–639
- [6] Bellairs R, Ede D A and Lash J W (eds) 1986 Somites in Developing Embryos (Nato Science Series A vol 118) (Springer US)
- [7] Tenin G, Wright D, Ferjentsik Z, Bone R, McGrew M J and Maroto M 2010 BMC Dev. Biol. 10 24
- [8] Gomez C, Özbudak E M, Wunderlich J, Baumann D, Lewis J and Pourquié O 2008 Nature 454 335–9
- [9] Cooke J and Zeeman E C 1976 J. Theor. Biol. 58 455–476
- [10] Lewis J 2003 Curr. Biol. 13 1398–1408
- [11] Schröter C, Ares S, Morelli L G, Isakova A, Hens K, Soroldoni D, Gajewski M, Jülicher F, Maerkl S J, Deplancke B and Oates A C 2012 PLoS Biol. 10 e1001364
- [12] Webb A B, Lengyel I M, Jörg D J, Valentin G, Jülicher F, Morelli L G and Oates A C 2016 eLife 5 e08438
- [13] Masamizu Y, Ohtsuka T, Takashima Y, Nagahara H, Takenaka Y, Yoshikawa K, Okamura H and Kageyama R 2006 Proc. Natl. Acad. Sci. USA 103 1313–1318
- [14] Giudicelli F, Ozbudak E M, Wright G J and Lewis J 2007 PLoS Biol. 5 e150
- [15] Aulehla A, Wiegraebe W, Baubet V, Wahl M B, Deng C, Taketo M, Lewandoski M and Pourquié O 2008 Nat. Cell Biol. 10 186–193
- [16] Morelli L G, Ares S, Herrgen L, Schröter C, Jülicher F and Oates A C 2009 HFSP J. 3 55
- [17] Ares S, Morelli L G, Jörg D J, Oates A C and Jülicher F 2012 Phys. Rev. Lett. 108(20) 204101
- [18] Delaune E A, François P, Shih N P and Amacher S L 2012 Dev. Cell 23 995–1005
- [19] Lauschke V M, Tsiairis C D, François P and Aulehla A 2013 Nature 493 101–105
- [20] Soroldoni D, Jörg D J, Morelli L G, Richmond D L, Schindelin J, Jülicher F and Oates A C 2014 Science 345 222–225
- [21] Niederreither K, McCaffery P, Dräger U C, Chambon P and Dollé P 1997 Mech. Dev. 62 67–78
- [22] Diez del Corral R, Olivera-Martinez I, Goriely A, Gale E, Maden M and Storey K 2003 Neuron 40 65–79
- [23] Moreno T A and Kintner C 2004 Dev. Cell 6 205–18
- [24] Aulehla A, Wehrle C, Brand-Saberi B, Kemler R, Gossler A, Kanzler B and Herrmann B G 2003 Dev. Cell 4 395–406
- [25] Dobbs-McAuliffe B, Zhao Q and Linney E 2004 Mech. Dev. 121 339–350
- [26] Dubrulle J and Pourquié O 2004 Nature 427 419–22
- [27] Goldbeter A, Gonze D and Pourquié O 2007 Dev. Dynam. 236 1495–1508
- [28] Shimozono S, Iimura T, Kitaguchi T, Higashijima S I and Miyawaki A 2013 Nature 496 363–6
- [29] Bajard L, Morelli L G, Ares S, Pecreaux J, Jülicher F and Oates A C 2014 Development 141 1381–1391
- [30] Bénazéraf B, François P, Baker R E, Denans N, Little C D and Pourquié O 2010 Nature 466 248–252
- [31] Ishimatsu K, Takamatsu A and Takeda H 2010 Development 137 1595–1599
- [32] Aulehla A and Pourquié O 2010 Cold Spring Harb. Perspect. Biol. 2 a000869
- [33] Murray P J, Maini P K and Baker R E 2011 J. Theor. Biol. 283 227–238
- [34] Baker R E, Schnell S and Maini P K 2006 Dev. Biol. 293 116–126
- [35] Baker R E, Schnell S and Maini P K 2006 J. Math. Biol. 52 458–482
- [36] Baker R E and Maini P K 2007 Math. Biosci. 209 30–50
- [37] Chisholm R H, Hughes B D, Landman K A, Mayer G and Whitington P M 2011 J. Theor. Biol. 279 150–160
- [38] Cotterell J, Robert-Moreno A and Sharpe J 2015 Cell Syst. 1 257–269
- [39] Tiedemann H B, Schneltzer E, Zeiser S, Rubio-Aliaga I, Wurst W, Beckers J, Przemeck G K H and Hrabě de Angelis M 2007 J. Theor. Biol. 248 120–9
- [40] Hester S D, Belmonte J M, Gens J S, Clendenon S G and Glazier J A 2011 PLoS Comp. Biol. 7 e1002155
- [41] Tiedemann H B, Schneltzer E, Zeiser S, Hoesel B, Beckers J, Przemeck G K H and Hrabě de Angelis M 2012 PLoS Comp. Biol. 8 e1002586
- [42] Jörg D J, Morelli L G, Soroldoni D, Oates A C and Jülicher F 2015 New J. Phys. 17 093042
- [43] Corless R, Gonnet G, Hare D, Jeffrey D and Knuth D 1996 Adv. Comput. Math. 5 329–359
- [44] Lecuit T, Brook W J, Ng M, Calleja M, Sun H and Cohen S M 1996 Nature 381 387–393
- [45] Chisholm R H, Hughes B D and Landman K A 2010 PLoS ONE 5 e12857