During development, cells are subject to stochastic fluctuations in their positions (i.e. cell-level noise) that can potentially lead to morphological noise (i.e. stochastic differences between morphologies that are expected to be equal, e.g. the right and left sides of bilateral organisms). In this study, we explore new and existing hypotheses on buffering mechanisms against cell-level noise. Many of these hypotheses focus on how the boundaries between territories of gene expression remain regular and well defined, despite cell-level noise and division. We study these hypotheses and how irregular territory boundaries lead to morphological noise. To determine the consistency of the different hypotheses, we use a general computational model of development: EmbryoMaker. EmbryoMaker can implement arbitrary gene networks regulating basic cell behaviors (contraction, adhesion, etc.), signaling and tissue biomechanics. We found that buffering mechanisms based on the orientation of cell divisions cannot lead to regular boundaries but that other buffering mechanisms can (homotypic adhesion, planar contraction, non-dividing boundaries, constant signaling and majority rule hypotheses). We also explore the effects of the shape and size of the territories on morphological noise.

The fine complexity of anatomy is reconstructed in each generation through the process of development (Gilbert and Barresi, 2019). This complexity develops in spite of the stochastic fluctuations that are inherent to the molecular and cellular levels. At the molecular level there is Brownian motion (Voet and Voet, 2010) and at the cellular level there are always minute fluctuations in cell position and movement (i.e. cell level noise) (Angelini et al., 2011; Devany et al., 2021; Sepúlveda et al., 2013).

In this article, morphology is understood as the distribution of cells, cell types and extracellular matrix (ECM) in space. By developmental stability we refer to the reliability of developmental outcomes in spite of the noise inherent at the cellular and molecular levels (Waddington, 1942; Klingenberg, 2003; Graham, 2021). Developmental stability can be measured as the inverse of morphological noise: the morphological dissimilarity between genetically identical individuals developing in the same environment or between the left and right side of bilateral organisms. The latter is also called bilateral asymmetry or fluctuating asymmetry (Van Valen, 1962; Swaddle and Witter, 1997; Savriama et al., 2016).

We refer to territory as a set of cells that express a given gene and have a specific continuous distribution over space (Salazar-Ciudad et al., 2003). As a result of processes such as cell signaling and cell movement, territories can change their shape, and new territories can appear and disappear during development (Gilbert and Barresi, 2019). It is widely recognized that having regular territory boundaries is important for developmental stability because it affects proper cell signaling and cell fate determination (Gilbert and Barresi, 2019; Umetsu and Dahmann, 2010; Dahmann et al., 2011). Similarly, the shape of territory boundaries can affect morphological noise directly because it affects where are cell behaviors regulated (e.g. cell contraction, cell adhesion, etc.) and, thus, morphogenesis.

In principle, cell-level noise and cell division can lead cells from one territory to protrude inside others. As a result, territory boundaries can become sinuous (Gilbert and Barresi, 2019; Umetsu and Dahmann, 2010; Dahmann et al., 2011). By sinuous boundaries we mean boundaries that are not regular, straight or gently curved but, instead, are irregular, with introgression of cells between territories at random locations along the boundary.

Based on some experimental evidence and computational models (Aliee et al., 2012; Calzolari et al., 2014; Canela-Xandri et al., 2011; Dahmann et al., 2011; Landsberg et al., 2009; Monier et al., 2010; Umetsu and Dahmann, 2010), previous research has proposed different hypothetical buffering mechanisms to explain how animals can reduce the sinuosity of territory boundaries in spite of cell-level noise and cell division. In this article, we propose additional buffering mechanisms and we study which of them actually hold when implemented in a general model of animal development: EmbryoMaker. Moreover, we also explore how the size and shape of territories affects the sinuosity they can acquire when they grow by cell division during development.

EmbryoMaker (Fig. 1) is a general model of animal development in the sense that it can simulate any developmental mechanism, i.e. any gene network and the cell behaviors (cell division, cell adhesion, cell contraction, etc.), cell mechanical properties (elasticity, plasticity, etc.) and extracellular signals, such as network regulators (Fig. 2) (Salazar-Ciudad et al., 2003). See the methods of Marin-Riera et al. (2016) and Hagolani et al. (2019) for a description of EmbryoMaker.

Fig. 1.

EmbryoMaker. (A) Types of nodes in the model. (B) Cell adhesion and repulsion. If the distance between two nodes is smaller than their equilibrium distances (EQDs), they repulse each other; if it is larger (but still in the range of their adhesion radii, ADD), they attract each other. (C) Bending forces. Rotational bending force precludes cells from moving up or down the epithelial layer; the radial bending force maintains epithelial cells parallel to each other. (D) Extracellular matrix secretion. A cell (in yellow) secretes a new extracellular matrix node (orange sphere). (E) Cell division. New cells initially appear as smaller cylinders that then grow over time. (F) Epithelial elastic spring. The upper and lower nodes of epithelial cells are connected by an elastic spring. (G) Planar contraction. Cells can change their shape by contraction in one direction and expansion in another within the plane of the epithelium. This polarization is in the direction of the gradient of a growth factor in space. (H) Cell contraction in the apical-basal axis. The equilibrium distance and adhesion radius of the yellow nodes become smaller (or larger) due to changes in the genes expressed in the node; as a result, the epithelium invaginates. (I) Extracellular growth factor diffusion. Some gene products are secreted to the extracellular space and diffuse in it (shown in different colors). (J) Gene regulation. Genes can activate or inhibit the expression of each other. In the example, gene 3 is activated by gene 1 and inhibited by gene 2.

Fig. 1.

EmbryoMaker. (A) Types of nodes in the model. (B) Cell adhesion and repulsion. If the distance between two nodes is smaller than their equilibrium distances (EQDs), they repulse each other; if it is larger (but still in the range of their adhesion radii, ADD), they attract each other. (C) Bending forces. Rotational bending force precludes cells from moving up or down the epithelial layer; the radial bending force maintains epithelial cells parallel to each other. (D) Extracellular matrix secretion. A cell (in yellow) secretes a new extracellular matrix node (orange sphere). (E) Cell division. New cells initially appear as smaller cylinders that then grow over time. (F) Epithelial elastic spring. The upper and lower nodes of epithelial cells are connected by an elastic spring. (G) Planar contraction. Cells can change their shape by contraction in one direction and expansion in another within the plane of the epithelium. This polarization is in the direction of the gradient of a growth factor in space. (H) Cell contraction in the apical-basal axis. The equilibrium distance and adhesion radius of the yellow nodes become smaller (or larger) due to changes in the genes expressed in the node; as a result, the epithelium invaginates. (I) Extracellular growth factor diffusion. Some gene products are secreted to the extracellular space and diffuse in it (shown in different colors). (J) Gene regulation. Genes can activate or inhibit the expression of each other. In the example, gene 3 is activated by gene 1 and inhibited by gene 2.

Fig. 2.

The ensemble approach. The initial condition used in this study is shown on the left. The expression of gene 1 is shown in yellow; expression of gene 2 is shown in orange. Two examples of developmental mechanisms are shown. Blue ellipses represent intracellular transcription factors; orange ellipses represent extracellular growth factors; green arrows indicate activation; red arrows indicate inhibition. The images in the central boxes represent the stable gene expression patterns obtained from the developmental mechanisms. Yellow represents high expression and blue represents low expression. The ‘Morphogenesis’ boxes represent how a specific gene in each developmental mechanism regulates a cell property or behavior [orange ellipses indicate cell contraction (COD)]. The boxes on the right show the morphologies resulting from this regulation. Whereas the morphology arising after mechanism 1 is very symmetrical (i.e. low morphological noise), the morphology arising after developmental mechanism 2 is not (i.e. high morphological noise).

Fig. 2.

The ensemble approach. The initial condition used in this study is shown on the left. The expression of gene 1 is shown in yellow; expression of gene 2 is shown in orange. Two examples of developmental mechanisms are shown. Blue ellipses represent intracellular transcription factors; orange ellipses represent extracellular growth factors; green arrows indicate activation; red arrows indicate inhibition. The images in the central boxes represent the stable gene expression patterns obtained from the developmental mechanisms. Yellow represents high expression and blue represents low expression. The ‘Morphogenesis’ boxes represent how a specific gene in each developmental mechanism regulates a cell property or behavior [orange ellipses indicate cell contraction (COD)]. The boxes on the right show the morphologies resulting from this regulation. Whereas the morphology arising after mechanism 1 is very symmetrical (i.e. low morphological noise), the morphology arising after developmental mechanism 2 is not (i.e. high morphological noise).

In this article, we used EmbryoMaker in two complementary approaches: the focal approach and the ensemble approach. In the focal approach, we ran two sets of simulations: the control simulations and the buffering mechanism simulations. Both started from a simple blastula-like spherical epithelium that contained a territory with initially straight boundaries. In the control simulations, the blastula simply grew by cell division and, as a result, the territory became sinuous over time. The buffering simulations were like the control simulations but implementing each hypothetical buffering mechanisms, one at a time. We then compared the sinuosity of the territories arising from each set of simulations. We repeated both sets of focal simulations from blastulas with different shapes and sizes in the initial territory. Using such a simple developmental setting allowed us to study how territory sinuosity is affected by buffering mechanisms, territory size and territory shape in isolation (i.e. without interference from other aspects of development).

In the ensemble approach, we studied the buffering mechanisms in a developmental setting that is more realistic and, thus, more complex. We considered developmental mechanisms in which, in addition to growth, there was cell signaling and other cell behaviors (e.g. cell adhesion, cell contraction, etc.). We did not simulate the development of any specific organ in any specific species. Instead, we simulated a large number of developmental mechanisms picked at random from the space of possible developmental mechanisms. In other words, we constructed a large ensemble of developmental mechanisms by randomly wiring gene products, cell behaviors, extracellular signals and mechanical properties (Fig. 2 and see Materials and Methods). From this large initial random sample, we selected those developmental mechanisms that were able to develop complex morphologies from a simple blastula-like initial condition (Fig. 2). We then re-ran these developmental mechanisms together with each of the buffering mechanisms (one at a time) to explore how the latter affect the sinuosity of territories and the resulting morphological noise.

Cell-level noise was implemented in each simulation as small random displacements in the position of cells (i.e. small increments or decrements in the x, y and z values of each cell) (see Materials and Methods for details). Cell-level noise is characterized by two parameters: its amplitude, i.e. the size of the random displacements applied to each cell; and its frequency, i.e. the proportion of cells that are subject to noise in each iteration.

The relationship between cell-level noise, territory sinuosity and morphological noise

In the focal approach simulations, the boundary between the territories (the yellow and the blue territories in Fig. 3A) was initially straight but, as development progressed, it became sinuous (Fig. 3A-E). This happened because cells dividing close to a territory boundary gave rise to daughter cells that, by mere chance, were sometimes displaced towards the opposite territory (Fig. 3A-E). Territory sinuosity increased with cell division rate for any values of the amplitude and frequency of the cell-level noise (Figs S1 and S2). Naturally, territory sinuosity increased with the amplitude of cell-level noise. However, the fold increase in sinuosity was relatively modest, even for noise amplitudes as large as of 3.5% of the radius of the cells in the initial condition size (Fig. S2). We also found that territory sinuosity is not affected by parameters related to the numerical precision of the calculations in EmbryoMaker (Fig. S3).

Fig. 3.

Cell division creates sinuous territory boundaries and these lead to morphological noise. (A) A blastula-like initial condition showing the yellow territory of a non-diffusible molecule (epithelial cells are indicated by spheres representing their apical sides). The concentration of this molecule is stable over the simulation and is inherited from mother to daughter cells. (B-C*) As the embryo grows, the sinuosity of the territory increases. Small green lines indicate the cell divisions that have the strongest effect on territory sinuosity. All of them are located at the boundary of the territory. (D,E) Territory sinuosity increases over time. Red arrows indicate daughter cells protruding between territories. The simulation starts with 372 cells and finishes with 1670. Cell division occurs homogeneously over space and time. (F) A blastula-like initial condition showing a gradient of concentration in the expression of an intracellular non-diffusible gene product. (G,H) As the embryo grows, cell divisions and cell level noise blur the gradient. (I) Scatter plot showing a positive linear regression (red line; see Table S6 for values) between territory sinuosity and bilateral asymmetry (in metric units of the model) for developmental mechanisms in the ensemble that only regulate cell contraction (see Materials and Methods for details on the measurements). Each point corresponds to the morphology arising from a different developmental mechanism that runs for 10 time units of the model (roughly 10,000 iterations). See Fig. S4 for the same plot but considering the developmental mechanisms regulating other cell behaviors. (J-O) Examples of how sinuous territories lead to noisy morphologies. (J) Upper left: the initial condition with a yellow territory. Lower right: the sinuous territory resulting after simulating cell division for 10 time units. (K,L) Anterior and dorsal views of the simulation J where the yellow territory regulates cell contraction (COD) for one time unit without cell division. (M) Initial condition that is the same size as that in the lower right of J but with straight boundaries. (N,O) Anterior and dorsal views of the morphology resulting from a simulation where the yellow territory in M regulates cell contraction (COD) for one time unit without cell division. In K,L,N,O, colors show the squared distance to the center of the embryo. (P,Q) Box-plots showing the territory sinuosity (P), measured as the ratio perimeter/area, and the bilateral asymmetry (Q) in the simulations with sinuous margins (K,L) and straight margins (N,O) in 10 simulations performed with different random seeds. Wilcoxon test: (P) P<0.001; W=0; (Q) P<0.001; W=1. Dots show outliers, values that are 1.5× the interquartile range (IQR) above the upper quartile or below the lower quartile. Whiskers show the highest and lowest values (excluding outliers). Boxes show the upper and lower quartiles. The horizontal line within the boxes represents the median.

Fig. 3.

Cell division creates sinuous territory boundaries and these lead to morphological noise. (A) A blastula-like initial condition showing the yellow territory of a non-diffusible molecule (epithelial cells are indicated by spheres representing their apical sides). The concentration of this molecule is stable over the simulation and is inherited from mother to daughter cells. (B-C*) As the embryo grows, the sinuosity of the territory increases. Small green lines indicate the cell divisions that have the strongest effect on territory sinuosity. All of them are located at the boundary of the territory. (D,E) Territory sinuosity increases over time. Red arrows indicate daughter cells protruding between territories. The simulation starts with 372 cells and finishes with 1670. Cell division occurs homogeneously over space and time. (F) A blastula-like initial condition showing a gradient of concentration in the expression of an intracellular non-diffusible gene product. (G,H) As the embryo grows, cell divisions and cell level noise blur the gradient. (I) Scatter plot showing a positive linear regression (red line; see Table S6 for values) between territory sinuosity and bilateral asymmetry (in metric units of the model) for developmental mechanisms in the ensemble that only regulate cell contraction (see Materials and Methods for details on the measurements). Each point corresponds to the morphology arising from a different developmental mechanism that runs for 10 time units of the model (roughly 10,000 iterations). See Fig. S4 for the same plot but considering the developmental mechanisms regulating other cell behaviors. (J-O) Examples of how sinuous territories lead to noisy morphologies. (J) Upper left: the initial condition with a yellow territory. Lower right: the sinuous territory resulting after simulating cell division for 10 time units. (K,L) Anterior and dorsal views of the simulation J where the yellow territory regulates cell contraction (COD) for one time unit without cell division. (M) Initial condition that is the same size as that in the lower right of J but with straight boundaries. (N,O) Anterior and dorsal views of the morphology resulting from a simulation where the yellow territory in M regulates cell contraction (COD) for one time unit without cell division. In K,L,N,O, colors show the squared distance to the center of the embryo. (P,Q) Box-plots showing the territory sinuosity (P), measured as the ratio perimeter/area, and the bilateral asymmetry (Q) in the simulations with sinuous margins (K,L) and straight margins (N,O) in 10 simulations performed with different random seeds. Wilcoxon test: (P) P<0.001; W=0; (Q) P<0.001; W=1. Dots show outliers, values that are 1.5× the interquartile range (IQR) above the upper quartile or below the lower quartile. Whiskers show the highest and lowest values (excluding outliers). Boxes show the upper and lower quartiles. The horizontal line within the boxes represents the median.

We also performed simulations from initial conditions in which gene expression varied gradually within the initial territories (i.e. gradients). In this case, we observed that cell division blurred the gradient over space (Fig. 3F-H). A gradient can be understood as a set of concentric narrow quasi-territories (e.g. the cells with the same color in Fig. 3F-H) and, thus, the blurring of a gradient is simply an increase in the sinuosity of its composing quasi-territories. We call these quasi-territories, instead of just territories, because cells within them not only express the same gene, they also express it at a similar level.

To explore whether these results apply to more-complex developmental settings, we simulated a large ensemble of randomly wired developmental mechanisms (the ensemble approach; see Materials and Methods). Among the 60,000 random developmental mechanisms simulated, we picked those that transformed the simple blastula-like initial condition (Fig. 2) into a qualitatively different morphology (i.e. a morphology that is not a sphere and is not broken). In these developmental mechanisms, the shape and size of territories changed over time due to signaling and morphogenesis.

The morphogenesis observed in the ensemble simulations happened because some genes acquired heterogeneous expression patterns and then regulated cell properties or behaviors that lead to forces and, thus, cell movement. This implies that the shape of the territories of these genes determines the spatial distribution of cell behavior regulation. If these shapes are sinuous, then the activation of cell behaviors over space is also sinuous and, thus, can lead to morphological noise. In fact, we did find a positive correlation between the sinuosity of the territories regulating cell behaviors in each developmental mechanisms and the morphological noise of the resulting morphologies (Fig. 3I and Fig. S4). As usual, morphological noise was measured as bilateral asymmetry and territory sinuosity was measured as the number of cells in the boundary versus the total number of cells (see Materials and Methods). The correlation was especially evident for developmental mechanisms that make extensive use of cell contraction or other cell behaviors that generate cellular forces. In the focal approach, it can also be seen that territory sinuosity leads to morphological noise (Fig. 3J-Q).

Hypotheses on how to minimize territory sinuosity

Several hypothetical buffering mechanisms have been proposed in the literature. These can be labeled as the homotypic adhesion, planar contraction, no-division boundary, oriented cell division, majority rule and constant signaling buffering mechanisms. The first three of these hypotheses have been modeled by previous work using a different mathematical model (Canela-Xandri et al., 2011; Dahmann et al., 2011), whereas the others are new.

To test the consistency of these hypotheses, we simulated again the growth of a blastula-like initial condition but implementing each of the buffering mechanisms, one at a time (focal approach, Fig. 4). In this section, we describe each buffering mechanism and its effects on territory sinuosity. Each buffering mechanism is implemented by adding some cell behavior and some gene regulatory network, each with an associated parameter (see Materials and Methods, and Table S6 for specific examples), to the simulation.

Fig. 4.

Some buffering mechanisms decrease the sinuosity of territories, whereas others do not. Focal approach. (A) Example simulation in which a blastula-like embryo with a yellow territory grows in size by cell division only. Cell division occurs homogeneously in time and space. (B-F) As in A but implementing each of the buffering mechanisms. (G-I) Box plots showing the sinuosity of the territory for different values of the parameter in each hypothesis. The red box represents the control (i.e. no buffering mechanism applied). Each simulation was repeated 10 times with different random seeds. (J-M) Example simulations of growth by cell division from initial conditions with different sizes of the rounded yellow territory. (N-Q) As in J-M, but for elongated yellow territories of different sizes. (R,S) Box plots showing the quantification of the sinuosity of territory boundaries for simulations implementing different hypotheses, and different initial sizes and shapes of the yellow territory. For each hypothesis and territory size, 10 simulations were run, each using a different random seed. Territory area and perimeter is calculated in number of cells (see Materials and Methods). All morphologies are at the same scale with the exception of the miniatures in G,H. All simulations were run for 10 time units (roughly 10,000 iterations) and went from 372 cells to 7500. Dots show outliers, values that are 1.5× the interquartile range (IQR) above the upper quartile or below the lower quartile. Whiskers show the highest and lowest values (excluding outliers). Boxes show the upper and lower quartiles. The horizontal line within the boxes represents the median.

Fig. 4.

Some buffering mechanisms decrease the sinuosity of territories, whereas others do not. Focal approach. (A) Example simulation in which a blastula-like embryo with a yellow territory grows in size by cell division only. Cell division occurs homogeneously in time and space. (B-F) As in A but implementing each of the buffering mechanisms. (G-I) Box plots showing the sinuosity of the territory for different values of the parameter in each hypothesis. The red box represents the control (i.e. no buffering mechanism applied). Each simulation was repeated 10 times with different random seeds. (J-M) Example simulations of growth by cell division from initial conditions with different sizes of the rounded yellow territory. (N-Q) As in J-M, but for elongated yellow territories of different sizes. (R,S) Box plots showing the quantification of the sinuosity of territory boundaries for simulations implementing different hypotheses, and different initial sizes and shapes of the yellow territory. For each hypothesis and territory size, 10 simulations were run, each using a different random seed. Territory area and perimeter is calculated in number of cells (see Materials and Methods). All morphologies are at the same scale with the exception of the miniatures in G,H. All simulations were run for 10 time units (roughly 10,000 iterations) and went from 372 cells to 7500. Dots show outliers, values that are 1.5× the interquartile range (IQR) above the upper quartile or below the lower quartile. Whiskers show the highest and lowest values (excluding outliers). Boxes show the upper and lower quartiles. The horizontal line within the boxes represents the median.

Homotypic adhesion hypothesis

If cells in a territory adhere more strongly to each other than to cells belonging to other territories (e.g. because they express cell-adhesion molecules with a high binding affinity for molecules of the same type), then cells from different territories are unlikely to intermix. This adhesion-based cell sorting has long been proposed based on in vitro evidence (Steinberg, 1970; Foty and Steinberg, 2005; Lecuit and Lenne, 2007), and on the lack of intermixing between territories expressing different adhesion molecules in in vivo experiments (García-Bellido et al., 1973; Rodriguez and Basler, 1997; Inoue et al., 2001). In our simulations, homotypic adhesion within territories does indeed reduce territory sinuosity (Fig. 4 and Tables S7-S9). Fig. 4H shows how territory sinuosity quantitatively decreases with the strength of homotypic adhesion.

Planar contraction hypothesis

According to this hypothesis, cells in a territory contract their contact area with cells from other territories, effectively forming a contraction chord along the boundary of the territory (Dahmann et al., 2011). This contraction precludes the intrusion of cells from one territory into another and, thus, leads to straight territory boundaries. These contraction chords have been found, experimentally, in the boundaries of some territories (e.g. in Drosophila wings; Landsberg et al., 2009; Monier et al., 2010; Umetsu et al., 2014) and in rhombomeres of the vertebrate hindbrain (Cayuso et al., 2019).

In our simulations, planar contraction was implemented through the contraction of the sides of cells that are in contact with cells from other territories (see Materials and Methods). As in previous studies (Canela-Xandri et al., 2011; Landsberg et al., 2009), we found that planar cell contraction does indeed reduce the sinuosity of territories (Fig. 4 and Tables S7-S9). Fig. 4G shows how territory sinuosity varies quantitatively with the strength of the planar contraction.

Non-division boundary hypothesis

According to this hypothesis, cells close to territory boundaries divide less than cells far from them (Canela-Xandri et al., 2011; Dahmann et al., 2011). As cell division leads some cells to protrude between territories, precluding cell division in territory boundaries should decrease territory sinuosity. If we preclude the division of the cells in territory boundaries, we observe, as in previous models (Canela-Xandri et al., 2011), low territory sinuosities (Fig. 4 and Tables S7-S9). Fig. 4I shows how territory sinuosity varies quantitatively with division rate in the margins of the territories.

Oriented division hypothesis

If the random displacements suffered by cells after cell division lead to protrusions between territories and, thus, to territory sinuosity, then one can then hypothesize that orienting cell division should reduce territory sinuosity. This could happen by either orienting all cell divisions in the same arbitrary direction or by orienting cell divisions parallel to the boundaries of territories. We found that neither of these hypotheses holds (Fig. 4 and Table S9). The margins of the territories are equally sinuous in the simulations with randomly oriented cell division as with those in the simulations with oriented cell division (Fig. 4).

To better understand why the orientation of cell division fails to affect territory sinuosity (Fig. 5), we tracked the trajectories of individual cells over space and time in the focal simulations. Each small random displacement applied to a cell results in it approaching or distancing itself from neighboring cells. In the former case, pushing forces are elicited between neighbors; in the latter case, pulling forces are elicited between neighbors. As cell-level noise is constantly being applied, cells are constantly pushing and pulling each other in this way. If cell divisions are oriented, daughter cells initially position themselves, or pile up, along the axis of cell division. Over time, however, the constant pushing and pulling between neighboring cells results in cells being displaced away from this axis and towards nearby areas where there is the most free space. Over time, therefore, the trajectory that cells follow in space is the same irrespective of whether cell division is oriented or not. As this also occurs in the territory boundaries, oriented cell division does not preclude the formation of sinuous boundaries.

Fig. 5.

Cell trajectories are not determined by the orientation of division. Focal approach. (A) Cell trajectories in the example control simulation shown in Fig. 4A (i.e. no buffering mechanism). (B) As in A but with cell divisions oriented parallel to the yellow territory boundary. (C) As in A, but with all cell division oriented in the same direction (to the right). (A-C) The colored lines show the trajectories of the first 36 new cells appearing by cell division in the margins of the territory of interest. Red arrows indicate the direction of cell division. Black lines indicate the perimeter of the yellow territory at 0, 3000 and 5000 iterations. (D-F) Box plots showing the angle (in radians) between the movement vectors of the 36 cells and a constant reference vector (1,0). Each box considers the movements occurring within windows of 100 iterations. The plot shows that, in all cases, cells move in all directions, except for the earliest iterations of the oriented cell division simulations.

Fig. 5.

Cell trajectories are not determined by the orientation of division. Focal approach. (A) Cell trajectories in the example control simulation shown in Fig. 4A (i.e. no buffering mechanism). (B) As in A but with cell divisions oriented parallel to the yellow territory boundary. (C) As in A, but with all cell division oriented in the same direction (to the right). (A-C) The colored lines show the trajectories of the first 36 new cells appearing by cell division in the margins of the territory of interest. Red arrows indicate the direction of cell division. Black lines indicate the perimeter of the yellow territory at 0, 3000 and 5000 iterations. (D-F) Box plots showing the angle (in radians) between the movement vectors of the 36 cells and a constant reference vector (1,0). Each box considers the movements occurring within windows of 100 iterations. The plot shows that, in all cases, cells move in all directions, except for the earliest iterations of the oriented cell division simulations.

Majority rule hypothesis

This hypothesis assumes that each territory secretes a different short-range signal, and that cells in territory boundaries express the genes of one territory or the other depending on which of these signals they receive the most. Thus, a cell from a territory A that becomes surrounded by cells from a territory B will receive more B signal than A signal, and would then change its gene expression to match that of territory B. Similar majority rule mechanisms have been proposed to maintain the identity of rhombomeres in the vertebrate hindbrain (Addison et al., 2018; Qiu et al., 2021; Wilkinson, 2018). Our simulations show that this buffering mechanism can indeed decrease the sinuosity of territories (Fig. 6 and Table S10). Fig. 6 explains how this mechanism is implemented and how territory sinuosity quantitatively varies with the buffering mechanism parameters.

Fig. 6.

The majority rule hypothesis. Focal approach. (A) Initial blastula-like condition with a yellow territory. (B) Control simulation with only cell division (homogeneous in time and space). (C) As in B, but implementing the majority rule hypothesis. (D-G) Zoomed screenshot of the simulation showing example cases of cells (marked with asterisks) changing their expression to match that of their neighbors. (H) Boxplot showing the quantitative results of territory sinuosity when modifying the concentration values of the receptor in the majority rule. We used territory sinuosity instead of area/perimeter because the gene expression in the territory is not binary. Each box represents 10 simulations performed with the same parameters but different random seed. The red box shows the control simulations with only cell division and no majority rule. (I) The buffering mechanism used to implement the majority rule in the simulations. Gene 1 encodes a membrane-tethered ligand (yellow) that can bind to receptor 2 in other cells (receptor 2 is assumed to be constitutively expressed in all cells). When ligand 1 and the receptor bind, they inhibit the expression of gene 3 (blue). Gene 3 inhibits gene 1. Gene 3 is activated by gene 4. Gene 4 is expressed homogeneously over the whole embryo. In cells that become surrounded by cells expressing gene 1 (i.e. the ligand), gene 1 signaling becomes so strong that it overcomes the activating effect of gene 4 over gene 3. As a result, gene 3 is inhibited in these cells and gene 1 becomes expressed (cells become yellow), as shown in D-G. Dots show outliers, values that are 1.5× the interquartile range (IQR) above the upper quartile or below the lower quartile. Whiskers show the highest and lowest values (excluding outliers). Boxes show the upper and lower quartiles. The horizontal line within the boxes represents the median.

Fig. 6.

The majority rule hypothesis. Focal approach. (A) Initial blastula-like condition with a yellow territory. (B) Control simulation with only cell division (homogeneous in time and space). (C) As in B, but implementing the majority rule hypothesis. (D-G) Zoomed screenshot of the simulation showing example cases of cells (marked with asterisks) changing their expression to match that of their neighbors. (H) Boxplot showing the quantitative results of territory sinuosity when modifying the concentration values of the receptor in the majority rule. We used territory sinuosity instead of area/perimeter because the gene expression in the territory is not binary. Each box represents 10 simulations performed with the same parameters but different random seed. The red box shows the control simulations with only cell division and no majority rule. (I) The buffering mechanism used to implement the majority rule in the simulations. Gene 1 encodes a membrane-tethered ligand (yellow) that can bind to receptor 2 in other cells (receptor 2 is assumed to be constitutively expressed in all cells). When ligand 1 and the receptor bind, they inhibit the expression of gene 3 (blue). Gene 3 inhibits gene 1. Gene 3 is activated by gene 4. Gene 4 is expressed homogeneously over the whole embryo. In cells that become surrounded by cells expressing gene 1 (i.e. the ligand), gene 1 signaling becomes so strong that it overcomes the activating effect of gene 4 over gene 3. As a result, gene 3 is inhibited in these cells and gene 1 becomes expressed (cells become yellow), as shown in D-G. Dots show outliers, values that are 1.5× the interquartile range (IQR) above the upper quartile or below the lower quartile. Whiskers show the highest and lowest values (excluding outliers). Boxes show the upper and lower quartiles. The horizontal line within the boxes represents the median.

Constant signaling hypothesis

In this mechanism there is a territory secreting an extracellular signal and an induced territory receiving it. The induced territory exists only as long as the signal is being received within some range of concentrations. As a result, the induced territory is always at a specific range of distances from the inducing territory and has a shape determined by that of the inducing territory. The size and shape of the induced territory are unaffected by cell division in it (or by the displacement of cells in general). Simply, cells from the induced territory that become pushed too far away from or too close to the inducing territory will stop expressing the induced gene and, thus, will no longer belong to the induced territory. In this way, the sinuosity of the territory can remain small. Fig. 7 shows that this hypothesis holds when implemented in EmbryoMaker (see also Fig. S6). Note that, as we later discuss, this mechanism requires that the shape of the inducing territory does not change much over development. In our simulations, this occurred because no cell division occurred in the territory secreting the extracellular signal.

Fig. 7.

Constant signaling can reduce morphological noise. Focal approach. (A,B) The dynamics of a developmental mechanism with (A) or without (B) constant signaling. With constant signaling (A), the final morphologies have lower morphological noise than without constant signaling. Both simulations start from the same initial conditions used for the other figures. Screenshots from different times during the simulations are shown. In B, first there is cell signaling (due to the gene network shown on the left). After this, there is cell division without cell signaling (t2). After this, gene 4 upregulates cell contraction, thus leading to a change in morphology (t3). In the left panel, cell signaling occurs continuously. Cell division starts at the beginning of the simulation (t0) and then stops (t2). From then, cell signaling continues but cell division does not. Gene 4, however, starts to activate cell contraction (t3). From t1 to t2, the colors represent the level of expression of gene 4. Later, they represent the expression of gene 4 (right) and the squared distance to the center (left).

Fig. 7.

Constant signaling can reduce morphological noise. Focal approach. (A,B) The dynamics of a developmental mechanism with (A) or without (B) constant signaling. With constant signaling (A), the final morphologies have lower morphological noise than without constant signaling. Both simulations start from the same initial conditions used for the other figures. Screenshots from different times during the simulations are shown. In B, first there is cell signaling (due to the gene network shown on the left). After this, there is cell division without cell signaling (t2). After this, gene 4 upregulates cell contraction, thus leading to a change in morphology (t3). In the left panel, cell signaling occurs continuously. Cell division starts at the beginning of the simulation (t0) and then stops (t2). From then, cell signaling continues but cell division does not. Gene 4, however, starts to activate cell contraction (t3). From t1 to t2, the colors represent the level of expression of gene 4. Later, they represent the expression of gene 4 (right) and the squared distance to the center (left).

The effect of territory shape and size on territory sinuosity

If cell division increases morphological noise by distorting the boundaries of territories, narrow or small territories should be more easily distorted by cell division. This is because small territories and territories with narrow shapes have, due to mere geometry, relatively larger boundaries (i.e. perimeters) with respect to their area.

To better characterize these effects, we repeated the simulations where a blastula-like initial condition grows by cell division, but using initial territories with different sizes and shapes (Fig. 4J-S). These simulations show that territories that are initially small or narrow develop very sinuous boundaries as they grow. Conversely, territories that are initially large or wide, grow to be larger but with fewer sinuous boundaries. The fold increase in cell number is the same in all simulations, so the crucial factor is not the amount of growth per se, but the size and shape of the territories before growth. We also ran these simulations under different buffering mechanisms (Fig. 4R,S and Fig. S5).

Very large territories are associated with morphological noise

From the previous section, it may seem that large territories ensure developmental stability. However, a previous study with EmbryoMaker shows that beyond certain size, large morphological noise arises if large territories occur during development (Hagolani et al., 2019). This effect is most notable when cells in large territories are dividing or contracting.

To acquire a better grasp of which range of territory sizes are the most compatible with developmental stability, we ran a set of simulations from initial conditions that have territories of different sizes, i.e. from embryos completely partitioned into different number of territories (Fig. 8). Apical contraction or expansion was then promoted in alternative contiguous territories (i.e. a territory with apical contraction is always surrounded by territories with apical expansion, and vice versa). These simulations, thus, represent an idealized simple development (i.e. no cell signaling, no territories with complex shapes) that is easy no analyze and yet leads to some non-trivial morphogenesis.

Fig. 8.

Medium-sized territories are associated with lower morphological noise. Focal approach. (A-E) Four different simulations starting from the same blastula-like initial condition that have been divided into different number of territories (4, 8, 16 and 32) and, thus, contain territories of different sizes. In each simulation, the territories in yellow regulate cell contraction positively, meaning that these cells will increase their apical-side area while decreasing the basal-side area. The opposite happens in the blue territories. Examples of resulting morphologies with 6000 cells are shown in B-E. The dashed red line represents the plane of symmetry. (A) Plot showing the relationship between bilateral asymmetry and the number of cells for each simulation. Ten replicas of each simulation were performed with different random seeds. The dots show the average bilateral asymmetry and the vertical lines are 2×s.d.

Fig. 8.

Medium-sized territories are associated with lower morphological noise. Focal approach. (A-E) Four different simulations starting from the same blastula-like initial condition that have been divided into different number of territories (4, 8, 16 and 32) and, thus, contain territories of different sizes. In each simulation, the territories in yellow regulate cell contraction positively, meaning that these cells will increase their apical-side area while decreasing the basal-side area. The opposite happens in the blue territories. Examples of resulting morphologies with 6000 cells are shown in B-E. The dashed red line represents the plane of symmetry. (A) Plot showing the relationship between bilateral asymmetry and the number of cells for each simulation. Ten replicas of each simulation were performed with different random seeds. The dots show the average bilateral asymmetry and the vertical lines are 2×s.d.

Our results (Fig. 8) show that morphological noise is higher in the simulations with relatively large territories and with many small territories than in the simulations with medium-size territories. All simulations were started from the same number of cells and finished with the same number of cells, so our results are not due to embryo size or relative growth itself, but to the territory sizes.

Buffering mechanisms acting in complex developmental settings

The simulations of buffering mechanisms in previous work (Canela-Xandri et al., 2011; Dahmann et al., 2011), and in the present article so far, are restricted to very simple developmental settings (i.e. flat epithelia, spheric blastulae, etc.). The question, then, is whether the studied buffering mechanisms would also increase developmental stability in more biologically realistic developmental settings. To explore this question, we implemented each buffering mechanism, one at a time, in the 2230 developmental mechanisms in the ensemble that lead to the development of non-trivial morphologies from the simple blastula-like initial condition. For the parameter values of each buffering mechanism (Table S5), we chose those that significantly reduced morphological noise in the focal approach (Fig. 3), but did not produce any appreciable morphological change.

We found that implementing the buffering mechanisms with such parameter values only reduced morphological noise in some of the developmental mechanisms in the ensemble (when compared with the developmental mechanisms without these buffering mechanisms). This happened for 45.34% of the developmental mechanisms when implementing homotypic adhesion, for 28.58% when implementing planar cell contraction and for 61.55% when implementing the non-division boundary. We then explored whether a reduction in morphological noise could be achieved by providing other values (i.e. fine-tuning) to the parameters of each buffering mechanisms. Fig. 9 shows example developmental mechanisms from the ensemble where morphological noise was further reduced through fine-tuning of the parameters in each buffering mechanism. We found that buffering mechanisms could work for many developmental mechanisms but the parameter values for which this happened were different for each combination of developmental mechanisms and buffering mechanisms. This diversity is perhaps not surprising as each developmental mechanism has its own dynamics (e.g. different cell division rates in different parts of the embryo, different cell behaviors involved, etc.) and, thus, each is differently affected by each buffering mechanism and by the specific parameter values thereof.

Fig. 9.

Different buffering mechanisms can decrease morphological noise in complex developmental settings if they are fine-tuned. Six example simulations of developmental mechanisms (from the ensemble) in which we have applied a specific buffering mechanism or not (control red boxes). The parameter values of the examples are in Table S5. The box plots show the bilateral asymmetry obtained in simulations for different values of the parameters of the buffering mechanism. Each box represents 10 simulations with the same parameters but changing the random seed. The results of the non-parametric statistical tests comparing each box with the control and the non-parametric linear regressions are in Tables S11 and S12. Dots show outliers, values that are 1.5× the interquartile range (IQR) above the upper quartile or below the lower quartile. Whiskers show the highest and lowest values (excluding outliers). Boxes show the upper and lower quartiles. The horizontal line within the boxes represents the median.

Fig. 9.

Different buffering mechanisms can decrease morphological noise in complex developmental settings if they are fine-tuned. Six example simulations of developmental mechanisms (from the ensemble) in which we have applied a specific buffering mechanism or not (control red boxes). The parameter values of the examples are in Table S5. The box plots show the bilateral asymmetry obtained in simulations for different values of the parameters of the buffering mechanism. Each box represents 10 simulations with the same parameters but changing the random seed. The results of the non-parametric statistical tests comparing each box with the control and the non-parametric linear regressions are in Tables S11 and S12. Dots show outliers, values that are 1.5× the interquartile range (IQR) above the upper quartile or below the lower quartile. Whiskers show the highest and lowest values (excluding outliers). Boxes show the upper and lower quartiles. The horizontal line within the boxes represents the median.

In this article, we have studied how cell-level noise and cell division can preclude the formation of straight boundaries between adjacent territories of gene expression (Fig. 3A-E) and lead to the blurring of gradients (Fig. 3F,G). We have also described how these distortions on territory boundaries correlate with morphological noise at the level of the whole embryo (Fig. 3I-O).

Most of the proposed hypotheses on how to decrease territory sinuosity are internally consistent and can, indeed, decrease territory sinuosity in our model. These includes the new hypotheses we proposed, the majority rule and constant signaling. The only exceptions are the hypotheses based on orienting cell division. Of all the studied buffering mechanisms, planar cell contraction seems to be the one that attains lower territory sinuosity (Fig. 4G-I,R,S).

Beyond these mechanistic insights, we also found that the shape and size that territories have at any given moment greatly influence the sinuosity they will exhibit after growing by cell division. Small and narrow territories growing by cell division lead to very sinuous territories. Large territories also lead to morphological noise, so it is territories of intermediate size that are the least prone to morphological noise (Fig. 8).

For the sake of simplicity, we have not included mesenchymal cells in our study. Except for the bending forces (see methods), mesenchymal and epithelial cells behave in very similar ways in our model. Thus, there is no reason to expect that our results would not also apply to mesenchymal tissues. However, there is at least one buffering mechanism that is more likely to work in the mesenchyme than in epithelia: the accumulation of extracellular matrix between territories (Dahmann et al., 2011). This mechanism can reduce morphological noise by creating a physical barrier that precludes territories from mixing, as it has been observed in the segmentation of zebrafish (Koshida et al., 2005).

Gradients and buffering mechanisms

Only some of the studied buffering mechanisms can be applied to gradients. As gradients can be conceptualized as being composed of thin quasi-territories of progressively high concentration, their sinuosity should be reducible, in principle, by buffering mechanisms affecting the sinuosity of their composing quasi-territories. However, most quasi-territory are likely to be thin and this, we think, may preclude most buffering mechanisms from effectively acting on them. For example, homotypic adhesion within quasi-territories would lead to cells from each quasi-territory to surround each other in order to maximize their contact area. This will lead the gradient to stop being composed of thin concentric quasi-territories and become a collection of juxtaposed round and thick territories. Similarly, the non-diving boundary buffering mechanism requires one or several layers of non-dividing cells in the boundary and, thus, cannot effectively act on thin quasi-territories made of very few cell layers.

In principle, the majority rule mechanism could act on narrow quasi-territories. However, that would require a specific majority-rule gene network for each quasi-territory. Similarly, the planar contraction buffering mechanisms would require that cells signal to their neighbors in order to determine which of their sides are in contact with cells from other quasi-territories and should, thus, contract. In both cases, many different regulatory interactions (i.e. a specific gene network per quasi-territory) would be required just to maintain the shape of gradients.

These shortcomings do not apply to the constant signaling buffering mechanism. This mechanism implies the formation of at least one extracellular signal gradient. This gradient can be interpreted by cells to lead to discrete territories or to lead directly to gradient territories. Thus, this mechanism has the advantage that it applies to discrete territories as well as to gradients. One shortcoming of this buffering mechanism is that it can act only on the gradients produced by itself. In addition, it requires that the shape of the inducing territory itself does not change. This can be achieved by any buffering mechanism (including constant signaling from other territories) or by precluding cell division in the inducing territories, which occurs in many signaling centers (Jernvall et al., 1998). For all these reasons, we think that the constant signaling buffering mechanism is especially well suited to stabilizing gradients.

Buffering mechanisms in complex developmental settings

Noise buffering mechanisms are usually studied in relatively simple settings such as a flat epithelium (Canela-Xandri et al., 2011; Dahmann et al., 2011; Landsberg et al., 2009). In this article, we found that in more-complex settings, where morphology changes in complex ways and there are several cell behaviors involved, buffering mechanisms can still decrease morphological noise. However, our simulations suggest that each developmental mechanism can require different values in the parameters of the buffering mechanisms for such a decrease to occur. This possible requirement suggests that buffering mechanisms may not be recruited in all territories or body parts during development. Simply minimizing morphological noise in all the territories of all developmental stages may require a lot of regulation and we do not think that such tight regulation is likely to have arisen in evolution. Instead, we expect that buffering mechanisms may only be deployed and fine-tuned in territories that are especially prone to noise or the regularity of which is especially crucial for later development, as we detail later.

Buffering mechanisms and territory size

We think that buffering mechanisms are much more likely to be found acting on small territories than on intermediate-size territories, as we found that the former become very sinuous when they grow. However, all the buffering mechanisms we have studied do not simply reduce the sinuosity of territories, they also alter their shape. This implies that small territories in embryos may tend to have shapes that are not as sensitive to cell-level noise, e.g. round, or be transformed into these shapes by buffering mechanisms. In other words, there are some territory shapes that may be unlikely to be found in animal development. For example, small narrow territories with complex shapes would be rare, as their shape would easily be distorted by either cell division or buffering mechanisms. Thus, we expect to find this kind of territory only in organs or in whole organisms, where cell division has already ceased, because it occurs in the adults of many species (Ruppert et al., 2004).

Buffering mechanisms and the diversity of animal development

We think our results on territory sizes and shapes may have some implications for the structure of development in different animal groups and their evolution. Different animal groups tend to have different territory sizes. Many invertebrates tend to have small territories during their development (Brozovic et al., 2018; Ma et al., 2021; Özpolat et al., 2021; Salvador-Martínez and Salazar-Ciudad, 2017). This is particularly the case in the so-called type I developers (Davidson, 1991; Davidson et al., 1995; Salazar-Ciudad, 2010). These include animal groups that have small bodies as adults (e.g. nematodes, tardigrades, rotifers and gastrotrichs) or in which there is a free larva stage made of a relatively small number of cells (e.g. trocophora larva of mollusks, annelids and other lophotrochozoan phyla; Nielsen, 2004, 2005). In these animals, the embryo is partitioned into small territories from very early developmental stages (i.e. most cells in the early embryo have different gene expression profiles). This partitioning usually involves short-range cell signaling between very small groups of cells, even between individual pairs of cells (Davidson, 1991). Moreover, the orientation of cell division and cell adhesion tend to be finely regulated, even at the level of single individual cells (Davidson, 1991; Davidson et al., 1995; Salazar-Ciudad, 2010). In the larva and in the adults of many type I developers, functional organs are made of one or a few cells, so it is likely that any noise in the positioning and signaling between cells can lead to the disappearance of a whole organ and, thus, be highly deleterious (Davidson, 1991). Because of this and the small size of territories, we expect that buffering mechanisms are extensively used in type I developers.

Almost the opposite of type I development is found in the development of many vertebrates, in which more growth occurs between (or during) signaling events than it does in type I developers (Davidson, 1991; Salazar-Ciudad, 2010). Thus, territories tend to be larger. In oviparous vertebrates, this is the case even in the earliest developmental stages: oocytes tend to be relatively large and there is not much cell signaling until the late blastula stage, when the embryo is already made of several thousands of cells (Gilbert and Barresi, 2019). Oocytes are already partitioned into a small number of simple intracellular territories (e.g. intracellular regions with specific long-lived mRNAs or proteins shed there by the mother; Gilbert and Barresi, 2019). These early territories are large in absolute size, because oocytes are also relatively large, at least compared with those found in type I development (Gilbert and Raunio, 1997). During cleavage, fast cell divisions occur without cell growth. This transforms the large intracellular territories into equally large cell territories. Owing to the large territory sizes, any possible intermixing of cells between territories is unlikely to have a major effect on the sinuosity of the territory boundaries. In this way, vertebrates have large territories with relatively straight boundaries from the earliest developmental stages. This may prevent the growth of large territories from small territories and the morphological noise that arises from that. Thus, owing to the large sizes that territories tend to have in vertebrates, we expect that they do not use buffering mechanisms as often as type I developers do.

EmbryoMaker

A detailed description of EmbryoMaker can be found in the original publications (Marin-Riera et al., 2016; Hagolani et al., 2019). The latest version of EmbryoMaker, the one used in this article, can be obtained from https://github.com/HugoCanoFernandez/EmbryoMaker. The general parameters used for all the simulations are shown in Tables S1-S4. Simulations were run using the fourth order Runge-Kutta method with a dynamic step size.

EmbryoMaker considers mesenchymal cells, epithelial cells, ECM and embryos made of these. In this study, we chose to assign each epithelial cell two nodes, an apical node and a basal node, and mesenchymal cells and extracellular matrix (ECM) one node (Fig. 1A).

EmbryoMaker calculates how an initial condition (i.e. an initial morphology) changes over time according to a developmental mechanism provided by the user. EmbryoMaker calculates how the variables provided in the initial conditions change over time. These variables include the position of each cell in 3D space (i.e. overall morphology), the expression of each gene in each cell and the mechanical properties of each cell (i.e. its size, its adhesivity, etc.). The calculations are carried out by iterating over time a set of differential equations and logical rules acting on the variables and governed by the developmental mechanism provided by the user (Fig. 1), and the quantitative values of the parameters associated with each developmental mechanism (e.g. how strongly a gene activates the expression of another or promotes cell contraction, etc.).

The equations and rules in EmbryoMaker are simple. They specify, for example, that if two cells become too close they will repulse each other. If two cells are close but not too close, a pulling force arises between them (Fig. 1B). Epithelial cells have additional forces against the bending of the epithelium (Fig. 1C).

A set of equations determines how the concentration of gene products changes as a consequence of regulation by other gene products within a cell (according to the provided developmental mechanism) (Fig. 1J). Some gene products can diffuse over extracellular space and affect gene expression in other cells in different locations in space (Fig. 1I).

Another set of equations determines how the behaviors and mechanical properties of a cell at a given time change, quantitatively, owing to regulation by the genes expressed by the cell. Each cell behavior in EmbryoMaker is implemented as a logical rule or as a variable in the equations that drive mechanical interactions (Fig. 1).

Ensemble

Sixty-thousand gene developmental mechanisms were created by randomly wiring 30 genes. When constructing the developmental mechanisms, each gene had a 50% probability of coding for a growth factor (i.e. extracellularly diffusible signal) and a 50% probability of coding for a non extracellularly diffusible gene product. Every gene had a 25% probability of regulating any other given gene in the network and this regulation could be either activation or inhibition with a probability of 50%. Notice that, in EmbryoMaker, a gene network simply specifies which gene products would regulate each other if they happened to be in the same cell. Whether this will happen, or not, for each given cell depends on the dynamics of each developmental mechanism, i.e. it is not specified from the beginning.

Each of the 60,000 developmental mechanisms were first simulated with no regulation of any cell property or behavior for 120 time units (around 1,200,000 iterations). From that initial set, we selected only the developmental mechanisms leading to at least two genes having heterogeneous expression in space, as long as these genes were expressed in more than 5% of cells and at a concentration higher than 10−8 (with respect to the theoretically maximum gene expression possible in a given developmental mechanism given its parameter values). We refer to these two genes as the pattern genes of a developmental mechanism. We also discarded the developmental mechanisms that formed territories of gene expression that were unstable over time or too similar to those in the initial condition (see Hagolani et al., 2019) for details).

To the resulting 52 signaling developmental mechanisms, we added cell behavior regulation to create complete developmental mechanisms. Each of these latter developmental mechanisms was simulated in two phases. In the first phase, we allowed cell signaling but no regulation of any other cell behavior. Once the patterns of gene expression reached a steady state, we started the second simulation phase. In the second phase, there was no cell signaling but a pattern gene product was allowed to regulate a cell behavior or cell property. All the possible combinations of pattern genes and regulation of cell properties or behaviors were simulated, obtaining a total number of 2320 complete developmental mechanisms and simulations. The intensity of the regulation of each cell behavior or cell property in a developmental mechanism was chosen to be 75% of the maximum intensity (see Table S4). In other words, the intensity of such regulation, eij, for a given cell behavior or node property i in a developmental mechanism j, was chosen to be eij=0.75*pimax gj. Where pimax is the maximal value a cell behavior or node property i could have without breaking the epithelium (see last paragraph of this section and Table S4). gj is the maximal expression that the patterning gene reached in any given cell in the developmental mechanism at the end of phase one. Notice that gj and pimax are not parameters, i.e. their values are a result of model dynamics. The equation for eij was chosen to ensure that each cell behavior or cell property will be regulated at a level likely to lead to morphogenetic changes but unlikely to lead to the embryo breaking (as the regulation of each cell behavior i or the value of any node property i depends on eij*g in each cell, see Marin-Riera et al., 2016).

The second phase was run for 10 time units, which corresponds to roughly 10,000 iterations. Simulations were prematurely stopped if any of the following conditions were met: a cell acquired more than 100 neighbors (i.e. cells close enough to interact); the embryo reached 10,000 cells; 4 h of computation time had passed; or if the total radius of the embryo exceeded one thousand space units (each cell in the initial conditions has a radius equal to 0.14 space units). The final morphologies were manually checked and the ones with a broken epithelium or no morphogenesis were discarded from further study.

The ranges of the parameters controlling the activation of cell properties and behaviors (Table S4) were obtained by performing an independent set of simulations where these cell properties and behaviors were activated homogeneously over space and at the same time the blastula grew by cell division. We tested a wide range of values for these parameters and we considered as biologically relevant the intervals in which morphogenesis occurred without the epithelium breaking.

Cell-level noise

In each iteration, a proportion of the nodes was chosen, which was designated the model parameter MNOI (Marin-Riera et al., 2016), and noise was applied to them. MNOI is what we call the frequency of noise in this article. The noise applied to each node consisted in a random 3D vector added to its position. This vector had three components: two angles and a radius. The angles were chosen at random so that displacements could occur in any direction with equal probability. The radius of each vector was also chosen at random with a uniform distribution from 0 to the amplitude of the noise.

Buffering mechanisms

Homotypic adhesion was implemented simply by changing the binding affinity between cells according to the territories to which they belonged. Between any two cells there is an adhesion force that has a generic component (the cell property pADH; Marin-Riera et al., 2016) and a component that depends on the adhesion molecules that the cells express and the binding affinity of these molecules. The binding affinity of each homotypic adhesion molecule is specified by a parameter bii in EmbryoMaker, where i is the adhesion molecule (Table S5). The adhesion force between two cells is then the generic adhesion plus the product of concentration of the adhesion molecules in each cell and by their binding affinity. In this buffering mechanism, we only consider molecules with homotypic adhesion, so that the adhesion within territories is larger than between territories. Without loss of generality, we set the expression level of each homotypic adhesion molecule to be 1 in all the cells within its territory of expression and 0 elsewhere. This way the only parameter regulating the homotypic adhesion is the binding affinity of each adhesion molecule to itself, bii (Table S5).

The planar contraction hypothesis was implemented by changing the shape of the cells in the margins of the territory. We modified EmbryoMaker to allow cells made of a single node to contract in the plane of the epithelium, becoming ellipses instead of cylinders (Fig. 1G) (see the planar cell contraction section below for details). In the previous version of EmbryoMaker, it was only possible to do this for cells made of many nodes. This planar cell contraction is controlled by two genes, one that determines the strength of the contraction and another one that determines its direction (based on its gradient in the space between cells). To implement this hypothesis, we performed simulations where there was a homogeneously expressed gene that activated planar contraction in all cells and a second gene expressed only in the territory of interest whose gradient determined the direction of contraction. Therefore, the cells that were in the boundary of this territory (i.e. had at least one neighbor with different gene expression) elongated towards their neighbors with different gene expression. This created an elastic belt of elongated cells around the territory of interest. The strength of this mechanism could be regulated by changing the parameter of cell planar contraction activation (Table S5).

The non-division boundary hypothesis was implemented by regulating the division rate of the cells in the territory boundary (i.e. cells in contact with cells from the other territory). In the extreme case, cells in the border were precluded from dividing at all.

The oriented division hypothesis was implemented in two alternative manners. In the first one, we forced all cells to have their division planes in the same direction. In the second one, we forced the division planes of the cells in the territory boundary to be parallel to that boundary.

The majority rule hypothesis was implemented by simulating a gene network with a membrane-tethered ligand that can bind to a homogeneously expressed receptor. When the ligand binds to its receptor, the expression of a transcription factor is inhibited. This transcription factor is expressed in cells outside the territory of interest and inhibits the expression of the ligand. Therefore, when cells expressing the receptor have some neighbors that express the ligand, they will activate the receptor. The activated receptor directly reduces the expression of the transcription factor, and the transcription factor does the same with the ligand. The expression of the receptor is constant and the transcription factor has a slight self-activation, meaning that, in absence of activated receptor, its expression increases. This mechanism has many parameters: one parameter for each gene interaction just described, and the binding affinity of the ligand and its receptor. However, all of them affect the efficiency of the majority rule signaling in a similar way, as it is only the efficiency of this signaling that matters. Because of that, we varied only the concentration of the receptor as the parameter of this buffering mechanism. This concentration was the same in all cells.

The constant signaling hypothesis was implemented by performing simulations with a signaling center in the animal pole of the blastula. In the animal pole there is a constitutively expressed transcription factor that activates the expression and secretion of an extracellular diffusible signal 2. This signal activates the expression of two transcription factors: gene 3 and gene 4. Gene 4 only becomes activated in cells that receive high concentrations of the signal, whereas gene 3 also becomes activated at lower concentrations. Gene 3 inhibits gene 4 and, as a result, gene 4 is expressed only at a distance from the source of the signal, forming a ring-like territory around the animal pole. The simulations with constant cell signaling involved a period of growth by cell division, together with cell signaling and a final period of morphogenesis (in which gene 4 promotes cell contraction in the cells where it is expressed) with simultaneous cell signaling but without cell division. The same initial condition was simulated without constant signaling as a control. In this case, we simulated a first period of cell signaling, obtaining the induced territory, and we fixed the gene expression and stopped cell signaling. We then simulated a second phase of growth by cell division and a third phase with only morphogenesis, in which gene 4 promotes cell contraction in the cells in which it is expressed.

The above specified gene network has a number of parameters specifying which genes regulate each other and at what strength. However, the buffering capacity of this mechanism is based on signaling being constant and, as such, it does not depend on the parameters of the signaling, as long as it is constant.

Bilateral asymmetry

The morphologies arising from the simulations of the ensemble mechanisms were expected to be bilaterally symmetrical because the initial condition was bilaterally symmetric. Before calculating the bilateral asymmetry, the embryo was re-centered around its centroid (i.e. its mass center) to account for any displacement of the whole morphology during development. The morphology was then split into two halves using the bilateral plane of symmetry and the coordinates of one half were mirrored into the other. Bilateral symmetry was measured by calculating the morphological distance between the two halves. For morphological distance we used the Euclidean mean distance:
formula
where n1 and n2 are the number of nodes in each half of the embryo, dk,min(k,2) is the Euclidean distance between node k in the first half and its closest node in the other, and dj,min(j,1) is the distance between node j in the second half and its closest node in the other (Salazar-Ciudad and Marín-Riera, 2013; Hagolani et al., 2019). EMD is simply a measure of morphological dissimilarity that can be applied to morphologies with no homologous landmarks. The bilateral asymmetry has units of space in the model (a cell diameter in the initial conditions has a value of 0.14).

Territory sinuosity

The sinuosity of gene territories was calculated using two methods. When the gene expression was binary, i.e. the gene was either expressed at given level or not at all, as in Fig. 4, sinuosity was calculated as the ratio between the area and the perimeter of the territory. The area was estimated as the number of cells in the territory and the perimeter was estimated as the number of cells in contact with cells from other territories (i.e. the territory boundary). When gene expression changed gradually over space (i.e. in a gradient) (Figs 3I and 5), sinuosity was calculated as the base 10 logarithm of the sum of absolute differences in gene concentration between all cells and their neighbors.

Statistical tests

The non-parametric linear correlations were performed using the package ‘Rfit’ and the Wilcoxon non-parametric test were performed with the function ‘wilcox.test’ in R (version 3.6.3).

Implementation of single-node planar cell contraction

Planar cell contraction was already implemented in the original EmbryoMaker model (Marin-Riera et al., 2016), but it was only effective when cells were composed of multiple nodes. In this study, we implemented tissue elongation by planar cell contraction in a way that can be applied to cells represented by only one (mesenchyme) or two (epithelium) nodes. Unlike other cell properties or behaviors, planar cell contraction requires at least two genes to be activated. First, a gene determining cell polarization must be expressed. Cells will ‘sense’ the gradient of this polarizing gene, i.e. the gradient vector of the comparison of the concentration of the gene in the cell with that in its neighboring cells. This gradient determines the direction of cell contraction. Second, the level of expression of another gene will determine how much cells contract (in the direction provided by the first gene).

The nodes of planarly contracted cells are simulated as ellipses instead of spheres, the major radius of that ellipse being equal to the radius of the node (pEQD) multiplied by the sum of the expression values of gene-activating elongation multiplied by the strength of their regulatory activity. The minor radius of the ellipse is equal to the radius of the node (pEQD) divided by the sum of the expression values of gene-activating planar contraction multiplied by the strength of their regulatory activity. Therefore, planarly contracting cells can elongate while keeping their volume roughly constant. There is, however, a maximum value of elongation allowed (MMAXELO), which is a global parameter of the model.

We thank Roland Zimm, Pascal Hagolani, Kevin Martinez-Añon, Antonio Loreto and Aleksa Ratarac for comments. We also thank the CERCA Programme/Generalitat de Catalunya for institutional support and acknowledge the CSC – IT Center for Science (Finland) for computational resources.

Author contributions

Conceptualization: H.C.-F., I.S.-C.; Methodology: H.C.-F., I.S.-C.; Software: M.B.-U., I.S.-C.; Investigation: H.C.-F., T.T., M.B.-U., I.S.-C.; Resources: I.S.-C.; Writing - original draft: I.S.-C.; Writing - review & editing: H.C.-F., T.T., M.B.-U.; Visualization: H.C.-F., T.T., M.B.-U.; Supervision: I.S.-C.; Project administration: I.S.-C.; Funding acquisition: I.S.-C.

Funding

We thank the Ministerio de Ciencia e Innovación for funding (PID2021-122930NB-I00 and CNS2022-135397 to I.S.-C., and FPU19/01428 to H.C.-F.). This work is also supported by the Agencia Estatal de Investigación, through the Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence in R&D (CEX2020-001084-M).

Data availability

The model and initial conditions are available in GitHub (https://github.com/HugoCanoFernandez/EmbryoMaker).

Addison
,
M.
,
Xu
,
Q.
,
Cayuso
,
J.
and
Wilkinson
,
D. G.
(
2018
).
Cell identity switching regulated by retinoic acid signaling maintains homogeneous segments in the hindbrain
.
Dev. Cell
45
,
606
-
620.e3
.
Aliee
,
M.
,
Röper
,
J.-C.
,
Landsberg
,
K. P.
,
Pentzold
,
C.
,
Widmann
,
T. J.
,
Jülicher
,
F.
and
Dahmann
,
C.
(
2012
).
Physical mechanisms shaping the Drosophila dorsoventral compartment boundary
.
Curr. Biol.
22
,
967
-
976
.
Angelini
,
T. E.
,
Hannezo
,
E.
,
Trepat
,
X.
,
Marquez
,
M.
,
Fredberg
,
J. J.
and
Weitz
,
D. A.
(
2011
).
Glass-like dynamics of collective cell migration
.
Proc. Natl Acad. Sci. USA
108
,
4714
-
4719
.
Brozovic
,
M.
,
Dantec
,
C.
,
Dardaillon
,
J.
,
Dauga
,
D.
,
Faure
,
E.
,
Gineste
,
M.
,
Louis
,
A.
,
Naville
,
M.
,
Nitta
,
K. R.
,
Piette
,
J.
et al. 
(
2018
).
ANISEED 2017: extending the integrated ascidian database to the exploration and evolutionary comparison of genome-scale datasets
.
Nucleic Acids Res.
46
,
D718
-
D725
.
Calzolari
,
S.
,
Terriente
,
J.
and
Pujades
,
C.
(
2014
).
Cell segregation in the vertebrate hindbrain relies on actomyosin cables located at the interhombomeric boundaries
.
EMBO J.
33
,
686
-
701
.
Canela-Xandri
,
O.
,
Sagués
,
F.
,
Casademunt
,
J.
and
Buceta
,
J.
(
2011
).
Dynamics and mechanical stability of the developing dorsoventral organizer of the wing imaginal disc
.
PLoS Comput. Biol.
7
,
e1002153
.
Cayuso
,
J.
,
Xu
,
Q.
,
Addison
,
M.
and
Wilkinson
,
D. G.
(
2019
).
Actomyosin regulation by Eph receptor signaling couples boundary cell formation to border sharpness
.
eLife
8
,
e49696
.
Dahmann
,
C.
,
Oates
,
A. C.
and
Brand
,
M.
(
2011
).
Boundary formation and maintenance in tissue development
.
Nat. Rev. Genet.
12
,
43
-
55
.
Davidson
,
E. H.
(
1991
).
Spatial mechanisms of gene regulation in metazoan embryos
.
Development
113
,
1
-
26
.
Davidson
,
E. H.
,
Peterson
,
K. J.
and
Cameron
,
R. A.
(
1995
).
Origin of bilaterian body plans: evolution of developmental regulatory mechanisms
.
Science
270
,
1319
-
1325
.
Devany
,
J.
,
Sussman
,
D. M.
,
Yamamoto
,
T.
,
Manning
,
M. L.
and
Gardel
,
M. L.
(
2021
).
Cell cycle–dependent active stress drives epithelia remodeling
.
Proc. Natl Acad. Sci. USA
118
,
e1917853118
.
Foty
,
R. A.
and
Steinberg
,
M. S.
(
2005
).
The differential adhesion hypothesis: a direct evaluation
.
Dev. Biol.
278
,
255
-
263
.
García-Bellido
,
A.
,
Ripoll
,
P.
and
Morata
,
G.
(
1973
).
Developmental compartmentalisation of the wing disk of Drosophila
.
Nat. New Biol.
245
,
251
-
253
.
Gilbert
,
S. F.
and
Barresi
,
M. J. F.
(
2019
).
Developmental Biology
.
Sunderland
:
Sinauer associates, Inc
.
Gilbert
,
S. F.
and
Raunio
,
A. M.
(
1997
).
Embryology: Constructing the Organism
.
Sunderland
:
Sinauer associates, Inc
.
Graham
,
J. H.
(
2021
).
Nature, nurture, and noise: Developmental instability, fluctuating asymmetry, and the causes of phenotypic variation
.
Symmetry
13
,
1204
.
Hagolani
,
P. F.
,
Zimm
,
R.
,
Marin-Riera
,
M.
and
Salazar-Ciudad
,
I.
(
2019
).
Cell signaling stabilizes morphogenesis against noise
.
Development
146
,
dev179309
.
Inoue
,
T.
,
Tanaka
,
T.
,
Takeichi
,
M.
,
Chisaka
,
O.
,
Nakamura
,
S.
and
Osumi
,
N.
(
2001
).
Role of cadherins in maintaining the compartment boundary between the cortex and striatum during development
.
Development
128
,
561
-
569
.
Jernvall
,
J.
,
Åberg
,
T.
,
Kettunen
,
P.
,
Keränen
,
S.
and
Thesleff
,
I.
(
1998
).
The life history of an embryonic signaling center: BMP-4 induces p21 and is associated with apoptosis in the mouse tooth enamel knot
.
Development
125
,
161
-
169
.
Klingenberg
,
C. P.
(
2003
).
A developmental perspective on developmental instability: theory
. In
Developmental Instability: Causes and Consequences
(ed.
M.
Polak
), pp.
14
-
34
.
New York
,
USA
:
Oxford University Press
.
Koshida
,
S.
,
Kishimoto
,
Y.
,
Ustumi
,
H.
,
Shimizu
,
T.
,
Furutani-Seiki
,
M.
,
Kondoh
,
H.
and
Takada
,
S.
(
2005
).
Integrin5-dependent fibronectin accumulation for maintenance of somite boundaries in zebrafish embryos
.
Dev. Cell
8
,
587
-
598
.
Landsberg
,
K. P.
,
Farhadifar
,
R.
,
Ranft
,
J.
,
Umetsu
,
D.
,
Widmann
,
T. J.
,
Bittig
,
T.
,
Said
,
A.
,
Jülicher
,
F.
and
Dahmann
,
C.
(
2009
).
Increased cell bond tension governs cell sorting at the Drosophila anteroposterior compartment boundary
.
Curr. Biol.
19
,
1950
-
1955
.
Lecuit
,
T.
and
Lenne
,
P.-F.
(
2007
).
Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis
.
Nat. Rev. Mol. Cell Biol.
8
,
633
-
644
.
Ma
,
X.
,
Zhao
,
Z.
,
Xiao
,
L.
,
Xu
,
W.
,
Kou
,
Y.
,
Zhang
,
Y.
,
Wu
,
G.
,
Wang
,
Y.
and
Du
,
Z.
(
2021
).
A 4D single-cell protein atlas of transcription factors delineates spatiotemporal patterning during embryogenesis
.
Nat. Methods
18
,
893
-
902
.
Marin-Riera
,
M.
,
Brun-Usan
,
M.
,
Zimm
,
R.
,
Välikangas
,
T.
and
Salazar-Ciudad
,
I.
(
2016
).
Computational modeling of development by epithelia, mesenchyme and their interactions: a unified model
.
Bioinformatics
32
,
219
-
225
.
Monier
,
B.
,
Pélissier-Monier
,
A.
,
Brand
,
A. H.
and
Sanson
,
B.
(
2010
).
An actomyosin-based barrier inhibits cell mixing at compartmental boundaries in Drosophila embryos
.
Nat. Cell Biol.
12
,
60
-
65
.
Nielsen
,
C.
(
2004
).
Trochophora larvae: cell-lineages, ciliary bands, and body regions. 1. Annelida and Mollusca
.
J. Exp. Zool. B Mol. Dev. Evol.
302
,
35
-
68
.
Nielsen
,
C.
(
2005
).
Trochophora larvae: Cell-lineages, ciliary bands and body regions. 2. Other groups and general discussion
.
J. Exp. Zool. B Mol. Dev. Evol.
304B
,
401
-
447
.
Özpolat
,
B. D.
,
Randel
,
N.
,
Williams
,
E. A.
,
Bezares-Calderón
,
L. A.
,
Andreatta
,
G.
,
Balavoine
,
G.
,
Bertucci
,
P. Y.
,
Ferrier
,
D. E. K.
,
Gambi
,
M. C.
,
Gazave
,
E.
et al. 
(
2021
).
The Nereid on the rise: Platynereis as a model system
.
EvoDevo
12
,
10
.
Qiu
,
Y.
,
Fung
,
L.
,
Schilling
,
T. F.
and
Nie
,
Q.
(
2021
).
Multiple morphogens and rapid elongation promote segmental patterning during development
.
PLoS Comput. Biol.
17
,
e1009077
.
Rodriguez
,
I.
and
Basler
,
K.
(
1997
).
Control of compartmental affinity boundaries by hedgehog
.
Nature
389
,
614
-
618
.
Ruppert
,
E. E.
,
Barnes
,
R. D.
and
Fox
,
R. S.
(
2004
).
Invertebrate Zoology: A Functional Evolutionary Approach
.
Florence
:
Cengage Learning
.
Salazar-Ciudad
,
I.
(
2010
).
Morphological evolution and embryonic developmental diversity in metazoa
.
Development
137
,
531
.
Salazar-Ciudad
,
I.
and
Marín-Riera
,
M.
(
2013
).
Adaptive dynamics under development-based genotype–phenotype maps
.
Nature
497
,
361
-
364
.
Salazar-Ciudad
,
I.
,
Jernvall
,
J.
and
Newman
,
S. A.
(
2003
).
Mechanisms of pattern formation in development and evolution
.
Development
130
,
2027
-
2037
.
Salvador-Martínez
,
I.
and
Salazar-Ciudad
,
I.
(
2017
).
How complexity increases in development: An analysis of the spatial-temporal dynamics of gene expression in Ciona intestinalis
.
Mech. Dev.
144
,
113
-
124
.
Savriama
,
Y.
,
Vitulo
,
M.
,
Gerber
,
S.
,
Debat
,
V.
and
Fusco
,
G.
(
2016
).
Modularity and developmental stability in segmented animals: Variation in translational asymmetry in geophilomorph centipedes
.
Dev. Genes Evol.
226
,
187
-
196
.
Sepúlveda
,
N.
,
Petitjean
,
L.
,
Cochet
,
O.
,
Grasland-Mongrain
,
E.
,
Silberzan
,
P.
and
Hakim
,
V.
(
2013
).
Collective cell motion in an epithelial sheet can be quantitatively described by a stochastic interacting particle model
.
PLoS Comput. Biol.
9
,
e1002944
.
Steinberg
,
M. S.
(
1970
).
Does differential adhesion govern self-assembly processes in histogenesis? Equilibrium configurations and the emergence of a hierarchy among populations of embryonic cells
.
J. Exp. Zool.
173
,
395
-
433
.
Swaddle
,
J. P.
and
Witter
,
M. S.
(
1997
).
On the ontogeny of developmental stability in a stabilized trait
.
Proc. R. Soc. Lond. Ser. B Biol. Sci.
264
,
329
-
334
.
Umetsu
,
D.
and
Dahmann
,
C.
(
2010
).
Compartment boundaries: sorting cells with tension
.
Fly
4
,
241
-
245
.
Umetsu
,
D.
,
Aigouy
,
B.
,
Aliee
,
M.
,
Sui
,
L.
,
Eaton
,
S.
,
Jülicher
,
F.
and
Dahmann
,
C.
(
2014
).
Local increases in mechanical tension shape compartment boundaries by biasing cell intercalations
.
Curr. Biol.
24
,
1798
-
1805
.
Van Valen
,
L.
(
1962
).
A study of fluctuating asymmetry
.
Evolution
16
,
125
-
142
.
Voet
,
D.
and
Voet
,
J. G.
(
2010
).
Biochemistry
.
New Jersey
,
USA
:
John Wiley & Sons
.
Waddington
,
C. H.
(
1942
).
Canalization of development and the inheritance of acquired characters
.
Nature
150
,
563
-
565
.
Wilkinson
,
D. G.
(
2018
).
Establishing sharp and homogeneous segments in the hindbrain
.
F1000Research
7
,
1268
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information