On the origins of developmental robustness: modeling buffering mechanisms against cell-level noise Free

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.

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.

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

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

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.

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.

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.

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.

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.

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.

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.

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

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.

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.

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.

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

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.

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.

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

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.

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

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⁻⁸ (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, e_ij, for a given cell behavior or node property i in a developmental mechanism j, was chosen to be e_ij=0.75*pⁱ_max g_j. Where pⁱ_max 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). g_j 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 g_j and pⁱ_max are not parameters, i.e. their values are a result of model dynamics. The equation for e_ij 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 e_ij^*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.

In each iteration, a proportion of the nodes was chosen, which was designated the model parameter M_NOI (Marin-Riera et al., 2016), and noise was applied to them. M_NOI 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.

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 p^ADH; 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 b_ii 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, b_ii (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.

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.

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

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 (p^EQD) 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 (p^EQD) 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 (M_MAXELO), 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.

The peer review history is available online at https://journals.biologists.com/dev/lookup/doi/10.1242/dev.201911.reviewer-comments.pdf.