Modularity of the segmentation clock and morphogenesis

eLife Assessment

This valuable manuscript uses mathematical modeling to address the synchrony of the vertebrate segmentation clock with the developmental processes. The authors use convincing arguments to support the idea that this would allow the evolution of flexible body plans and a variable number of segments. This manuscript could be of interest to developmental biologists and systems biologists.

[Editors' note: this paper was reviewed by Review Commons.]

https://doi.org/10.7554/eLife.106316.2.sa0

Significance of the findings:

Valuable: Findings that have theoretical or practical implications for a subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Convincing: Appropriate and validated methodology in line with current state-of-the-art

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Methods
Data availability
References
Article and author information
Metrics

Abstract

Vertebrates have evolved great diversity in the number of segments dividing the trunk body, however, the developmental origin of the evolvability of this trait is poorly understood. The number of segments is thought to be determined in embryogenesis as a product of morphogenesis of the pre-somitic mesoderm (PSM) and the periodicity of a molecular oscillator active within the PSM known as the segmentation clock. Here, we explore whether the clock and PSM morphogenesis exhibit developmental modularity, as independent evolution of these two processes may explain the high evolvability of segment number. Using a computational model of the clock and PSM parameterised for zebrafish, we find that the clock is broadly robust to variation in morphogenetic processes such as cell ingression, motility, compaction, and cell division. We show that this robustness is in part determined by the length of the PSM and the strength of phase coupling in the clock. As previous studies report no changes to morphogenesis upon perturbing the clock, we suggest that the clock and morphogenesis of the PSM exhibit developmental modularity.

Introduction

Vertebrates exhibit great diversity in the number of segments that divide the skeleton, musculature, and nervous system of the body, along the rostral-caudal axis (Richardson et al., 1998). However, the developmental basis for this evolvability remains poorly understood. Segmentation of the body is established by paired blocks of mesodermal tissue, known as somites, that form periodically on either side of the developing embryo as it elongates posteriorly, in a process known as somitogenesis (Morin-Kensicki et al., 2002; Gomez and Pourquié, 2009; Oates et al., 2012). The total number of segments corresponds to the number of somites formed in the embryo, which is thought to be an emergent property of the morphogenesis of the pre-somitic mesoderm (PSM) and the dynamics of a molecular oscillator known as the segmentation clock (Morin-Kensicki et al., 2002; Cooke and Zeeman, 1976; Gomez et al., 2008; Gomez and Pourquié, 2009; Harima et al., 2013; Schröter and Oates, 2010). Here, we hypothesise that the evolvability of vertebrate segment number may be underpinned by developmental modularity of the segmentation clock and morphogenesis of the PSM (Raff, 1996), and that the potential of these two processes to evolve independently from one another may explain the diversity observed in vertebrate segment number.

The segmentation clock is a complex gene regulatory network thought to be driven by cell-autonomous oscillations of transcription factors in the Hes/Her family (Palmeirim et al., 1997; Lewis, 2003; Harima et al., 2013; Schröter et al., 2012; Webb et al., 2016). Noisy oscillations are synchronised across cells by delta-notch signalling (Horikawa et al., 2006; Lewis, 2003; Venzin and Oates, 2020), creating synchronous travelling waves of gene expression that pulse from the posterior to the anterior of the PSM. In the anterior PSM somite boundaries are patterned by the interaction of clock oscillations with a posterior-anterior decreasing gradient of FGF signalling that is thought to act as a ‘wavefront’, reading out the phase of the clock to create a spatially periodic pattern of somite boundaries (Cooke and Zeeman, 1976; Soroldoni et al., 2014; Simsek and Özbudak, 2018; Simsek et al., 2023). The clock synchronises the differentiation of PSM cells as they adopt somite fates at the anterior of the PSM (Cooke and Zeeman, 1976), and so controls both the tempo at which somites are formed and the anterior-posterior polarity of the somites (Simsek et al., 2023). Thus, clock synchrony controls the accuracy of somite patterning, and the frequency of clock oscillations determines the frequency of somite formation, which (together with the total duration of somitogenesis) determines the total number of somites formed.

Concurrent with somitogenesis, the PSM undergoes elongation. Comparative studies have shown that a diverse range of mechanisms are responsible for elongation of the PSM (Gomez et al., 2008; Bénazéraf et al., 2010; Steventon et al., 2016; Mongera et al., 2018; Thomson et al., 2021; Michaut et al., 2022). Importantly, the elongation dynamics of the PSM are thought to control the total duration of somitogenesis, which terminates once the PSM becomes sufficiently short in the AP direction (Gomez et al., 2008; Gomez and Pourquié, 2009; Steventon et al., 2016). Therefore, the total number of segments formed in the developing vertebrate is thought to be the product of both the dynamics of the clock (which controls the rate of somite formation) and the morphogenesis of the PSM, via its control on the total duration of somitogenesis.

Developmental modularity is a property of two or more developmental processes where their uncoupling in space or time permits their evolution independently of one another (Raff, 1996). This property is thought to give rise to increased phenotypic diversity by enhancing the evolvability of the system (Raff, 1996). In species examined thus far, such as the Corn snake Pantherophis guttatus, it appears that the evolution of segment number is driven by changes in both the dynamics of the clock and the elongation of the PSM (Gomez et al., 2008), and it is possible that independent evolution of these two processes is responsible for the high degree of diversity observed in vertebrate segment number. However, it is not obvious whether the two processes are modular, and if so, to what degree. Perturbing the clock’s periodicity and function does not appear to affect elongation of the body axis (Schröter and Oates, 2010; Lleras Forero et al., 2018), so morphogenesis of the PSM is likely to be robust to changes in the segmentation clock through evolution. However, it is possible that many of the cellular- and tissue-level processes which drive PSM elongation could exert an effect on the dynamics of the clock.

For instance, cell rearrangements thought to drive elongation of the PSM in zebrafish (Danio rerio) (Lawton et al., 2013; Mongera et al., 2018) and in chicken (Gallus gallus) (Bénazéraf et al., 2010; Michaut et al., 2022), promote synchronisation of the segmentation clock (Uriu et al., 2013; Uriu et al., 2010; Uriu and Morelli, 2014; Uriu et al., 2017). Additionally, arrest of transcription during chromatin condensation is known to cause PSM cells to lag their clock expression relative to neighbours after cell division (Horikawa et al., 2006; Delaune et al., 2012), causing defects in clock synchrony (Murray et al., 2013) and suggesting clock synchrony may depend on the degree of proliferative growth in the PSM. The ingression of PSM progenitor cells from dorso-posterior and lateral donor tissues (Kanki and Ho, 1997; Steventon et al., 2016; Xiong et al., 2020) can also create clock asynchrony, as incoming progenitor cells do not appear to show clock gene expression (Mara et al., 2007) and thus may be asynchronous with their neighbours when they enter the PSM. Finally, tissue convergence movements associated with the elongation of the PSM (Thomson et al., 2021) have been proposed to negate the effect of random cell mixing and be deleterious for clock synchronisation (Uriu and Morelli, 2014). It is thus non-trivial to determine whether the duration and rate of somitogenesis are modular since, in order to do so, one must examine the effect of varying morphogenesis on the dynamics of the clock.

To do this, we use a computational approach to simulate clock dynamics and cell movements within the PSM, and study how the clock responds to changes in morphogenesis. We use this approach as it is much quicker and more flexible than experimentally manipulating morphogenesis in vivo. As we are constrained by computational complexity and cannot simulate all possible means by which the PSM can elongate (this would require simulation of the growth and dynamics of cells in surrounding tissues which exert forces on the PSM; Xiong et al., 2020), we limit our study to simulating regimes of cell movement and growth within the PSM that are known to be associated with, or causal in, PSM growth and elongation across vertebrates.

We use a previously established three-dimensional model of cell movements and clock dynamics within the zebrafish PSM (see Figure 1A; Uriu et al., 2021) and adapt this model to simulate various regimes of cell movement and growth. Briefly (a more detailed description of the model can be found in the Methods), the model assumes the PSM is in an inertial frame of reference, with elongation being encoded by the advection of cells towards the anterior PSM (the anterior limit of which is denoted $x = x_{a}$ , see Figure 1A). Cells are modelled as point particles and are subject to random cell mixing, a cell-cell repulsion force that enforces volume exclusion of cells, and a boundary force confining the cells within a hollow horseshoe-shaped domain that reflects the geometry of the PSM. This model uses a phase-oscillator approach to describe segmentation clock dynamics and each cell is assigned a clock phase which is frozen when the cell exits the PSM at $x = x_{a}$ , crossing the wavefront and patterning a segment (Figure 1A). To replenish cells lost at the anterior, new cells are added with random phase and random position to keep the tissue at a constant density (Figure 1B).

Figure 1

Download asset Open asset

Computational model of the clock and the pre-somitic mesoderm (PSM).

(A) Stills of a simulation of the model of Uriu et al., 2021. Kinematic phase ( $θ_{i}$ ) waves emerge in the posterior (right) and travel towards the tissue anterior (left, $x = x_{a}$ ), where phase is arrested. The model is parameterised to data from zebrafish, and accordingly the clock oscillates every 30 min. (B) Insets illustrating the key processes driving cell movements in the PSM within the model. Top: Cells advect towards the anterior of the tissue, simulating elongation of the PSM. Middle: New cells are added to replenish the loss of cellular material as cells advect towards the anterior. Bottom: Cells undergo motility-driven rearrangements. (C) Functions in the model describing (top) the intrinsic oscillation frequency, (middle) the advection velocity, (bottom) and the motility, of each cell depending on its normalised position along the anterior-posterior axis, $(x - x_{a}) / L$ . Plots were generated using the parameters given by Uriu et al., 2021.

This model is appealing for several reasons, the first being that the model is three-dimensional and so can quantitatively recapture the rates of cell mixing that we observe within the PSM (Uriu et al., 2017; Uriu et al., 2021). The second is that the segmentation clock is described by an abstract phase-oscillator, broadly agnostic to the molecular details which are known to vary between vertebrates (Krol et al., 2011). Finally, the model makes extensive use of experimentally derived parameters for cell movements, tissue dimension, and the phase-oscillator model of the clock in zebrafish (Uriu et al., 2021). This allows us to make quantitative predictions as to how the zebrafish segmentation clock would react to changes in PSM morphogenesis.

Using this model we test the clock’s response to changes in the position of cell ingression into the PSM, changes in the anterior-posterior profile of random cell motility, changes in the length and density of the PSM, and to mitosis. We find that clock synchrony and frequency are robust to these changes except when mitosis is introduced, where we find that clock synchrony is negatively impacted except in cases where mitosis is confined to the posterior PSM. We find that this robustness is underpinned by the tissue length and cell density in the model, as well as the rate of cell mixing and strength of clock phase coupling. Together, these results suggest that segmentation clock dynamics and PSM morphogenesis are modular components of somitogenesis that can evolve independently from one another, conferring evolvability and helping to explain the diversity in segment number across the vertebrates.

Results

Varying the position of cell ingression has minimal effect on clock dynamics

During somitogenesis new cells enter the PSM posteriorly via ingression from dorsal tissues (Kanki and Ho, 1997; Banavar et al., 2021). In the early stages of zebrafish somitogenesis, there is also a contribution of cells from lateral tissues (Steventon et al., 2016). As ingressing cells do not appear to express segmentation clock genes (Mara et al., 2007), the position at which cells ingress into the PSM can create challenges for clock patterning, as only in the ‘off’ phase of the clock will ingressing cells be in-phase with their neighbours. Therefore, continuous ingression of PSM progenitor cells is likely to create local asynchrony of oscillations where it occurs, and so varying the position of cell ingression may have an effect on clock dynamics at the wavefront.

To test this hypothesis, we alter how cells are added to the tissue in the model. To model the ingression of cells from dorsal tissues into the posterior tailbud (Kanki and Ho, 1997; Banavar et al., 2021), we define a density-dependent cell addition process where cells are added onto the dorsal surface of the posterior half-toroid subdomain (see Figures 2A and 9) if the average density in this subdomain falls below $ρ_{0}$ , with a constant initial phase $θ = 0$ . To maintain equal density across the PSM, we found it necessary to add cells in the two cylindrical PSM subdomain as well, as when cells were only added to the posterior half-toroid subdomain we observed a loss in density in the two cylindrical subdomains (data not shown). To avoid bias we add cells at random positions if the average density falls below $ρ_{0}$ . These cells are initialised with a random phase drawn from the uniform distribution $[0, 2 π)$ . To model the early stages of zebrafish somitogenesis where, concurrent with dorso-posterior ingression of cells, the PSM also experiences a contribution of cells from lateral tissues, we modify the above simulation so that cell density in the lateral cylindrical subdomain is maintained by the addition of cells with constant phase ( $θ = 0$ ) on the ventral surface of these two lateral subdomain (see Figures 2A and 9). These two ingression scenarios are termed ‘dorso-posterior’ (‘DP’) and ‘dorso-posterior+lateral-ventral’ (‘DP+LV’) ingression (Figure 2A), respectively, and for brevity we refer to them as ‘DP’ and ‘DP+LV’ in the rest of the text. We compare the effect of these scenarios of cell ingression with ’Random’ ingression (Figure 2A), where cell density is maintained by the addition of cells at random positions within the tissue with random phase ( $θ \in [0, 2 π)$ ).

Figure 2 with 1 supplement see all

Download asset Open asset

The effect of cell ingression position on clock frequency and synchrony.

(A) Diagram highlighting the position of cell addition across the three ingression scenarios tested here. Magenta shading shows where cells are added onto the tissue surface. The pale magenta and blue planes respectively correspond to the anterior limit of the pre-somitic mesoderm (PSM), $x = x_{a}$ , and the anterior limit of cell addition, $x = x_{d}$ . In the ‘Random’ condition cells are added at random positions within the PSM posterior to the plane $x = x_{d}$ , and accordingly no magenta surface is shown. Similarly in the dorso-posterior case (‘DP’) cells are added at random positions in the two lateral cylinders to maintain density, and no magenta surface is shown there. In the dorso-posterior+lateral-ventral (‘DP+LV’) case, cells are only added on the tissue surface at the positions shown by magenta shading. (B) Oscillation synchrony ( $r$ ) at the PSM anterior after 1000 min, for the three tested scenarios of cell ingression. $N = 100$ simulations. (C) Mean frequency of oscillations for cells at the PSM anterior after 1000 min of simulation, for the three tested scenarios of cell ingression. $N = 100$ simulations. (D) Kymographs of synchrony along the $x$ -axis on the left-hand side of the PSM for single simulations from each of the three scenarios of cell ingression tested. Black arrowheads highlight strongly asynchronous populations of cells being transported to the tissue anterior by advection.

We test these scenarios by initialising each simulation with a synchronous initial condition for the clock and simulating cell movements and clock dynamics in the presence of each cell ingression scenario for 1000 min (see Methods for further justification). After 1000 min we measure the synchrony ( $r$ ) and mean frequency ( $d θ / d t$ ) of a cylindrical region of tissue at the left-hand PSM anterior, one cell diameter ( $d_{c}$ ) in length (see Methods). We restrict our measurement to the PSM anterior as it is here that the clock patterns somites (Simsek and Özbudak, 2018; Simsek et al., 2023) and, thus, dynamics elsewhere in the tissue are unlikely to present a phenotype in the animal. Using this method, we find that DP ingression has no effect on synchrony at the PSM anterior after 1000 min when compared to Random ingression, but that DP+LV ingression creates a minor decrease in anterior synchrony (Figure 2B). However, in neither scenario is the anterior frequency changed when compared to Random ingression (Figure 2C). We see equivalent results when ingressing cells are initiated with a random phase (Figure 2—figure supplement 1). Kymographs of synchrony along the anterior-posterior axis reveal that cell addition onto the tissue surface creates asynchronous populations of cells that travel towards the tissue anterior with cell advection (Figure 2D, black arrowheads). These are more frequent than in the DP case, likely explaining the lower anterior synchrony observed in the DP+LV simulations.

The clock is robust to changes in mode of cell ingression regardless of the cell motility profile

As cell mixing promotes clock synchrony (Uriu et al., 2013; Uriu et al., 2010; Uriu and Morelli, 2014), and an anterior-posterior profile of increasing cell mixing is present in the PSM of zebrafish and chicken (Kanki and Ho, 1997; Bénazéraf et al., 2010; Lawton et al., 2013; Mongera et al., 2018; Michaut et al., 2022) (and presumably most vertebrate taxa), we speculated that the shape of this gradient may confer robustness to clock dynamics against cell ingression. For instance, a steep profile (e.g. as reported from the early somitogenesis stages of zebrafish; Mongera et al., 2018) might confer robustness to lateral ingression as mixing is higher along a greater length of the PSM than a more graded profile (e.g. as reported from chicken embryos; Bénazéraf et al., 2010; Michaut et al., 2022). Therefore, the shape of the cell motility profile in the PSM may constrain the evolution of cell ingression position and vice versa.

To test this, we systematically enumerated the steepness ( $h$ ) and inflexion point ( $X_{v}$ ) of the function $v_{0} (x)$ (Figure 3B), which controls the speed of random cell movement in the PSM (see Methods, Figure 3A). For each parameter pair we performed $N = 100$ simulations and recorded the median anterior synchrony and median anterior mean frequency ( $d θ / d t$ , see Methods) after 1000 min. The results for each ingression scenario are shown in Figure 3C. We find no trend in median anterior synchrony and frequency across different parameter pairs for Random ingression, DP ingression, or DP+LV ingression (Figure 3C). Similar to before, results for DP ingression are indistinguishable from Random ingression. For DP+LV ingression, anterior synchrony is lower across all parameter sets than DP ingression or Random ingression (Figure 3C) but the effect is still very minor in this case. We do not observe any $X_{v}$ - or $h$ -dependent trend when we examine synchrony along the anterior-posterior axis (Figure 3—figure supplement 1). We therefore predict that, at least with respect to effects on the zebrafish segmentation clock, the evolution of the motility profile and the position of cell ingression are not constrained by one another.

Figure 3 with 1 supplement see all

Download asset Open asset

Effect of cell motility profile on clock frequency and synchrony.

(A) Overview of how intrinsic cell motion is encoded in the model. Each cell is given a random direction vector $v_{0} (x_{i}) n_{i} (t)$ (black arrows) whose magnitude $v_{0} (x_{i})$ increases towards the pre-somitic mesoderm (PSM) posterior. (B) Magnitude of intrinsic cell motion $v_{0} (x)$ for the parameters used by Uriu et al., 2021 (blue), and how the shape of the function can change when increasing (yellow) or decreasing (green) the inflexion point and curve steepness parameters, $X_{v}$ and $h$ , respectively. (C) Clock synchrony $r$ (top) and difference from expected mean frequency $Δ d θ / d t$ (bottom) for different $v_{0} (x)$ specified by combinations of $X_{v}$ and $h$ . The corresponding pixels display the synchrony or frequency at the PSM anterior after 1000 min of simulation using the specified motility profile, averaged across $N = 100$ simulations. A black + corresponds to the parameter pair used elsewhere in this paper, unless otherwise stated ( $X_{v} = 0.4, h = 3$ ). Parameter ranges used are $X_{v} \in {0, 0.1, 0.2, \dots, 1}$ and $h \in {0, 0.5, 1, \dots, 5}$ .

Length and density of the PSM confer robustness to the clock

The PSM is a transient tissue that shrinks in length over time. While this feature is conserved across vertebrate species, the absolute length and rate of shrinkage are known to vary between taxa (Gomez et al., 2008; Steventon et al., 2016; Thomson et al., 2021). We postulate that PSM length is an important factor in determining the clock’s response to different morphogenetic modes; e.g., intuitively cells will have more time to synchronise oscillations with their neighbours before reaching the PSM anterior in a longer PSM, or if cell ingression is spread out along a longer tissue the degree of asynchrony imparted by ingression will be lessened. Similarly, we predict that the density of the tissue confers robustness to clock dynamics as increasing the average number of neighbours for each cell may act to ‘correct’ against clock noise. Overall, tissue density has been observed to increase over time in zebrafish (Thomson et al., 2021) and has been observed to vary in space in chicken, with density decreasing towards the PSM posterior (Bénazéraf et al., 2010; Michaut et al., 2022).

We thus sought to explore whether the length and density of the PSM confer robustness to changing morphogenesis. We test this by co-varying the position of cell ingression (Figure 2A) with either tissue density or length. When varying the average tissue density $ρ_{0}$ and studying model dynamics after 1000 min we see that increasing cell density marginally increases anterior synchrony (Figure 4B) and decreases the variability in frequency at the PSM anterior (Figure 4C). As expected, the number of neighbours for each cell is increased (Figure 4—figure supplement 1). However, the increase in density does not appear to be sufficient to ‘correct’ the decrease in synchrony in the DP+LV case (Figure 4B). Increasing the length of the tissue does however appear to ‘correct’ this decrease (Figure 4E). Here, the PSM length $L$ is increased but the anterior limit of cell addition $x_{d}$ is held constant at $x_{d} = 100 µm$ . Anterior synchrony after 1000 min positively correlates with $L$ , and decreasing $L$ reveals a decrease in anterior synchrony for DP ingression (Figure 4E). No obvious trend for mean anterior frequency is observed (Figure 4F). Kymographs of synchrony along the $x$ -axis suggest that while much of this correlation is driven by cells having more time to synchronise their neighbours with increasing $L$ , a significant part of this may be due to cell addition occurring along a longer domain, and in so doing ‘diluting’ the noise-giving effect of cell ingression (Figure 4—figure supplement 2). To control for this, we restrict cell addition to only occur along a fixed length of the tissue with $x_{d} = L - R - r_{T} - 100$ (where $R$ represents the major radius of the half-toroid subdomain and $r_{T}$ the radius of the PSM tube, see Figure 9), and vary the length of the tissue $L$ . We observe similar results as to before, with synchrony positively correlating with $L$ and no clear trend for clock frequency (Figure 4—figure supplement 3). Kymographs of synchrony along the $x$ -axis confirm that increasing $L$ increases the effective time cells have to synchronise oscillations with their neighbours (Figure 4—figure supplement 4).

Figure 4 with 5 supplements see all

Download asset Open asset

Effect of tissue density and length.

(A) Stills from exemplar simulations illustrating the impact of changes in tissue density $ρ_{0}$ . (B) Anterior synchrony after 1000 min for changing tissue density $ρ_{0}$ and varying position of cell ingression. $N = 100$ simulations. (C) Anterior mean frequency after 1000 min for changing tissue density $ρ_{0}$ and varying position of cell ingression. $N = 100$ simulations. (D) Stills from exemplar simulations illustrating changes in tissue length. (E) Anterior synchrony after 1000 min for changing tissue length $L$ and varying position of cell ingression. $N = 100$ simulations. (F) Anterior mean frequency after 1000 min for changing tissue length $L$ and varying position of cell ingression. $N = 100$ simulations.

We note that when either one of density ( $ρ_{0}$ ) or PSM length ( $L$ ) increases, the differences in synchrony between the three scenarios of cell ingression become less pronounced (Figure 4B and E). This suggests that density and length are capable of independently conferring robustness to the segmentation clock against changes in morphogenesis, however, co-varying density and length suggests that the relative effect of length is stronger (Figure 4—figure supplement 5). In zebrafish density and length are known to co-vary after the 16-somite stage (ss), with the PSM increasing in density as it shrinks in length (Thomson et al., 2021). We would expect decreasing the length of the PSM to decrease the clock’s robustness to noise, but it is possible that the concurrent increase in density is sufficient to counteract this and preserve clock dynamics. As the dynamics of somitogenesis are consistent across individuals at these late stages (Schröter et al., 2008), the dynamics of the zebrafish segmentation clock are likely unaffected by compaction-extension of the PSM. It is not clear however if this can be explained by the concurrent changes in PSM length and density that occur during compaction-extension or if clock properties such as phase coupling need to vary in order to maintain robust clock dynamics.

Here, we simulate compaction-extension by shrinking the PSM radius $r$ and length $L$ , while increasing its density $ρ$ over time, using published rates and initial conditions for these values (Figure 5A and B, Methods) (Thomson et al., 2021). We also decrease the cell diameter $d_{c}$ , as has been observed in zebrafish (Thomson et al., 2021), as we expect this to be important in determining clock dynamics by changing the number of neighbouring cells that a given cell can couple its phase to.

Figure 5 with 3 supplements see all

Download asset Open asset

Effect of compaction-extension on clock synchrony and frequency.

(A) Snapshots of an exemplar simulation showing how the pre-somitic mesoderm (PSM) shrinks in length and diameter as time progresses. A 100 µm scalebar is shown in the lower right-hand side of the figure. (B) Functions for PSM length $L$ , radius $r_{T}$ , density $ρ$ , and cell diameter $d_{c}$ , derived from Thomson et al., 2021 (yellow), and the constant functions (blue) with which the effect of these functions is compared. (C) Anterior synchrony and (D) mean anterior frequency over time, for $N = 100$ simulations. The solid line indicates the median and the inter-quartile range is given by a shaded band either side of this line. Blue shows simulations where the tissue does not undergo compaction-extension after $t_{s h r i n k}$ , and yellow shows simulations where after $t_{s h r i n k}$ the tissue undergoes compaction-extension according to the functions shown in B. Results are plotted until the time at which at least one of the replicate simulations encounters a gap in the tissue at the tissue anterior (see Methods).

As these changes are dynamic, we depart from previous methodology for measuring clock dynamics and plot the anterior synchrony and frequency over time for $N = 100$ simulations. We also no longer simulate for 1000 min as in previous simulations, rather terminating the simulation when the length of the tissue $L (t)$ is less than $x_{d}$ (i.e. when $x_{a} (t) + x_{d} \geq X_{c} + R + r_{T} (t)$ ). We impose this constraint so as to prevent asynchrony induced by random addition of new cells biasing our results. It is important to note that we therefore only simulate for a part of the total duration of somitogenesis in zebrafish, and that there may be dynamics in the terminal stages of somitogenesis not captured by this model.

The anterior synchrony and mean anterior frequency over time for $N = 100$ simulations of compaction-extension are shown in Figure 5C and D, respectively. They are compared with simulations that do not shrink the PSM and maintain constant $L$ , $r_{T}$ , $d_{c}$ , and $ρ$ . We find that towards the end of these simulations the clock experiences more noise and fluctuations in synchrony and frequency (Figure 5C and D) than the non-shrinking equivalent simulations, however, for the majority of the simulation the dynamics of the compacting tissue are broadly comparable with those of the non-compacting tissue, but with synchrony minorly decreasing and frequency becoming noisier towards the end of the simulation (Figure 5C and D). As model parameters that are fixed throughout the simulation, such as advection velocity or intrinsic frequency, are thought to change towards the end of somitogenesis (Schröter et al., 2008; Uriu et al., 2021), it is possible that for the latter stages of somitogenesis the results of our simulations are inaccurate. Data for the terminal stages of somitogenesis are extremely limited and we cannot make estimates for model parameters at this stage. Therefore, we conclude that for much of somitogenesis in zebrafish the clock is not majorly impacted by the compaction-extension of the tissue.

Beyond dynamic changes in density, we also attempted to encode a spatial gradient of decreasing cell density towards the posterior, as reported from chicken (Bénazéraf et al., 2010; Michaut et al., 2022). However, we could not generate a gradient that quantitatively recaptures that reported by Bénazéraf et al., 2010 (data not shown). Based on our results for changing average cell density $ρ_{0}$ (Figure 4B and C), it is highly probable that this gradient would affect clock dynamics. We suggest that another model, such as one encoding cell-cell adhesion, may be a more powerful tool for exploring the effect of this gradient and dynamics at the terminal stages of somitogenesis.

Clock arrest during cell division creates asynchrony

During M-phase of the cell cycle, segmentation clock gene expression in zebrafish is known to arrest (Horikawa et al., 2006), causing a lag in nascent daughter cells relative to their neighbours (Delaune et al., 2012). If cell divisions are asynchronous, as appears to be the case in zebrafish (Kanki and Ho, 1997; Steventon et al., 2016), then such lags can create asynchrony of oscillations in the tissue (Murray et al., 2013).

As previous models of mitosis and the segmentation clock have been two-dimensional and therefore may not quantitatively recapture the cell mixing and geometry present in zebrafish (Murray et al., 2013; Murray et al., 2019), we sought to investigate the effect of mitosis in the presence of the more accurate three-dimensional geometry and cell mixing. To do so, we simulated cell division by assigning each cell a cell cycle phase $τ$ that increases at a constant rate of $1 {min}^{- 1}$ , and spawning a new cell adjacent to that cell once that its $τ \geq T_{G} + T_{M}$ , where $T_{G}$ denotes the total time in minutes to complete G1, S, and G2 phases of the cell cycle, and $T_{M}$ denotes the time spent in M-phase, in minutes. Once a daughter cell is generated, the $τ$ for both cells is reset to zero. As daughter cells share the same clock phase after division in vivo (Delaune et al., 2012), the daughter cell shares the same clock phase $θ$ as its sibling. To simulate the arrest of the clock during M-phase, if $τ \in [T_{G}, T_{G} + T_{M})$ oscillations cease, i.e., $d θ_{i} / d t = 0$ . For simplicity, we allow neighbouring cells to couple their phase to that of a neighbour with arrested phase, as transcriptional arrest is presumably equivalent to the ‘off’ phase of the segmentation clock. For further details, we refer the reader to the Methods section below.

To avoid measuring trivial changes in frequency and synchrony due to arrest of the clock, we measure synchrony and frequency of the clock only in cells where $τ < T_{G}$ . Therefore, our results represent a lower bound on the effect of mitosis on the clock. Using the estimated values of $T_{G}$ and $T_{M}$ for zebrafish (see Methods for calculation), we find that the presence of cell division creates asynchrony as reported previously (Murray et al., 2013) and has no identifiable effect on frequency (Figure 6B and C). Fixing $T_{G} + T_{M} = 187.5 min$ and varying $T_{M}$ shows a $T_{M}$ -dependent trend where synchrony decreases with increasing $T_{M}$ (Figure 6E). After $T_{M} = 15$ min the synchrony recovers slightly but remains noisy (Figure 6E). No noticeable trend is observed for the frequency (Figure 6—figure supplement 1A). We observe similar results when neighbouring cells are forbidden from coupling their phase to neighbours in M-phase (Figure 6—figure supplement 2).

Figure 6 with 3 supplements see all

Download asset Open asset

Effect of cell division on clock frequency and synchrony.

(A) Diagram showing how during cell division, clock phase $θ$ arrests, causing a cell to fall out of phase with its neighbour. On the right-hand side a still from a simulation for $T_{M}$ = 15 min is shown, illustrating how this creates asynchrony of oscillations. The effect of cell division on anterior synchrony and frequency is shown in figures **B and C**, respectively, for $T_{M} = 15$ min after 1000 min. $N = 100$ . (D) The effect on anterior synchrony when division is restricted to only occur posterior to $x = x_{d i v}$ , for $T_{M} = 15$ min. $N = 100$ . (E) The effect on anterior synchrony after 1000 min when $T_{M}$ is varied. In each case, the total length of the cell cycle is maintained at a constant length, i.e., $T_{M} + T_{G} = 187.5$ min. $N = 100$ . To rule out trivial changes in synchrony and frequency, in all analysis here we restrict measurement to non-dividing cells, i.e., cells such that $τ \in [0, T_{G})$ .

In previous simulations (Figure 6), cell division occurs up to the wavefront where segments are pre-patterned ( $x = x_{a}$ ). However, if division is confined to the posterior (Bouldin et al., 2014), cells may have time to re-synchronise before reaching the wavefront. To test whether this is possible within the timescales of zebrafish cell movement, division, and clock dynamics, we define a point in space $x = x_{a} + x_{d i v}$ anterior to which (i.e. for $x < x_{a} + x_{d i v}$ ) we halt cell cycle progression by fixing $τ = 0 min$ in all cells. Performing simulations with $x_{d i v} = 25 µm$ , $175 µm$ , $325 µm$ , we see that while for $x_{d i v} = 25 µm$ we observe defects in synchrony, for $x_{d i v} = 175 µm$ , $x_{d i v} = 325 µm$ , synchrony is rescued (Figure 6D). This suggests that there exists some point along the AP axis of the PSM, posterior to which cell division can occur without affecting the anterior dynamics of the clock. To try and resolve the value of this point and how it depends on the length of M-phase $T_{M}$ , we systematically enumerated values of $x_{d i v}$ from the tissue anterior to the tissue posterior and $T_{M}$ across the values tested in Figure 6E. We find that for $T_{M} < 15$ min, division can occur up to the tissue anterior without a major defect in synchrony (Figure 6—figure supplement 3A). For $T_{M} > 15$ min, simulations with division show high synchrony at the anterior although not so high as those with $T_{M} < 15$ min, where the value of $x_{d i v}$ such that clock synchrony is preserved is displaced to the posterior (Figure 6—figure supplement 3A). However, only for $T_{M} = 15$ min is this value majorly displaced to the posterior, with $x_{d i v} \approx 150 µm$ being the anterior-most value of $x_{d i v}$ such that synchrony is preserved ( $r \approx 1$ ) (Figure 6—figure supplement 3A). This result is perhaps intuitive, as the parameters in the model set the clock period to $2 π / ω_{0} = 30$ min, and a clock arrest lasting $T_{M} = 15$ min would be sufficient to move a cell into antiphase relative to its neighbours. Indeed, we see that mitosis is most deleterious to synchrony when M-phase lasts approximately half the length of a clock cycle (Figure 6—figure supplement 3B), however, for long clock cycles this relationship becomes nonlinear and the most deleterious values of $T_{M}$ lie between the length of one or two clock cycles (Figure 6—figure supplement 3B).

Zebrafish clock coupling and cell mixing confer robustness to changes in cell ingression

Across many of the morphogenetic scenarios tested here, we have observed a high degree of robustness in clock dynamics with respect to changing cell movements and processes. As shown above, some of this robustness can be attributed to the choice of tissue length $L$ in simulations and choices of the parameters $v_{a}$ and $v_{p}$ governing the advection velocity of cells. However, we also expect that the rate of cell mixing (Uriu et al., 2013; Uriu et al., 2017; Uriu et al., 2010) and the strength of clock coupling confer robustness too. Within the present model, the global rate of cell mixing is controlled by the parameter $v_{s}$ and the clock coupling strength is denoted by $κ$ . Notably, the values taken for these parameters have been experimentally measured in zebrafish (Riedel-Kruse et al., 2007; Herrgen et al., 2010; Uriu et al., 2017; Uriu et al., 2021). Therefore, we were interested in where these parameters lie within the space of parameter pairs that generate clock dynamics robust to changes in morphogenesis - do these results reflect something unique about the zebrafish parameters or can many ( $κ$ , $v_{s}$ ) pairs confer robustness?

To determine this, we take changing cell ingression as an example of varying morphogenesis and for each ingression scenario, systematically enumerate the parameters $v_{s}$ and $κ$ , studying the anterior synchrony after 1000 min with N=100 replicates per parameter pair. This reveals a threshold curve of ( $κ$ , $v_{s}$ ) pairs below which the clock begins to exhibit asynchrony (Figure 7). Comparing across the three scenarios of cell ingression, we see that the position of this curve changes with the cell ingression scenario (Figure 7), indicating that the space of ( $κ$ , $v_{s}$ ) pairs that achieve clock synchrony is constrained by the cell ingression mode present within the tissue. Notably, the experimentally derived parameter pair ( $κ = 0.07 {min}^{- 1}$ , $v_{s} = 1 µm \cdot {min}^{- 1}$ ) lies very close to the threshold for the DP+LV case of cell ingression in the $κ$ direction, suggesting that the robustness we have observed in response to changes in cell ingression is sensitive to the choice of $κ$ . However, in each case our experimental choices for $κ$ and $v_{s}$ lie above the threshold, implying that at least part of the robustness we observed with respect to cell ingression is due to the experimental values for zebrafish being sufficient to confer this robustness.

Figure 7 with 1 supplement see all

Download asset Open asset

Synchrony for varying coupling strength $κ$ and magnitude of intrinsic cell motion.

$v_{s}$ . Anterior synchrony after 1000 min, for varying maximum magnitude of intrinsic cell motion $v_{s}$ and clock phase coupling strength $κ$ , for three different scenarios of cell ingression.Each pixel corresponds to the median value of anterior synchrony across $N = 100$ simulations. A black + marks the experimental values for zebrafish, $κ = 0.07 {min}^{- 1}$ , $v_{s} = 1 µm \cdot {min}^{- 1}$ , derived in Riedel-Kruse et al., 2007, and Uriu et al., 2017, respectively, that are used elsewhere in this paper. All other parameters are held constant at their normal values (see Table 1).

Finally, we test if this robustness holds in the presence of a phase coupling delay between cells. We add to our simulations the time taken for information about the transcriptional state of a cell’s neighbour to be transduced by cell-cell signalling and impact the cell’s transcriptional profile. Such delays are important to consider as they can be significant for clock synchrony and frequency (Lewis, 2003; Morelli et al., 2009), however, they significantly increase the computational complexity of such a model. The value of the cell-cell phase coupling delay has been estimated to be $21 min$ in zebrafish (Herrgen et al., 2010), and so to test the effect of delays we repeat the analysis above (Figure 7) in the presence of a cell-cell coupling delay ( $τ$ ) of $21 min$ .

We observe that, for $τ = 21$ min, oscillations are asynchronous in all three scenarios of cell ingression with the experimentally measured values of $κ$ and $v_{s}$ , however, for larger values of $κ$ (and, to a lesser extent, $v_{s}$ ), we observe synchronous oscillations (Figure 8B). Notably, we cannot achieve synchronous oscillations for any of the chosen parameter pairs for the DP+LV scenario of cell ingression (Figure 8B), suggesting that ingression of progenitor cells is a far greater source of noise for the clock in the presence of coupling delay than in the instantaneous coupling case. To test this, we varied the length of the PSM $L$ and either restricted the addition of cells to occur within a region of length $100 + R + r_{T}$ at the tissue posterior or allowed cell addition to occur at any point in the tissue posterior to $x_{d}$ = $100 µm$ . In either case the rate of cell addition is assumed to be constant and thus in the latter case we expect the longer tissue domain to ‘dilute’ the noise effect of cell ingression. We observe no difference between either scenario and find no value of $L$ such that synchronous oscillations are recovered for $τ = 21 min$ (Figure 8—figure supplements 2 and 3).

Figure 8 with 3 supplements see all

Download asset Open asset

Clock dynamics in the presence of a coupling delay between cells.

(A) Diagram illustrating phase coupling with delay. Cells couple their phase to the phase of their neighbours $τ$ min ago. (B) Anterior synchrony in the presence of coupling delay ( $τ = 21 m i n$ ) after 1000 min, for varying maximum magnitude of intrinsic cell motion $v_{s}$ and clock phase coupling strength $κ$ , for three different scenarios of cell ingression. Each pixel corresponds to the median value of anterior synchrony across $N = 100$ simulations. A black + marks the experimental values for zebrafish, $κ = 0.07 {min}^{- 1}$ , $v_{s} = 1 µm \cdot {min}^{- 1}$ , derived in Riedel-Kruse et al., 2007, and Uriu et al., 2017, respectively, that are used elsewhere in this paper. All other parameters are held constant at their normal values (see Table 1). (C) Anterior synchrony (top) $r$ after 1000 min, for varying phase coupling delay $τ$ and clock phase coupling strength $κ$ , for three different scenarios of cell ingression. Each pixel corresponds to the median value of anterior synchrony across $N = 100$ simulations. Bottom: Difference from expected mean frequency ( $Δ d θ / d t$ ) at the pre-somitic mesoderm (PSM) anterior after 1000 min, for three different scenarios of cell ingression. Each pixel corresponds to the mean value across $N = 100$ simulations. A black + marks the experimental values for zebrafish, $κ = 0.07 {min}^{- 1}$ , $τ = 21 min$ , derived in Riedel-Kruse et al., 2007, and Herrgen et al., 2010, respectively. All other parameters are held constant at their normal values (see Table 1).

To explore how the co-evolution of the clock coupling strength $κ$ and the coupling delay $τ$ might constrain the position of cell ingression, we systematically enumerated values of $κ$ as before while co-varying $τ$ , and analysed clock synchrony across the three different scenarios of cell ingression (Figure 8C). In agreement with previous modelling work (Morelli et al., 2009; Herrgen et al., 2010), we see that synchrony depends on the value of $τ$ in a repeating manner, with values of $τ$ closest to half of a clock cycle (~ 30 min) being more asynchronous (Figure 8C). These results also predicted that frequency decreases with increasing $τ$ within regions of stability (Morelli et al., 2009; Herrgen et al., 2010). We see evidence of this when we examine the collective frequency at the PSM anterior, where within each region of synchrony the frequency decreases as $τ$ increases (Figure 8C). We note that for small values of $τ$ the DP+LV scenario of cell ingression is synchronous, however, for larger values of $τ$ oscillations are largely asynchronous, unlike the DP scenario of ingression (Figure 8C).

We note that, despite our model failing to achieve synchronous oscillations for $τ = 21$ min, only relatively minor increases of $τ$ are required to achieve synchrony in the Random and DP+LV scenarios of cell ingression (Figure 8C). This closely resembles the analytical results of Herrgen et al., 2010, who found that the measured delay $τ = 21$ min lies towards the extreme end of values of $τ$ for which stable oscillations can occur (Herrgen et al., 2010). Together these results suggest that the zebrafish segmentation clock might exist in a state of criticality and close to a region of parameter space where oscillations are asynchronous. Overall our results suggest that coupling strength and coupling delay strongly constrain the evolvability of cell ingression, and that robustness to morphogenesis is conferred only for specific combinations of these two parameters.

Discussion

Here, we investigated whether morphogenesis of the PSM and the segmentation clock exhibit developmental modularity, as this could explain the high degree of evolvability observed in vertebrate segment number (Raff, 1996). We tested a broad range of PSM cellular behaviours and mechanisms that are associated with the elongation of the PSM across vertebrate species, to see if they had an effect on the dynamics of the segmentation clock. We predict that the dynamics of the zebrafish segmentation clock, specifically synchrony and frequency, are generally robust to changes in these mechanisms such as cell ingression, motility, and (under certain conditions) division. As PSM morphogenesis is independent of clock function (Schröter and Oates, 2010; Lleras Forero et al., 2018), our results suggest that the clock and morphogenesis of the PSM exhibit developmental modularity.

Our results suggest that a major determinant of this robustness is the clock coupling $κ$ , which in zebrafish appears to be sufficiently strong to confer robustness to changes in cell ingression (Figure 7). While synchrony of the clock shows a dependence on cell mixing rate ( $v_{s}$ ) and tissue length ( $L$ ), the dependence of synchrony on these parameters is weak compared to its dependence on coupling strength (Figure 7, Figure 4—figure supplement 1), suggesting that this robustness may be more dependent on properties of the segmentation clock than those of the tissue. Due to a lack of quantitative understanding of PSM morphogenesis, it is difficult to predict whether the requirements for robustness imposed on cell mixing or tissue length and density constrain the space of possible morphogenetic processes that can be altered without perturbing the clock. However, our results highlight that these requirements are themselves constrained by the strength of clock phase coupling, so it is equally possible that vertebrate segmentation clocks have strong phase coupling generally, and these requirements on morphogenesis are sufficiently relaxed to not constrain the evolution of PSM elongation. Further work, investigating the evolution of phase coupling strength in vertebrates, and developing accurate quantitative models of PSM morphogenesis, will be necessary to determine this.

It is also worth noting that while in the model formulation $κ$ represents the strength of delta-notch signalling between cells, the experimental value for this parameter is measured at a tissue level, incorporating coupling effects from other pathways, e.g., cell mixing (Riedel-Kruse et al., 2007). It is therefore possible that the value of $κ$ measured by Riedel-Kruse et al., 2007, represents an overestimate of the strength of delta-notch signalling and therefore our results present an overestimate of robustness to morphogenesis. However, our results appear less sensitive to cell mixing ( $v_{s}$ ) than $κ$ (Figure 7), so we suggest that the contribution of cell mixing to the estimate of $κ$ may be minor. Furthermore, as the precise value of $r$ necessary for correct segment patterning is unknown, and as downstream processes in segmentation may act to correct any defects in clock patterning (Naganathan et al., 2022), somitogenesis may be more robust than we assume here. Therefore, we suggest that our findings are robust to any overestimates of $κ$ by incorporation of cell mixing.

Our results when considering a phase coupling delay between cells (Figure 8) suggest that evolution of cell ingression is constrained by both the clock coupling delay between cells and vice versa. This suggests that the robustness of the segmentation clock to changes in morphogenesis is mediated by the combination of the cell-cell coupling delay and the coupling strength of the clock. While little is known about what determines the coupling strength of the clock, coupling delay is thought to be at least in part controlled by the number of intermediary components involved in signal transduction (Lewis, 2003; Yoshioka-Kobayashi et al., 2020). Therefore, the morphogenesis of the PSM may constrain the ability of the segmentation clock to undergo developmental system drift, at least within the components of the delta-notch pathway.

We note that our model fails to recover synchrony for the experimentally measured values of $κ$ and $τ$ , and thus may not accurately model the in vivo condition in some way. It is possible that the way in which we simulate cell ingression introduces too much noise into the clock, for instance the rate of progenitor cell ingression may be much higher than that in vivo, and indeed we note that elongation is driven in vivo by a combination of processes, rather than solely ingression, as we implicitly assume here. If this is true, this implies that clock imposes a constraint on the degree to which ingression is capable of driving PSM elongation. This possibility is worthy of further study and is perhaps better suited to models where PSM elongation is not dependent on cell ingression, as we assume here.

Furthermore, as discussed above, we note that for relatively small increases in the value of $τ$ , e.g., for $τ = 25 min$ , it is possible to achieve synchronous oscillations in both the random and DP scenarios of cell ingression, and partially synchronous oscillations in the DP+LV scenario of cell ingression (Figure 8C). This could suggest that the biological parameters controlling delay and cell ingression exist in a state of criticality, i.e., lying close to a region of parameter space where dynamical behaviour abruptly changes or becomes chaotic. Such states can have important implications for evolvability of a biological system (Kauffman, 1993; Verd et al., 2019), however, it is not clear if such a state is functional in this case. Further study, perhaps with models unlike the mean-field approximation we assume here, may reveal what implication, if any, this has for the evolvability of PSM morphogenesis.

Our analysis is restricted to the modes of PSM morphogenesis that are known to occur in vertebrates, and we cannot rule out the possibility that the robustness we observe merely reflects evolution of mechanisms of PSM elongation under selective pressure to minimise their effect on the clock, and that the clock may not be robust to all possible modes of PSM elongation. Nevertheless, we argue that our results demonstrate there exists a qualitatively distinct set of processes (cell ingression, division, motility, and compaction-extension) that underpin PSM elongation and do not affect clock dynamics. As vertebrate species appear to employ a combination of these processes to varying degrees (Bénazéraf et al., 2010; Steventon et al., 2016; Mongera et al., 2018; Xiong et al., 2020; Thomson et al., 2021), it is possible that evolution alters each of these processes independently to generate diversity in the elongation dynamics of the PSM (Gomez et al., 2008; Steventon et al., 2016). Therefore, this relatively small set of processes may be sufficient to explain the high degree of evolvability in vertebrate segment number.

As clock dynamics appear to be robust to changes in morphogenesis, we suggest that the clock and morphogenesis of the PSM comprise two developmental modules within zebrafish, and that this is true more generally whenever the phase coupling of the clock is sufficiently strong, cell mixing sufficiently rapid, or the PSM sufficiently dense and long. We suggest that this may have contributed to the evolution of broad diversity in vertebrate segment number (Richardson et al., 1998) and be a present source of evolvability in vertebrates. We note that the evolution of modularity is an active area of research (Wagner and Altenberg, 1996; Wagner et al., 2007) and suggest that studying the evolution of delta-notch signalling in vertebrates, and the evolutionary origin of the PSM and somitogenesis, may be an illuminating model paradigm for this field. We predict that the modularity of the clock and morphogenesis of the PSM permits a synergy between these two processes that heightens the evolvability of segment number in vertebrates, and that this may be responsible for the high diversity we observe in this trait (Richardson et al., 1998).

Our results suggest that the zebrafish segmentation clock is robust to most changes in PSM morphogenesis, and that this robustness is determined by the strength of clock phase coupling and the phase coupling delay, as well as morphogenetic parameters such as the rate of cell mixing, and the length and density of the PSM. This suggests that, at least in the regions of parameter space where these conditions for robustness hold, PSM morphogenesis is not constrained by the clock and free to evolve, and the clock and morphogenesis of the PSM exhibit developmental modularity. The precise topology of this space remains unclear and poses an intriguing question for further study. We suggest that this modularity may be responsible for the broad diversity of segment number that we observe in the vertebrates (Richardson et al., 1998) and indeed is a present source of evolvability in the clade. We note that the evolution of modularity is an active area of research (Wagner and Altenberg, 1996; Wagner et al., 2007), and while several well-characterised examples of modularity exist (Monteiro et al., 2003; Mallarino et al., 2011; Verd et al., 2019), there are few examples to our knowledge where the existence of such modularity relies on a handful of experimentally tractable parameters. We therefore suggest that studying the evolution of delta-notch signalling in vertebrates, and the evolutionary origin of the PSM and somitogenesis, may be an illuminating model paradigm for this field.

Methods

Computational model of the segmentation clock and PSM

To simulate the movements of cells within the PSM and the dynamics of the segmentation clock, we used the model of Uriu et al., 2021, who sought to explain defective somite patterning after loss of Notch-Delta signalling. Here, an abstract phase-oscillator Kuramoto model is used to describe intracellular clock oscillations and cell-cell phase coupling. Cells are confined within a three-dimensional tissue domain whose geometry approximates that of the PSM.

This model is suitable for our purposes because it provides a three-dimensional model of cell movements within the PSM, allowing us to alter morphogenesis in order to test hypotheses. Additionally, as the segmentation clock here is described as an abstract phase-oscillator, the model has broad applicability to different vertebrate species which are known to differ in the structure and regulation of their segmentation clock gene regulatory networks (Krol et al., 2011). Finally, the validity of the model has been shown in several ways via comparison with data for zebrafish (Webb et al., 2016; Liao et al., 2016; Uriu et al., 2021) and uses parameters experimentally inferred from zebrafish embryos (Riedel-Kruse et al., 2007; Soroldoni et al., 2014; Uriu et al., 2021). This gives an accurate read-out of the effect of changing morphogenesis on the segmentation clock of a vertebrate.

The model is described elsewhere (Uriu et al., 2021), but for completeness we include an account of it here. For clarity, vector variables are shown in bold, and scalar variables in normal font.

Tissue geometry and frame of reference

The model assumes that the posterior tip of the PSM is held in an inertial frame of reference and models the posterior-ward elongation of the PSM by movement of cells towards the tissue anterior.

Cells are confined within a horseshoe-shaped domain resembling the PSM. This domain is comprised of two cylindrical subdomains and one half-toroid subdomain, each representing the lateral sides of the PSM and tailbud, respectively. A schematic can be seen in Figure 9. Cells are free to move between each subdomain and leave the tissue at the anterior at $x = x_{a}$ .

Figure 9

Download asset Open asset

Geometry of the pre-somitic mesoderm (PSM) assumed in the present model.

Top left: Major axes ( $x$ , $y$ , $z$ ) used in the model. $x$ corresponds to the anterior-posterior axis of the embryo and increases towards the tissue posterior, $y$ corresponds to the left-right axis and increases to the right-hand side of the tissue, and $z$ corresponds to the dorsal-ventral axis of the embryo. Top right: Schematic of the PSM in the xy plane. The PSM is comprised of two cylinders, centred at $y = r_{T}$ and $y = 2 R + r_{T}$ , respectively, with radius $r_{T}$ . The ‘tailbud’ is represented as a half-torus subdomain centred at $x = (X_{c}, Y_{c}, Z_{c})^{T}$ , with minor radius $r_{T}$ and major radius $R$ . Bottom: Cross sections of the tissue showing how a point $x_{i}$ is assigned the polar coordinates $r_{i}$ and $q_{i}$ , in both the PSM cylinders ( $x_{i} \leq X_{c}$ ) and the half-toroid tailbud ( $x_{i} > X_{c}$ ). Adapted from Uriu et al., 2021.

The centre of the toroid subdomain is denoted by the coordinate vector $(X_{c}, Y_{c}, Z_{c})^{T}$ , and the anterior of the PSM is demarcated by the x-coordinate $x = x_{a}$ (see Figure 9). The major radius of the torus is denoted by $R$ and the minor radius (and the radius of the two cylindrical subdomains) is denoted by $r_{T}$ (see Figure 9). The resulting length of the PSM in the $x$ direction, $L$ , is therefore given by $L = r_{T} + R + X_{c} - x_{a}$ .

Cell movement

To approximate cell mixing and the advection of cells out of the PSM towards the anterior (in the chosen tailbud-inertial frame of reference), cell movement in the model is described by an equation of motion that incorporates cell advection, cell volume exclusion, random cell motion, and a boundary force confining cells within the tissue domain (Equation 1).

Specifically, each cell is assigned a position vector $x_{i} \in R^{3}$ , where $x_{i} = (x_{i}, y_{i}, z_{i})^{T}$ . This position $x_{i}$ corresponds to the cell centre. Cells are assumed to be spheres with diameter $d_{c}$ . The equation of motion for cell $i$ is

\frac{d x_{i}}{d t} = v_{d} (x_{i}) + v_{0} (x_{i}) n_{i} (t) + \sum_{j = 1, j \neq i}^{N} F (x_{i}, x_{j}) + F_{b} (x_{i}),

where $v_{d} (x_{i})$ represents the advection of the cell towards the PSM anterior, $v_{0} (x_{i}) n_{i} (t)$ the intrinsic random motion of the cell, $F (x_{i}, x_{j})$ the repulsion force between cell $i$ and cell $j$ , and $F_{b} (x_{i})$ the boundary force confining the cell within the tissue.

The advection velocity $v_{d} (x_{i})$ is given by

v_{d} (x_{i}) = (- v_{d} (x_{i}), 0, 0)^{T},

where the function $v_{d} (x_{i})$ describes the local strain rate along the AP axis:

v_{d} (x_{i}) = {\begin{cases} - \frac{v_{a} - v_{p} (1 - x_{q})}{x_{q}} \cdot \frac{x_{i} - x_{a}}{L} + v_{a} & \frac{x_{i} - x_{a}}{L} \leq x_{q}, \\ - \frac{v_{p} (x_{i} - x_{a})}{L} + v_{p} & \frac{x_{i} - x_{a}}{L} > x_{q}, \end{cases}

The model assumes that the random cell motility observed within the PSM is intrinsic, i.e., cells would exhibit random motility in isolation from the tissue. This is modelled using the direction vector $n_{i} (t)$ , which evolves according to a random walk on the surface of the unit sphere specified by the following differential equation:

\frac{d n_{i}}{d t} = \sqrt{2 D} ξ_{ϕ_{i}} (t) m_{x} + \sqrt{2 D} ξ_{φ_{i}} (t) m_{y} - 2 D n_{i} (t),

where the orthogonal vectors (for $e_{z} = (0, 0, 1)^{T}$ )

m_{x} = \frac{n_{i} (t) \times e_{z}}{| n_{i} (t) \times e_{z} |}, m_{y} = \frac{m_{x} \times n_{i} (t)}{| m_{x} \times n_{i} (t) |},

together define a plane tangent to the unit sphere, and $ξ_{ϕ_{i}} (t)$ , $ξ_{φ_{i}} (t)$ are white Gaussian noise terms that satisfy, for any $t$ , $t^{'}$ , $⟨ ξ_{ϕ_{i}} (t) ⟩ = ⟨ ξ_{φ_{i}} (t) ⟩ = ⟨ ξ_{ϕ_{i}} (t) ξ_{φ_{j}} (t) ⟩ = 0$ , and $⟨ ξ_{ϕ_{i}} (t) ξ_{ϕ_{j}} (t) ⟩ = ⟨ ξ_{φ_{i}} (t) ξ_{φ_{j}} (t) ⟩ = δ_{i j} δ (t - t^{'})$ . A derivation for Equation 4 is given by Uriu et al., 2021.

To model the increase in random cell motion towards the PSM posterior, the magnitude of intrinsic random cell movement, $v_{0} (x)$ , increases towards the posterior according to Equation 6:

v_{0} (x_{i}) = \frac{v_{s}}{1 + (\frac{1 - \frac{x_{i} - x_{a}}{L}}{X_{v}})^{h}},

The parameters $v_{s}$ , $X_{v}$ , and $h$ determine the maximum magnitude of intrinsic cell motion, profile inflexion point, and steepness of the profile, respectively. By varying these parameters, one can create a variety of motility profiles with different shapes and amplitudes (see Figure 3B). Here, the x-position of cell $i$ is normalised according to the length of the PSM $L$ , and so the profile of cell motility scales with the length of the PSM. This is a valid assumption as the length of the PSM is thought to be partially specified by FGF signalling (Simsek and Özbudak, 2018), and the random motility of PSM cells is also under the control of FGF (Bénazéraf et al., 2010).

Cell volume exclusion is encoded by repulsion of two cells if their centres are within a cell diameter $d_{c}$ of one another, according to the cell-cell repulsion force $F (x_{j} - x_{i})$ :

F (x_{j} - x_{i}) = F (x_{j} - x_{i}) \frac{x_{j} - x_{i}}{| x_{j} - x_{i} |},

where the magnitude of the cell-cell repulsion force $F (x_{j} - x_{i})$ is given by

F (x_{j} - x_{i}) = {\begin{cases} μ (\frac{| x_{j} - x_{i} |}{d_{c}} - 1) & | x_{j} - x_{i} | \leq d_{c}, \\ 0, & | x_{j} - x_{i} | > d_{c}, \end{cases}

and the parameter µ controls the strength of cell-cell repulsion throughout the tissue. For computational simplicity, cell-cell adhesion is not encoded within this model.

The form of the boundary force $F_{b} (x_{i})$ varies depending on whether the cell $i$ is in one of the lateral PSM cylindrical subdomains or in the half-toroid tailbud (Figure 9). If cell $i$ lies within the cylinders, one can use the cylindrical coordinates for cell $i$ , $r_{i}$ , and $q_{i}$ (see Figure 9), to naturally define the boundary force:

F_{b} (x_{i}) = (\begin{matrix} 0 \\ - μ_{b} e^{- \frac{δ y}{r_{b}}} \cos q_{i} \\ - μ_{b} e^{- \frac{δ z}{r_{b}}} \sin q_{i} \end{matrix}),

where

δ y = | (r_{T} - r_{i}) \cos q_{i} |,

δ z = | (r_{T} - r_{i}) \sin q_{i} | .

The magnitude of the boundary force is denoted by $μ_{b}$ and the parameter $r_{b}$ determines the length-scale of the boundary force. In the half-toroid tailbud the position $x_{i}$ of each cell $i$ can be described according to the coordinates $r_{i}$ , $p_{i}$ , and $q_{i}$ :

x_{i} = (\begin{matrix} X_{c} + R \cos p_{i} + r_{i} \cos p_{i} \cos q_{i} \\ Y_{c} + R \sin p_{i} + r_{i} \sin p_{i} \cos q_{i} \\ Z_{c} + r_{i} \sin q_{i} \end{matrix}) .

One can then naturally define

δ x = | (r_{T} - r_{i}) \cos p_{i} \cos q_{i} |,

δ y = | (r_{T} - r_{i}) \sin p_{i} \cos q_{i} |,

δ z = | (r_{T} - r_{i}) \sin q_{i} |,

describing the distance of cell $i$ from the tissue surface in the $x$ , $y$ , and $z$ directions, respectively. The boundary force then becomes

F_{b} (x_{i}) = (\begin{matrix} - μ_{b} e^{- \frac{δ x}{r_{b}}} \cos p_{i} \cos q_{i} \\ - μ_{b} e^{- \frac{δ y}{r_{b}}} \sin p_{i} \cos q_{i} \\ - μ_{b} e^{- \frac{δ z}{r_{b}}} \sin q_{i} \end{matrix}) .

Note that in both definitions of $F_{b} (x_{i})$ , the tissue subdomains are open-ended, i.e., cells are free to move between subdomains and out of the PSM at the anterior.

For further details we direct the reader to consult Uriu et al., 2021.

Cell addition

As cells are removed from the tissue by advection, new cells must be added to replenish the PSM. Loss of cells is monitored by comparing the total density of cells in each of the subdomains with a target density, $ρ_{0}$ . If the density in any of these subdomains falls below $ρ_{0}$ , a new cell is added to that subdomain. If the density of more than one of the subdomains falls below $ρ_{0}$ then a cell is only added to the subdomain with the lowest density. The position and phase of the nascent cell both vary across different simulations discussed here. Below is a description of how these quantities vary.

In simulations where cell addition is ‘Random’, cells are added with a random position within each subdomain, with a random phase, i.e., if a cell is added within one of the cylindrical subdomains, it is given a random $x_{i} \in [x_{a} + x_{d}, X_{c}]$ , a random $r_{i} \in [0, r_{T}]$ , and $q_{i} \in [0, 2 π)$ . Alternatively, if the cell is added within the half-toroid subdomain, it is given a random $p_{i} \in [- π / 2, π / 2]$ , and random $r_{i}$ and $q_{i}$ , defined as above. The cell is also assigned a random phase $θ_{i} \in [0, 2 π)$ . This method of cell addition is unbiased in terms of its effect on the model but is a biologically implausible scenario of cell ingression, so we use it here as a negative control.

In simulations where we seek to more accurately model cell ingression, cells are added on the cell surface, i.e., $r_{i} = r$ . Additionally, the surface on which cells can be added is restricted. For the half-toroid subdomain, we restrict cell addition to the dorsal quarter of the subdomain, i.e., $q_{i} \in [π / 4, 3 π / 4]$ , $p_{i} \in [- π / 2, π / 2]$ . For the cylindrical subdomains, we restrict addition to the bottom quarter of the surface, i.e., $q_{i} \in [5 π / 4, 7 π / 4]$ . To more accurately model conditions of cell ingression, combinations of these three modes of cell addition are employed (Figure 2).

In newly added cells the direction vector n is initialised with a random vector of unit length. For simulations where cell division is present, cells are added with a random cell cycle phase $τ \in [0, T_{G} + T_{M})$ .

Coupled oscillator model of the segmentation clock

Here, only the phase of each cell’s segmentation clock is considered, and oscillation amplitudes are neglected. This simplification is appropriate when oscillator coupling is weak (Kuramoto, 1984). Here, oscillators have a spatially dependent intrinsic frequency $ω (x_{i})$ and couple to adjacent cells (defined as those within $d_{c}$ , a cell’s diameter) with a magnitude controlled by the scalar value $κ$ :

\frac{d θ_{i}}{d t} = ω (x_{i}) + \frac{κ}{N_{i} (t)} \sum_{| x_{j} - x_{i} | \leq d_{c}} \sin (θ_{j} - θ_{i}) + \sqrt{2 D_{θ}} ξ_{θ_{i}} (t) .

Oscillations are noisy, with phase noise intensity $D_{θ}$ . The term $ξ_{θ_{i}} (t)$ describes Gaussian noise as before. $N_{i} (t)$ is the number of cells within a radius of length $d_{c}$ of cell $i$ at time $t$ , i.e., $N_{i} (t) = | {j | | x_{j} (t) - x_{i} (t) | \leq d_{c}} | .$

The intrinsic frequency $ω (x_{i})$ decreases towards the PSM anterior, creating a frequency gradient and travelling phase waves that move towards the anterior (Soroldoni et al., 2014). The intrinsic frequency is defined as:

ω (x_{i}) = ω_{0} (σ + (1 - σ) \cdot \frac{1 - e^{- \frac{k (x - x_{a})}{L}}}{1 - e^{- k}}),

where $ω_{0} > 0$ corresponds to the frequency of cells at the posterior of the tissue, $σ > 0$ the fold-change in frequency at the anterior compared to the posterior, and $k > 0$ controls the shape of the gradient. For simplicity we assume that $ω (x_{i})$ is constant over time, as cells at the PSM posterior are thought to maintain a constant frequency over time (Giudicelli et al., 2007; Webb et al., 2016).

To simulate the effect of the wavefront, where the clock phase is interpreted by cells to pre-pattern somite boundaries (Cooke and Zeeman, 1976), beyond $x = x_{a}$ oscillations arrest, forming stable patterns of phase (Figure 2A), i.e., $d θ_{i} / d t = 0$ if $x_{i} < x_{a}$ .

Modelling the segmentation clock with a phase coupling delay

For simulations where we seek to model the effects of a coupling delay between cells, the phase of cells evolves according to

\frac{d θ_{i}}{d t} = ω (x_{i}) + \frac{κ}{N_{i} (t)} \sum_{| x_{j} - x_{i} | \leq d_{c}} \sin (θ_{j} (t - τ) - θ_{i} (t)) + \sqrt{2 D_{θ}} ξ_{θ_{i}} (t),

where $τ$ is a constant representing the temporal delay in exchange of phase information between cells, after Morelli et al., 2009. In such simulations we initialise the phase of each newly added cell with the initial conditions $θ_{i} (t \leq t_{a d d}) = 0$ , where $t_{a d d}$ is the time at which cell $i$ was added to the tissue.

Cell division and segmentation clock arrest

During M-phase of the cell cycle, the segmentation clock arrests (Horikawa et al., 2006; Delaune et al., 2012). To simulate this, we assign each cell $i$ a cell cycle phase $τ \in [0, T_{G} + T_{M})$ , where $T_{G}$ and $T_{M}$ are the durations of G1-G2 and M phases in minutes, respectively, and evolve $τ$ according to

τ = (t - t_{a d d i t i o n}) \mod (T_{G} + T_{M}),

where $t_{a d d i t i o n}$ denotes the time at which the cell was added.

To simulate the arrest of segmentation clock gene expression during M-phase of the clock, we evolve $θ_{i}$ according to

\frac{d θ_{i}}{d t} = H (T_{G} - τ_{i}) \cdot (ω (x_{i}) + \frac{κ}{N_{i} (t)} \sum_{| x_{j} - x_{i} | \leq d_{c}} \sin (θ_{j} - θ_{i}) + \sqrt{2 D_{θ}} ξ_{θ_{i}} (t)),

where $H (x)$ is the Heaviside step function:

H (x) = {\begin{cases} 1 & x > 0, \\ 0 & x \leq 0. \end{cases}

For simulations where we forbid cells from coupling their phase to neighbours undergoing mitosis, the evolution of $θ_{i}$ is given by

\frac{d θ_{i}}{d t} = H (T_{G} - τ_{i}) \cdot (ω (x_{i}) + \frac{κ}{N_{i} (t)} \sum_{| x_{j} - x_{i} | \leq d_{c}, τ_{j} < T_{G}} \sin (θ_{j} - θ_{i}) + \sqrt{2 D_{θ}} ξ_{θ_{i}} (t)) .

If $τ_{i} \geq T_{G} + T_{M}$ then we consider the cell cycle complete and a new cell is added at a random position of length $d_{n e w}$ away from cell $i$ , i.e.,

x_{N_{t o t a l} + 1} = x_{i} + d_{n e w} (\begin{matrix} \cos φ \sin ϕ \\ \sin φ \sin ϕ \\ \cos ϕ \end{matrix}),

for uniformly distributed random $ϕ \in [0, π]$ , $φ \in [0, 2 π)$ , where $N_{t o t a l}$ is the total number of cells within the tissue. Other variables are inherited from cell $i$ , i.e., $n_{N_{t o t a l} + 1} = n_{i}$ and $θ_{N_{t o t a l} + 1} = θ_{i}$ . The choice of $d_{n e w}$ is somewhat arbitrary, however, as this distance decreases the magnitude of the cell-cell repulsion force (Equation 8) grows very large, and one must choose a $d_{n e w}$ sufficiently small to be biologically accurate while large enough to avoid implausible cell velocities due to a very large cell-cell repulsion force. We find that $d_{n e w} = d_{c} / 10$ is a suitable compromise between these two scenarios.

We are not aware of any direct measurements of $T_{G}$ within the zebrafish PSM, and so have calculated a value for this parameter based on some assumptions. M-phase within the PSM lasts approximately 15 min, and approximately 8% of cells in the PSM at a single timepoint are in M-phase (Horikawa et al., 2006; Kanki and Ho, 1997). Cells in M-phase appear to be distributed homogenously in space (Kanki and Ho, 1997), and, therefore, cell cycles in the PSM may be independent, in which case the total duration $T_{G} + T_{M}$ is approximately equal to $15 / 0.08 = 187.5$ min and $T_{G} = 172.5$ min. Unless otherwise stated, we fix $T_{G} + T_{M} = 187.5$ min.

Modelling compaction-extension

To model the compaction-extension of the PSM (Thomson et al., 2021), we shrink the PSM radius $r_{T}$ , and length $L$ , over time, concurrent with an increase in density $ρ_{0}$ , and a stepwise decrease in cell diameter $d_{c}$ (see Figure 5). This is derived from the observations of Thomson et al., 2021, who observed changes in the PSM height, width, length, density, and cell diameter, during the latter stages of zebrafish somitogenesis. We sought to model this process using the rates of change in these parameters described by Thomson et al., 2021, however, as their published rates are in units of change per somite stage (ss), and the rate of somitogenesis is nonlinear in zebrafish (Schröter et al., 2008), we needed an interpolation of somite number over time in order to describe the according change in rate over time.

To this, we extracted the data for somite number over time from Schröter et al., 2008, using the web app WebPlotDigitizer (Rohatgi, 2022) and fitted a function of the form

s (t) = 6 + \frac{t}{24.7} - a e^{b (t - 300)},

for unknown parameters $a$ , $b$ , using the Julia package LsqFit.jl. We obtained values of $a = 0.5001$ and $b = 0.0049$ using this process. Using this interpolation of somite number over time, we then could change the tissue radius, length, cell diameter, and density, using the values published in Thomson et al., 2021. As this only occurs in the latter stages of somitogenesis (measured only from 16ss onwards), we hold these variables constant before the time at which 16 somites are formed in zebrafish.

The equations describing change in tissue length ( $L$ ), wavefront position ( $x_{a}$ ), and anterior limit of cell addition therefore become:

L (t) = {\begin{cases} L_{0} & s (t) < 16, \\ L_{0} - m_{a} (s (t) - 16) & s (t) \geq 16, \end{cases}

x_{a} (t) = {\begin{cases} x_{a 0} & s (t) < 16, \\ x_{a 0} + m_{a} (s (t) - 16) & s (t) \geq 16, \end{cases}

x_{d} (t) = x_{a} (t) + 100,

where $m_{a} > 0$ denotes the rate at which the PSM shrinks in length, and $L_{0} > 0$ and $x_{a 0} > 0$ the initial values for $L$ and $x_{a}$ before shrinking, respectively. The equation describing the change in radius $r_{T}$ is

r_{T} (t) = {\begin{cases} r_{0} & s (t) < 16, \\ r_{0} - m_{r} (s (t) - 16) & s (t) \geq 16, \end{cases}

where $m_{r} > 0$ denotes the rate at which the PSM shrinks in radius, and $r_{0} > 0$ the initial value for $r_{T}$ before shrinking. In their work Thomson et al., 2021, reported that the PSM shrinks in height more rapidly at the PSM posterior than at the anterior, however, the rate of shrinkage of the PSM in width (along the medio-lateral axis) was the same across the anterior and posterior halves of the PSM (Thomson et al., 2021). As our model assumes the tissue is as wide as it is tall, we take only the rate of shrinkage in the medio-lateral axis in our value for $m_{r}$ , and calculate $m_{r}$ by taking the mean of the two (similar) rates published by Thomson et al., 2021, for the anterior and posterior halves of the PSM.

The functions describing change in tissue density $ρ$ and cell diameter $d_{c}$ are similarly

ρ (t) = {\begin{cases} ρ_{t 0} & s (t) < 16, \\ ρ_{t 0} + m_{ρ} (s (t) - 16) & s (t) \geq 16, \end{cases}

and

d_{c} (t) = {\begin{cases} d_{c 0} & s (t) < 16, \\ d_{c 0} - m_{d_{c}} (s (t) - 16) & s (t) \geq 16, \end{cases}

where $m_{ρ} > 0$ and $m_{d_{c}} > 0$ represent the rates of change in density and cell diameter, respectively, and $ρ_{t 0} > 0$ and $d_{c 0} > 0$ the initial values for density and cell diameter, respectively, before shrinking. In their work Thomson et al., 2021, only report the cell diameter for 16ss and 26ss zebrafish embryos, but we choose to shrink the diameter of cells at the rate they describe, measured between these two timepoints. Our results are insensitive to choice of a function for $d_{c} (t)$ where cells continue shrinking in diameter at the same rate past 26ss, or not (see Figure 5—figure supplement 1).

Initial conditions

To initialise the simulation, $N$ cells are generated with random positions

x_{i} = (\begin{matrix} x_{i} \\ r_{T} + r_{i} \cos q_{i} \\ Z_{c} + r_{i} \sin q_{i} \end{matrix}),

and a further $N$ with positions

x_{i} = (\begin{matrix} x_{i} \\ r_{T} + 2 R + r_{i} \cos q_{i} \\ Z_{c} + r_{i} \sin q_{i} \end{matrix}),

and then a further $N_{T B}$ with positions

x_{i} = (\begin{matrix} X_{c} + R \cos p_{i} + r_{i} \cos p_{i} \cos q_{i} \\ Y_{c} + R \sin p_{i} + r_{i} \sin p_{i} \cos q_{i} \\ Z_{c} + r_{i} \sin q_{i} \end{matrix}),

for $x_{i} \in [x_{a}, X_{c}]$ , $p_{i} \in [- π / 2, π / 2]$ , $q_{i} \in [0, 2 π]$ , and $r_{i} \in [0, r_{T}]$ , where $N = ⌊ ρ_{0} π r_{T}^{2} X_{c} ⌋$ and $N_{T B} = ⌊ ρ_{0} π^{2} r_{T}^{2} R ⌋$ . As this generates more cells near the tissue mid-line than at the periphery, in order to homogenise local density the tissue is ‘relaxed’ for 10 min, whereby $x_{i}$ evolves according to

\frac{d x_{i}}{d t} = v_{0} (x_{i}) n_{i} (t) + \sum_{j = 1, j \neq i}^{N} F (x_{i}, x_{j}) + F_{b} (x_{i}),

and $n_{i} (t)$ evolves according to Equation 4. This evenly distributes cells throughout the tissue.

When cells are initialised prior to relaxation, $n_{i} (t)$ is assigned for each cell at random, i.e., $n_{i} (0) = (\cos φ_{i} \sin ϕ_{i}, \sin φ_{i} \sin ϕ_{i}, \cos ϕ_{i})^{T}$ for random $φ_{i} \in [0, 2 π)$ , $ϕ_{i} \in [0, π]$ . The segmentation clock phase $θ_{i}$ is also assigned as a constant value $θ_{i} (0) = 3 π / 2$ . For simulations with a coupling delay, we initialise clock phase with $θ_{i} (t \leq 0) = 3 π / 2$ .

Parameter values

Parameter values, and their source, are shown in Table 1.

Table 1

Parameter values used for simulations, unless otherwise stated.

Symbol	Value	Units	Source(s)	Symbol	Value	Units	Source(s)
$X_{c}$	300	µm	Uriu et al., 2021	$L$	$X_{c} - x_{a} + R + r_{T}$	µm	Uriu et al., 2021
$Y_{c}$	85	µm	Uriu et al., 2021	$x_{a}$	0	µm	Uriu et al., 2021
$Z_{c}$	25	µm	Uriu et al., 2021	$x_{d}$	100	µm	Uriu et al., 2021
$R$	60	µm	Uriu et al., 2021	$x_{q}$	0.7	–	Uriu et al., 2021
$r_{T}$	25	µm	Uriu et al., 2021	$d_{c}$	11	µm	Uriu et al., 2021
$v_{a}$	1.67	µm⋅min^–1	Uriu et al., 2021	$D$	0.026	min^-1	Uriu et al., 2021
$v_{p}$	3	µm⋅min^–1	Uriu et al., 2021	µ	8.71	µm⋅min^–1	Uriu et al., 2021
$X_{v}$	0.4	–	Uriu et al., 2021	$μ_{b}$	20	µm⋅min^–1	Uriu et al., 2021
$h$	3	–	Uriu et al., 2021	$r_{b}$	1	µm	Uriu et al., 2021
$v_{s}$	1	µm⋅min^–1	Uriu et al., 2021	$ρ_{0}$	0.0015	µm^–3	Uriu et al., 2021
$L_{0}$	325	µm	Thomson et al., 2021	$T_{G}$	187.5–15	min	Kanki and Ho, 1997, Horikawa et al., 2006
$m_{a}$	14.1	µm ⋅somite^–1	Thomson et al., 2021	$T_{M}$	15	min	Horikawa et al., 2006
$r_{0}$	27.6	µm	Thomson et al., 2021	$κ$	0.07	min^–1	Riedel-Kruse et al., 2007
$m_{r}$	1.1625	µm ⋅somite^–1	Thomson et al., 2021	$ω_{0}$	0.2094	min^–1	Soroldoni et al., 2014, Uriu et al., 2021
$d_{c 0}$	9.2	µm	Thomson et al., 2021	$σ$	0.66	–	Uriu et al., 2021
$m_{d_{c 0}}$	0.2	µm ⋅somite^–1	Thomson et al., 2021	$k$	3.07	–	Uriu et al., 2021
$ρ_{t 0}$	0.002123	µm^–3	Thomson et al., 2021	$D θ$	0.0013	min^–1	Riedel-Kruse et al., 2007
$m_{ρ}$	0.000121	µm ⋅somite^–1	Thomson et al., 2021	$τ$	21	min	Herrgen et al., 2010

Implementation

Equations for cell position (Equation 1) and phase (Equation 23) are solved from the initial conditions outlined above across a 1000 min time span using an forward Euler scheme in Julia v1.8.2, using a time step $d t$ = 0.01 min. Code for simulation and analysis can be found at https://github.com/jewh/ModularityPSMClock (copy archived at Hammond, 2024).

Analysis and metrics

Here, we measure clock dynamics using two metrics. The first is the synchrony of the segmentation clock across cells, denoted $r$ , calculated by

r e^{i ψ} = \frac{1}{N} \sum_{j = 1}^{N} e^{i θ_{j}},

which is sometimes referred to as the Kuramoto phase order parameter. When $r \approx 0$ , oscillations are asynchronous, and when $r \approx 1$ , oscillations are synchronous.

The second metric is the average (mean) frequency of cells in a given domain, where the frequency $d θ / d t$ is simply the value of Equation 17. In cases where we are interested in changes in frequency upon varying model parameters, we find it convenient to define

Δ \frac{d θ}{d t} = \frac{d θ}{d t} - \frac{1}{d_{c}} \int_{x_{a}}^{x_{a} + d_{c}} ω (x) d x,

in order to highlight changes in frequency.

In order to analyse the output of a simulation in terms of one value, for many analyses we calculate $r$ , $d θ / d t$ , or $Δ d θ / d t$ , for the cells within one cell diameter of the anterior (i.e. $x_{a} \leq x_{i} \leq x_{a} + d_{c}$ ). We do this because at the wavefront oscillations arrest and pre-pattern somites, and therefore the frequency and synchrony of oscillations immediately adjacent to the PSM anterior will determine the somite length and precision of pre-patterning. We thus regard the values of $r$ , $d θ / d t$ , and $Δ d θ / d t$ as proxies for phenotype exhibited by the simulation. As discussed above, unless stated otherwise we measure clock dynamics after 1000 min, as we find this sufficient to ensure that the dynamics of the clock have reached a steady state that we deem would be reached in vivo (see Figure 10).

Figure 10

Download asset Open asset

Exemplar dynamics of synchrony and mean anterior frequency.

Top: Trace of synchrony $r$ for a $d_{c}$ -wide domain of cells at the left-hand anterior edge of the pre-somitic mesoderm (PSM) over 1000 min. Data drawn from a single simulation with random cell addition, using parameters as per Uriu et al., 2021. Bottom: Trace of mean frequency $d θ / d t$ for a $d_{c}$ -wide domain of cells at the left-hand anterior edge of the PSM over 1000 min. Data drawn from the same simulation as above. Plotted with a black dashed line is the average intrinsic frequency $ω (x)$ across the domain, calculated using the formula shown.

Data availability

All code is available at https://github.com/jewh/ModularityPSMClock (copy archived at Hammond, 2024).

References

(2021) Mechanical control of tissue shape and morphogenetic flows during vertebrate body axis elongation
Scientific Reports 11:8591.

https://doi.org/10.1038/s41598-021-87672-3
- PubMed
- Google Scholar
(2010) A random cell motility gradient downstream of FGF controls elongation of an amniote embryo
Nature 466:248–252.

https://doi.org/10.1038/nature09151
- Google Scholar
(2014) Restricted expression of cdc25a in the tailbud is essential for formation of the zebrafish posterior body
Genes & Development 28:384–395.

https://doi.org/10.1101/gad.233577.113
- Google Scholar
1. Cooke J
2. Zeeman EC
(1976) A clock and wavefront model for control of the number of repeated structures during animal morphogenesis
Journal of Theoretical Biology 58:455–476.

https://doi.org/10.1016/S0022-5193(76)80131-2
- Google Scholar
(2012) Single-cell-resolution imaging of the impact of notch signaling and mitosis on segmentation clock dynamics
Developmental Cell 23:995–1005.

https://doi.org/10.1016/j.devcel.2012.09.009
- Google Scholar
(2007) Setting the tempo in development: an investigation of the zebrafish somite clock mechanism
PLOS Biology 5:e150.

https://doi.org/10.1371/journal.pbio.0050150
- PubMed
- Google Scholar
(2008) Control of segment number in vertebrate embryos
Nature 454:335–339.

https://doi.org/10.1038/nature07020
- Google Scholar
1. Gomez C
2. Pourquié O
(2009) Developmental control of segment numbers in vertebrates
Journal of Experimental Zoology. Part B, Molecular and Developmental Evolution 312:533–544.

https://doi.org/10.1002/jez.b.21305
- PubMed
- Google Scholar
Software
1. Hammond J
(2024) ModularityPSMClock, version swh:1:rev:9f117989683029027e75c507d2bab88dc2f297d8
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:40ab17ba13fc80efd545f24e841f214c8fabbfa6;origin=https://github.com/jewh/ModularityPSMClock;visit=swh:1:snp:83ee8d001d4d37c297cd1e03d7559c1629964b3b;anchor=swh:1:rev:9f117989683029027e75c507d2bab88dc2f297d8
1. Harima Y
2. Takashima Y
3. Ueda Y
4. Ohtsuka T
5. Kageyama R
(2013) Accelerating the tempo of the segmentation clock by reducing the number of introns in the Hes7 gene
Cell Reports 3:1–7.

https://doi.org/10.1016/j.celrep.2012.11.012
- Google Scholar
1. Herrgen L
2. Ares S
3. Morelli LG
4. Schröter C
5. Jülicher F
6. Oates AC
(2010) Intercellular coupling regulates the period of the segmentation clock
Current Biology 20:1244–1253.

https://doi.org/10.1016/j.cub.2010.06.034
- Google Scholar
(2006) Noise-resistant and synchronized oscillation of the segmentation clock
Nature 441:719–723.

https://doi.org/10.1038/nature04861
- Google Scholar
1. Kanki JP
2. Ho RK
(1997) The development of the posterior body in zebrafish
Development 124:881–893.

https://doi.org/10.1242/dev.124.4.881
- Google Scholar
Book
1. Kauffman SA
(1993) The Origins of Order: Self-Organization and Selection in Evolution
Oxford University Press.

https://doi.org/10.1093/oso/9780195079517.001.0001
- Google Scholar
1. Krol AJ
2. Roellig D
3. Dequéant M-L
4. Tassy O
5. Glynn E
6. Hattem G
7. Mushegian A
8. Oates AC
9. Pourquié O
(2011) Evolutionary plasticity of segmentation clock networks
Development 138:2783–2792.

https://doi.org/10.1242/dev.063834
- PubMed
- Google Scholar
Book
1. Kuramoto Y
(1984) Chemical turbulence
In: Kuramoto Y, editors. Chemical Oscillations, Waves, and Turbulence. Springer. pp. 111–140.

https://doi.org/10.1007/978-3-642-69689-3_7
- Google Scholar
1. Lawton AK
2. Nandi A
3. Stulberg MJ
4. Dray N
5. Sneddon MW
6. Pontius W
7. Emonet T
8. Holley SA
(2013) Regulated tissue fluidity steers zebrafish body elongation
Development 140:573–582.

https://doi.org/10.1242/dev.090381
- PubMed
- Google Scholar
1. Lewis J
(2003) Autoinhibition with transcriptional delay: A simple mechanism for the zebrafish somitogenesis oscillator
Current Biology 13:1398–1408.

https://doi.org/10.1016/s0960-9822(03)00534-7
- PubMed
- Google Scholar
(2016) Faster embryonic segmentation through elevated Delta-Notch signalling
Nature Communications 7:11861.

https://doi.org/10.1038/ncomms11861
- PubMed
- Google Scholar
(2018) Segmentation of the zebrafish axial skeleton relies on notochord sheath cells and not on the segmentation clock
eLife 7:e33843.

https://doi.org/10.7554/eLife.33843
- PubMed
- Google Scholar
1. Mallarino R
2. Grant PR
3. Grant BR
4. Herrel A
5. Kuo WP
6. Abzhanov A
(2011) Two developmental modules establish 3D beak-shape variation in Darwin’s finches
PNAS 108:4057–4062.

https://doi.org/10.1073/pnas.1011480108
- PubMed
- Google Scholar
(2007) Priming, initiation and synchronization of the segmentation clock by deltaD and deltaC
Nature Cell Biology 9:523–530.

https://doi.org/10.1038/ncb1578
- Google Scholar
Preprint
1. Michaut A
2. Mongera A
3. Gupta A
4. Serra M
5. Rigoni P
6. Lee JG
7. Duarte F
8. Hall AR
9. Mahadevan L
10. Guevorkian K
11. Pourquié O
(2022) Activity-driven extracellular volume expansion drives vertebrate axis elongation
bioRxiv.

https://doi.org/10.1101/2022.06.27.497799
- Google Scholar
1. Mongera A
2. Rowghanian P
3. Gustafson HJ
4. Shelton E
5. Kealhofer DA
6. Carn EK
7. Serwane F
8. Lucio AA
9. Giammona J
10. Campàs O
(2018) A fluid-to-solid jamming transition underlies vertebrate body axis elongation
Nature 561:401–405.

https://doi.org/10.1038/s41586-018-0479-2
- PubMed
- Google Scholar
1. Monteiro A
2. Prijs J
3. Bax M
4. Hakkaart T
5. Brakefield PM
(2003) Mutants highlight the modular control of butterfly eyespot patterns
Evolution & Development 5:180–187.

https://doi.org/10.1046/j.1525-142X.2003.03029.x
- Google Scholar
1. Morelli LG
2. Ares S
3. Herrgen L
4. Schröter C
5. Jülicher F
6. Oates AC
(2009) Delayed coupling theory of vertebrate segmentation
HFSP Journal 3:55–66.

https://doi.org/10.2976/1.3027088
- PubMed
- Google Scholar
(2002) Segmental relationship between somites and vertebral column in zebrafish
Development 129:3851–3860.

https://doi.org/10.1242/dev.129.16.3851
- PubMed
- Google Scholar
Book
(2013) Modelling oscillator synchronisation during vertebrate axis segmentation
In: Murray PJ, editors. Pattern Formation in Morphogenesis. Springer. pp. 95–105.

https://doi.org/10.1007/978-3-642-20164-6_9
- Google Scholar
(2019) Cell cycle regulation of oscillations yields coupling of growth and form in a computational model of the presomitic mesoderm
Journal of Theoretical Biology 481:75–83.

https://doi.org/10.1016/j.jtbi.2019.05.006
- Google Scholar
(2022) Left-right symmetry of zebrafish embryos requires somite surface tension
Nature 605:516–521.

https://doi.org/10.1038/s41586-022-04646-9
- PubMed
- Google Scholar
(2012) Patterning embryos with oscillations: structure, function and dynamics of the vertebrate segmentation clock
Development 139:625–639.

https://doi.org/10.1242/dev.063735
- PubMed
- Google Scholar
(1997) Avian hairy gene expression identifies a molecular clock linked to vertebrate segmentation and somitogenesis
Cell 91:639–648.

https://doi.org/10.1016/S0092-8674(00)80451-1
- Google Scholar
Book
1. Raff RA
(1996) The Shape of Life – Genes, Development and Evolution of Animal Forms
University of Chicago Press.

https://doi.org/10.7208/chicago/9780226256573.001.0001
- Google Scholar
(1998) Somite number and vertebrate evolution
Development 125:151–160.

https://doi.org/10.1242/dev.125.2.151
- Google Scholar
(2007) Synchrony dynamics during initiation, failure, and rescue of the segmentation clock
Science 317:1911–1915.

https://doi.org/10.1126/science.1142538
- Google Scholar
Software
1. Rohatgi A
(2022) Webplotdigitizer
Version 4.6.

https://automeris.io/WebPlotDigitizer
(2008) Dynamics of zebrafish somitogenesis
Developmental Dynamics 237:545–553.

https://doi.org/10.1002/dvdy.21458
- PubMed
- Google Scholar
1. Schröter C
2. Oates AC
(2010) Segment number and axial identity in a segmentation clock period mutant
Current Biology 20:1254–1258.

https://doi.org/10.1016/j.cub.2010.05.071
- Google Scholar
1. Schröter C
2. Ares S
3. Morelli LG
4. Isakova A
5. Hens K
6. Soroldoni D
7. Gajewski M
8. Jülicher F
9. Maerkl SJ
10. Deplancke B
11. Oates AC
(2012) Topology and dynamics of the zebrafish segmentation clock core circuit
PLOS Biology 10:e1001364.

https://doi.org/10.1371/journal.pbio.1001364
- PubMed
- Google Scholar
1. Simsek MF
2. Özbudak EM
(2018) Spatial fold change of fgf signaling encodes positional information for segmental determination in zebrafish
Cell Reports 24:66–78.

https://doi.org/10.1016/j.celrep.2018.06.023
- PubMed
- Google Scholar
(2023) Periodic inhibition of Erk activity drives sequential somite segmentation
Nature 613:153–159.

https://doi.org/10.1038/s41586-022-05527-x
- PubMed
- Google Scholar
(2014) A Doppler effect in embryonic pattern formation
Science 345:222–225.

https://doi.org/10.1126/science.1253089
- Google Scholar
1. Steventon B
2. Duarte F
3. Lagadec R
4. Mazan S
5. Nicolas JF
6. Hirsinger E
(2016) Species tailoured contribution of volumetric growth and tissue convergence to posterior body elongation in vertebrates
Development 143:1732–1741.

https://doi.org/10.1242/dev.126375
- Google Scholar
(2021) The zebrafish presomitic mesoderm elongates through compaction-extension
Cells & Development 168:203748.

https://doi.org/10.1016/j.cdev.2021.203748
- PubMed
- Google Scholar
(2010) Random cell movement promotes synchronization of the segmentation clock
PNAS 107:4979–4984.

https://doi.org/10.1073/pnas.0907122107
- PubMed
- Google Scholar
1. Uriu K
2. Ares S
3. Oates AC
4. Morelli LG
(2013) Dynamics of mobile coupled phase oscillators
Physical Review E 87:032911.

https://doi.org/10.1103/PhysRevE.87.032911
- Google Scholar
1. Uriu K
2. Morelli LG
(2014) Collective cell movement promotes synchronization of coupled genetic oscillators
Biophysical Journal 107:514–526.

https://doi.org/10.1016/j.bpj.2014.06.011
- PubMed
- Google Scholar
1. Uriu K
2. Bhavna R
3. Oates AC
4. Morelli LG
(2017) A framework for quantification and physical modeling of cell mixing applied to oscillator synchronization in vertebrate somitogenesis
Biology Open 6:1235–1244.

https://doi.org/10.1242/bio.025148
- PubMed
- Google Scholar
1. Uriu K
2. Liao B-K
3. Oates AC
4. Morelli LG
(2021) From local resynchronization to global pattern recovery in the zebrafish segmentation clock
eLife 10:e61358.

https://doi.org/10.7554/eLife.61358
- PubMed
- Google Scholar
1. Venzin OF
2. Oates AC
(2020) What are you synching about? Emerging complexity of Notch signaling in the segmentation clock
Developmental Biology 460:40–54.

https://doi.org/10.1016/j.ydbio.2019.06.024
- Google Scholar
1. Verd B
2. Monk NA
3. Jaeger J
(2019) Modularity, criticality, and evolvability of a developmental gene regulatory network
eLife 8:e42832.

https://doi.org/10.7554/eLife.42832
- PubMed
- Google Scholar
1. Wagner GP
2. Altenberg L
(1996) Perspective: Complex adaptations and the evolution of evolvability
Evolution; International Journal of Organic Evolution 50:967–976.

https://doi.org/10.1111/j.1558-5646.1996.tb02339.x
- Google Scholar
(2007) The road to modularity
Nature Reviews Genetics 8:921–931.

https://doi.org/10.1038/nrg2267
- Google Scholar
1. Webb AB
2. Lengyel IM
3. Jörg DJ
4. Valentin G
5. Jülicher F
6. Morelli LG
7. Oates AC
(2016) Persistence, period and precision of autonomous cellular oscillators from the zebrafish segmentation clock
eLife 5:e08438.

https://doi.org/10.7554/eLife.08438
- PubMed
- Google Scholar
(2020) Mechanical coupling coordinates the co-elongation of axial and paraxial tissues in avian embryos
Developmental Cell 55:354–366.

https://doi.org/10.1016/j.devcel.2020.08.007
- PubMed
- Google Scholar
(2020) Coupling delay controls synchronized oscillation in the segmentation clock
Nature 580:119–123.

https://doi.org/10.1038/s41586-019-1882-z
- Google Scholar

Article and author information

Author details

James E Hammond

Biology Department, University of Oxford, Oxford, United Kingdom

Contribution
Conceptualization, Data curation, Software, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9603-3410
Ruth E Baker

Mathematical Institute, University of Oxford, Oxford, United Kingdom

Contribution
Supervision, Funding acquisition, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-6304-9333
Berta Verd

Biology Department, University of Oxford, Oxford, United Kingdom

Contribution
Conceptualization, Supervision, Funding acquisition, Project administration, Writing – review and editing

For correspondence
berta.verdfernandez@biology.ox.ac.uk

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-9835-009X

Funding

Natural Environment Research Council (Environmental Research DTP grant NE/S007474/1)

James E Hammond

John Fell Fund, University of Oxford (0009780)

Berta Verd

European Research Council (101163722 StG Counts)

Berta Verd

Simons Foundation (MP-SIP-00001828)

Ruth E Baker

Royal Society (RGS\R1\211324)

Berta Verd

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Koichiro Uriu for sharing his code. We thank Usha Kadiyala for comments regarding estimates for T_G. We thank Ben Steventon for feedback on the project. The authors would like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility (http://dx.doi.org/10.5281/zenodo.22558) in carrying out this work. JH was supported by a Natural Environment Research Council DTP studentship (grant NE/S007474/1). This work was supported by a grant from the Simons Foundation (MP-SIP-00001828, REB).

Version history

Preprint posted: February 13, 2024
Sent for peer review: February 9, 2025
Reviewed Preprint version 1: March 11, 2025
Version of Record published: April 1, 2025

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.106316. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.