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 PSM posterior to the anterior. 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); 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. (2010), 2013); 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). Besides advection, cells are subject to random cell mixing, cell-cell repulsion, and a boundary force confining the cells within a hollow horseshoe-shaped domain. 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 position to keep the tissue at a constant density (figure 1B).

Figure 1.

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), 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 the clock’s synchrony and frequency are robust to these changes except when mitosis is introduced, where we find that synchrony is negatively impacted except 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 synchronous 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 domain (see figures 2A 8) if the average density in this domain falls below ρ₀, 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 domains as well, as when cells were only added to the posterior half-toroid domain we observed a loss in density in the two cylindrical domains (data not shown). To avoid bias we add cells at random positions if the average density falls below ρ₀. 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 domains is maintained by addition of cells with constant phase (θ = 0) on the ventral surface of these two lateral domains (see figures 2A 8). 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 addition of cells at random positions within the tissue with random phase (θ ∈ [0, 2π)).

Figure 2.

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 mins (see methods for further justification). After 1000 mins we measure the synchrony (r) and mean frequency (dθ/dt) 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 mins 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). 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. (2010, 2013); 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₀(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θ/dt, see methods) after 1000 mins. 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 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.

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; for example, 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 ρ₀ and studying model dynamics after 1000 mins we see that increasing cell density marginally increases anterior synchrony (figure 4B), and decreases the variability in frequency at the PSM anterior (figure 4C). 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 mins 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). Note that the results for the DP and DP+LV cases are exactly equal for L = 185 µm, as for x_d = 100 µm cell addition is confined to the dorsal-posterior domain and the two ingression methods are numerically equivalent in the model (see methods). While much of this recovery in synchrony may be due to cell addition occurring along a longer domain (and in so doing ‘diluting’ the effect of cell ingression), we observe similar results when cell addition can only occur along a fixed length of the tissue and the total length of the tissue varies (figure 4 supplementary figure 1). Here, the anterior limit of cell addition x_d is fixed in relation to the length of the PSM, by x_d = L − R − r − 100 (where R represents the major radius of the half-toroid domain and r the radius of the PSM tube, see figure 8), and we observe a recovery in synchrony across all three cell ingression scenarios as the total length of the tissue increases (figure 4 supplementary figure 1).

Figure 4.

Figure 4—figure supplement 1.

Figure 4—figure supplement 2.

However 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-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.

Figure 5—figure supplement 1.

Figure 5—figure supplement 2.

Figure 5—figure supplement 3.

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. For small PSM lengths and widths the model generates voids in the tissue as cells do not exit the toroid domain at a sufficient rate to maintain tissue homogeneity, which we deem to be biologically inaccurate. To exclude such cases we only plot synchrony and frequency up until the time at which one of the N=100 replicate simulations encounters such a gap. Thus it is important to note that the simulations presented here do not cover the whole timespan of zebrafish somitogenesis, and we cannot examine dynamics towards the end of the process.

The anterior synchrony and mean anterior frequency over time for N = 100 simulations of compaction-extension are shown in figures 5C and 5D, respectively. They are compared with simulations that do not shrink the PSM and maintain constant L, r, 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-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-D). We treat the effects observed towards the end of the simulation with caution in light of the irregularities in tissue structure mentioned above. In addition, other 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)), and therefore for the latter stages of somito-genesis the results of our simulations may be inaccurate. However, data for the final stages of zebrafish somitogenesis is 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, but that there may be effects at the terminal stages of somitogenesis undetected by this model.

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 ρ₀ (figure 4B-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), 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⁻¹, and spawning a new cell adjacent to that cell once that its τ ≥ 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 τ ∈ [T_G, T_G + T_M) oscillations cease, i.e. dθ_i/dt = 0. 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-C). Fixing T_G +T_M = 187.5 mins and varying T_M shows a T_M -dependent trend where synchrony decreases with increasing T_M (figure 6E). After T_M = 15 mins the synchrony recovers slightly but remains noisy (figure 6E). No noticeable trend is observed for the frequency (figure 6 supp 1A).

Figure 6.

Figure 6—figure supplement 1.

Figure 6—figure supplement 2.

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_div anterior to which (i.e. for x < x_a +x_div) we halt cell cycle progression by fixing τ = 0 mins in all cells. Performing simulations with x_div = 25µm, 175µm, 325µm, we see that while for x_div = 25µm we observe defects in synchrony, for x_div = 175µm, x_div = 325µm, synchrony is rescued (figure 6D). This suggests that there exists some point along the PSM’s AP axis 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_div 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 mins division can occur up to the tissue anterior without a major defect in synchrony (figure 6 supp 2A). For T_M > 15 mins simulations with division show high synchrony at the anterior although not so high as those with T_M < 15 mins, where the value of x_div such that clock synchrony is preserved is displaced to the posterior (figure 6 supp 2A). However, only for T_M = 15 mins is this value majorly displaced to the posterior, with x_div ≈ 150µm being the anterior-most value of x_div such that synchrony is preserved (r ≈ 1) (figure supp 2A). This result is perhaps intuitive, as the parameters in the model set the clock period to 2π/ω₀ = 30 mins, and a clock arrest lasting T_M = 15 mins 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 supp 2B), 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 supp 2B).

Zebrafish clock coupling and cell mixing confers 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. (2010, 2013, 2017)), and the strength of clock coupling, confer robustness too. Within the present model, the global rate of cell mixing is controlled by the parameter υ_s and the clock coupling strength is denoted by k. 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, 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 (κ, υ_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 υ_s and k, studying the anterior synchrony after 1000 mins with N=100 replicates per parameter pair. This reveals a threshold curve of (κ, υ_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 (k, 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⁻¹, υ_s = 1 µm · min⁻¹) lies very close to the threshold for the DP+LV case of cell ingression in the k 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 υ_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.

Figure 8.

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 the clock Schröter et al. (2008); 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 (υ_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 supp 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 sufficiently 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 our choice of κ represents an over-estimate of the strength of delta-notch signalling and, when we combine this with cell mixing in the model, our results present an over-estimate of robustness. However, the precise value of r at which a somite is correctly patterned is unknown, and indeed there exist processes independent of the segmentation clock that can act to correct against inaccuracies in somite length (Naganathan et al. (2022)), so it is plausible that somitogenesis may be more tolerant of clock asynchrony than we suppose here and that somites are correctly patterned for r < 1. As our results appear less sensitive to cell mixing (υ_s) than κ (figure 7), we suggest that the contribution of cell mixing to an experimental estimate of κ may be minor. Therefore, combined with a tolerance of somitogenesis for values of r < 1, we suggest that our findings are robust to any overestimates of κ incorporating cell mixing.

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)).

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 (Uriu et al. (2021)). 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 to-wards the tissue anterior.

Cells are confined within a horseshoe-shaped domain resembling the PSM. This domain is comprised of two cylindrical sub-domains and one half-toroid subdomain, each representing the lateral sides of the PSM and tailbud, respectively. A schematic can be seen in figure 8. Cells are free to move between each domain and leave the tissue at the anterior at x = x_a.

The centre of the toroid domain 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 8). The major radius of the torus is denoted by R and the minor radius (and the radius of the two cylindrical domains) is denoted by r (see figure 8). The resulting length of the PSM in the x direction, L, is therefore given by L = r + 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-cell repulsion, 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 ∈ ℝ³, 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

where v_d (x_i) represents the advection of the cell towards the PSM anterior, υ₀(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

where the function υ_d (x_i) describes the local strain rate along the AP axis:

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 a random walk on the surface of the unit sphere specified by the following differential equation:

where the orthogonal vectors (for e_z = (0, 0, 1)^T)

together define a plane tangent to the unit sphere, and ξ_ϕ (t), ξ_φ (t) are white Gaussian noise terms that satisfy, for any t, t^′, , and . 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, υ₀(x), increases towards the posterior according to equation (6):

The parameters υ_s, X_υ, 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-cell repulsion 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):

where the magnitude of the cell-cell repulsion force F (x_j − x_i) is given by

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 domains or in the half-toroid tailbud (figure 8). If cell i lies within the cylinders, one can use the cylindrical coordinates for cell i, r_i and q_i (see figure 8) to naturally define the boundary force:

where

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:

One can then naturally define

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

Note that in both definitions of F_b(x_i), the tissue domains 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, ρ₀. If the density in any of these subdomains falls below ρ₀, a new cell is added to that domain. If the density of more than one of the domains falls below ρ₀ 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 domain, with a random phase, i.e. if a cell is added within one of the cylindrical domains, it is given a random x_i ∈ [x_a + x_d, X_c], a random r_i ∈ [0, r], and q_i ∈ [0, 2π). Alternatively, if the cell is added within the half-toroid domain, it is given a random p_i ∈ [−π/2, π/2], and random r_i and q_i, defined as above. The cell is also assigned a random phase θ_i ∈ [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 domain on which cells can be added is restricted. For the half-toroid domain, we restrict cell addition to the dorsal quarter of the domain, i.e. q_i ∈ [π/4, 3π/4], p_i ∈ [−π/2, π/2]. For the cylindrical domains we restrict addition to the bottom quarter of the surface, i.e. q_i ∈ [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 τ ∈ 0, T_G + T_M.

Coupled oscillator model of the segmentation clock

Here only the phase of each cell’s segmentation clock cycle 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 κ:

Oscillations are noisy, with phase noise intensity is a Gaussian noise term with for all t, t′. 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)| ≤ d_c }|.

The intrinsic frequency ω(x_i) decreases toward the PSM anterior, creating a frequency gradient and travelling phase waves that move toward the anterior (Soroldoni et al. (2014)). ω(x_i) is defined as:

where ω₀ corresponds to the frequency of cells at the posterior of the tissue, σ the fold-change in frequency at the anterior compared to the posterior, and κ controls the shape of the gradient.

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/dt = 0 if x_i < x_a.

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 τ ∈ [0, T_G + T_M), where T_G and T_M are the duration of G1-G2 and M phases in minutes respectively, and evolve τ according to

where t_addition 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

where H(x) is the Heaviside step function:

If τ_i ≥ T_G + T_M then we consider the cell cycle complete and a new cell κ is added at a random position of length d_new away from cell i, i.e.

for uniformly distributed random ϕ ∈ [0, π], φ ∈ [0, 2π). Other variables are inherited from cell i, i.e. n_κ = n_i and θ_κ = θ_i. The choice of d_new 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_new 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_new = 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 mins, and approximately 8% of cells in the PSM at a single time point 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 there-fore 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 mins, and T_G = 172.5 mins. Unless otherwise stated, we fix T_G + T_M = 187.5.

Modelling compaction-extension

To model the compaction-extension of the PSM (Thomson et al. (2021)), we shrink the PSM radius r, and length L over time, concurrent with an increase in density ρ₀ and a step-wise 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, 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

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:

where m_a denotes the rate at which the PSM shrinks in length, and L₀ and x_a0 the initial values for L and x_a before shrinking, respectively. The equation describing the change in radius r is

where m_r denotes the rate at which the PSM shrinks in radius, and r₀ the initial value for r 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

and

where m_ρ and represent the rates of change in density and cell diameter, respectively, and ρ_t0 and d_c0 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 supp 1).

Initial conditions

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

and a further N with positions

and then a further N_{T B} with positions

for x_i ∈ [x_a, X_c], p_i ∈ [−π/2, π/2], q_i ∈ [0, 2π], and r_i ∈ [0, r], where N = ⌊ρ₀πr²X_c⌋ and N_{T B} = ⌊ρ₀π²r²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 mins, whereby x_i evolves according to

and n_i(t) evolves according to equation (4). This causes cells to be evenly distributed 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 ∈ [0, 2π), ϕ_i ∈ [0, π]. The segmentation clock phase θ_i is also assigned as a constant value θ_i(0) = 3π/2.

Parameter values

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

Table 1.

Implementation

Equations 1 and 17 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 dt = 0.01 mins. Code for simulation and analysis can be found at https://github.com/jewh/ModularityPSMClock.

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

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

The second metric is the average (mean) frequency of cells in a given domain, where the frequency dθ/dt 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

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θ/dt, or Δdθ/dt, for the cells within one cell’s diameter of the anterior (i.e. x_a ≤ x_i ≤ 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θ/dt, and Δdθ/dt as proxies for phenotype exhibited by the simulation. As discussed above, unless stated otherwise we measure clock dynamics after 1000 mins, 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 9).

Figure 9.

Acknowledgements

We thank Koichiro Uriu for his help in providing code and first establishing the model. We thank Usha Kadiyala for her comments regarding estimates for T_G. We thank Ben Steventon for his feedback on the project. The authors would like to acknowledge the use of the University of Oxford Advanced Research Computing (ARC) facility in carrying out this work. http://dx.doi.org/10.5281/zenodo.22558. This work was supported by a grant from the Simons Foundation (MP-SIP-00001828, REB).