Cell signaling stabilizes morphogenesis against noise

There is no consensus definition of complexity, yet it is evident that organisms are complex and explaining such complexity is one of the most fundamental questions of biology. In the case of the distribution of cell types in space (i.e. morphology), complexity is generated within each generation through a process of development starting from a simple initial condition, e.g. a zygote. Morphological complexity can also change between generations as part of evolution, but has not increased in the evolution of all lineages (Bonner, 2004; Williams, 1996; McCoy, 1977; Hinegardner and Engelberg, 1983; Gould, 2002; Arendt, 2008; Canning and Okamura, 2004; Arthur, 2010) and it is unclear whether there is a general trend of increasing complexity in evolution (Fisher, 1986; Ruse, 1996; Gould, 2002; McShea, 1996; but see Fleming and Mcshea, 2013). Yet, the mechanisms by which such complexity has increased in some lineages remain to be determined. How complexity increases during evolution is necessarily related to development: any evolutionary change in morphology is first a change in the developmental processes that produce such morphology.

The process of development can be described as a sequence of transformations of specific distributions of cell types in space – what we call a developmental pattern – into other, usually more complex, developmental patterns (Salazar-Ciudad et al., 2003). The first such developmental pattern would be, for example, the zygote, and the last would be the adult phenotype.

It has been argued that, despite the remarkable complexity of organisms, their development is achieved through a limited number of cell behaviors and types of cell interactions (Salazar-Ciudad et al., 2003; Davies, 2013; Newman and Bhat, 2009). These cell behaviors would be cell division, cell adhesion, cell death, cell growth, cell contraction, extracellular signal and matrix secretion, extracellular signal reception and cell differentiation (Salazar-Ciudad et al., 2003; Davies, 2013; Newman and Bhat, 2009). One may consider, in addition, cell migration and cell shape changes resulting from specific patterns of cell contraction and adhesion. Similarly, it has been argued that the diversity of animal morphology can be understood, to some extent, from a small set of ‘generic’ morphogenetic transformations that are inherent to the physical properties of animal cells and tissues (Newman and Comper, 1990; Newman and Müller, 2000).

In addition to cell behaviors, development involves interactions between cells. The two main types of cell interactions in development are cell signaling, and cell mechanical interactions arising from forces generated by cell behaviors (e.g. cell contraction). Cell signaling typically occurs through the secretion of extracellular diffusible molecules by some cells and the reception of those by other cells but it can also occur through membrane-bound signals or in other ways (Gilbert and Barresi, 2016). Both types of interactions can lead to gene expression changes. Gene expression changes lead, in turn, to changes in the behaviors of cells (Gilbert and Barresi, 2016). Which cells express which genes is also affected by cell behaviors because these affect the distribution of cells in space and this, in turn, affects the distribution of extracellular signals and forces in space (Salazar-Ciudad et al., 2003).

In this article, as previously (Salazar-Ciudad et al., 2003), we define a ‘developmental mechanism’ as a network of gene interactions, cell interactions and cell behaviors required for the transformation of one developmental pattern into another (see Fig. 1). Here, the network of gene interactions, or gene network, is defined as the set of interactions that can occur between a set of genes as a result of their genetically encoded structure. This is a specification of which gene products will interact with each other if they coincide in time and space. Whether two gene products coincide in time and space is determined by the dynamics of the developmental mechanisms. The dynamics depend on the developmental mechanism but also on things like the environment or the initial developmental pattern.

The question we want to approach in this study is: how are these interactions and cell behaviors coordinated such that they produce complex and robust morphologies, i.e. are there some logical requirements that developmental mechanisms should fulfill in order to lead to complex robust morphologies? Are there, for example, some requirements at the level of gene network topology or at the level of cell behaviors and their coordination during development?

If, as suggested above, pattern transformations in development involve a limited set of cell behaviors and cell interactions, then any mathematical model implementing those and intracellular gene networks should be able to reproduce, to a large extent, the range of pattern transformations possible in animal development. In this work we use one such model, EmbryoMaker (Marin-Riera et al., 2015), to simulate many possible developmental mechanisms and try to discover what, if anything, the mechanisms leading to robust complex morphologies have in common. To that end, we randomly wired genes and cell behaviors into developmental mechanisms and simulated, with EmbryoMaker, the morphologies they produce from a simple initial developmental pattern. This initial condition consists of a small flat epithelium with an underlying layer of mesenchymal cells (see Fig. 1A). We did not consider the mesenchyme when measuring complexity and developmental instability because in our simulations the morphology of the mesenchyme tended to mimic that of the epithelium or was very noisy.

Real developmental mechanisms are not random but the result of millions of years of evolution. However, the study of random developmental mechanisms allows us to identify general requirements without being conditioned by our current understanding of development. Despite its statistical nature, our approach should also be informative about evolution as the developmental mechanisms found in current animals may still need to fulfill these general logical requirements.

EmbryoMaker is a general mathematical model of animal development in the sense that it can simulate all the basic behaviors of animal cells, extracellular signal-mediated interactions and mechanical interactions between cells as well as any arbitrary intracellular gene network (see Figs 1 and 2). EmbryoMaker represents cells as a set of parts (herein called nodes) with specific mechanical properties. Mesenchymal cells are represented by single spherical nodes and epithelial cells by a cylinder consisting of two nodes (one basal and one apical bound by an elastic link). Nodes touching each other experience adhesion forces, but, if they get closer than a given distance, they experience a repulsive force (Fig. S1). For mesenchymal cell biomechanics, our model is based on the subcellular element model (Newman, 2005). Epithelial cells exert additional forces into neighboring epithelial nodes that reflect their specific mechanical properties and their organization in epithelial sheets (Fig. S1 and Materials and Methods).

Within each node there are gene products that can affect the expression of other gene products and the mechanical properties and cell behaviors of their containing nodes (e.g. their size, distance at which repulsion forces apply, adhesion, mitosis rates, elasticity, etc.). Some of these molecules can diffuse in extracellular space and affect cells other than the cell in which they are produced, thus leading to extracellular signaling. Nodes' mechanical properties determine how they respond to forces (for example, node size, node elasticity). At the mathematical level, the concentration of each molecule in a node, a node's mechanical properties and a node's 3D position are continuous variables that are calculated by differential equations that take into account some of these same variables in neighboring nodes [see Materials and Methods and Marin-Riera et al. (2015) for details].

EmbryoMaker is not a model of the development of any specific system. Models for specific systems, e.g. a model of limb development, are built in EmbryoMaker by specifying a concrete developmental mechanism: which gene products regulate other gene products, mechanical properties and cell behaviors in a specific system. EmbryoMaker has been applied, so far, to spiral cleavage (Brun-Usan et al., 2017) and tooth morphogenesis (Marin-Riera et al., 2018). As a result of an initial developmental pattern and a developmental mechanism, EmbryoMaker simulates development and outputs how the spatial distribution of cells and gene expression in the initial pattern change over time until some final developmental pattern (a morphology) is reached (Fig. 1B). On these virtual morphologies we can measure complexity and robustness.

There is no consensus definition of complexity and any study using one is likely to be controversial. Given this, we use two different quantitative measures of morphological complexity, and stress that our results do not necessarily apply to complexity in general but to complexity as defined by these measures. There are other possible measures of complexity (Saunders et al., 1999; McShea, 1996; Fleming and Mcshea, 2013) but they are not as easy to quantify in an objective way as we attempt here or are not easy to apply directly to 3D morphologies. Roughly, our two measures reflect the likelihood of randomly guessing the position of a cell in 3D space knowing the position of its neighbors at different distances but knowing nothing about the developmental mechanism that produced such morphology. The first measure of complexity we use is angle variance, or AV, which is the variance of the angles between the polarity axis of epithelial cells at different distance intervals (see Materials and Methods and Fig. S2 for a full description and Fig. S3 for examples). The second measure we use is the so-called ‘orientation patch count’ or OPC (Evans et al., 2007), which is the number of regions with the same slope orientation (see Materials and Methods and Fig. S4 for a full description and S3 for examples).

Robustness is also a concept that is defined and understood in many ways in the literature (Arjan et al., 2003; Salazar-Ciudad, 2007). Here, we restrict ourselves to robustness understood as the suppression of developmental instability (Waddington, 1942; Shapiro, 1971; Klingenberg, 2002). Developmental instability reflects how different are the morphologies of genotypically identical individuals that develop in exactly the same environment. This is morphological variation arising from noise in the developmental process itself. To measure developmental instability each developmental mechanism is simulated several times from the same initial developmental pattern but with noise. Noise is implemented as small random displacements of the positions of all cells in each iteration of the model. The morphological distance between the resulting embryos is then a measure of developmental instability. We use two different complementary measures of morphological distance: AV and OPC (see Fig. S5 for examples and Materials and Methods and supplementary Materials and Methods, section 3.1, for details).

We ran 100,000 random developmental mechanisms, in what we call the ‘broad ensemble’, and found very few morphologies differing from the initial conditions in a non-trivial way. We found also that nearly all random developmental mechanisms were unable to change gene expression patterns over space. Most genes were expressed only in the most central cell, in its immediate neighbors or not at all. As a result, most cells did not activate any cell behaviors nor did they change any mechanical property over space and, thus, there were no changes in cell positions or morphogenesis.

To circumvent this problem, we made a simpler ensemble, the ‘signaling-only ensemble’, in which cells were not allowed to move, grow or divide. As cells do not move in this ensemble, no morphogenesis occurs. In this ensemble, we identified which developmental mechanisms can produce temporally stable changes in gene expression over space, as in a previous publication (Salazar-Ciudad et al., 2000), so that genes are not only expressed in the most central cell or its immediate surroundings but more widely over the embryo. We then used the networks identified this way to construct another ensemble, the ‘signaling ensemble’, by making some of the genes in the networks regulate some randomly chosen node properties or cell behaviors at randomly chosen intensities (see Materials and Methods and supplementary Materials and Methods, section 2). In this and subsequent ensembles, cells are allowed to move due to the regulation of cell behaviors and mechanical properties. As a result, morphogenesis can occur and a morphology different from that of the initial pattern (a flat epithelium) can arise. We constructed 20,000 developmental mechanisms in such ensembles. In the signaling ensemble, an initial developmental pattern with a gene expressed in a gradient was used (see Materials and Methods and Fig. 3). A planar cell polarity (PCP) ensemble was also built. This ensemble is exactly like the signaling ensemble but includes a set of cells that secrete, constitutively, an extracellular signal that diffuses over space. In this PCP ensemble cells polarize the direction of their cell division and movement in the direction of the gradient of this signal. This ensemble was built to explore the effect of PCP on morphological complexity and developmental instability.

Additionally, the developmental mechanisms in the signaling ensemble were used to build three other ensembles. These ensembles were built to explore the morphological complexity and developmental instability that should be possible without cell-cell signaling. These were then compared with the signaling ensemble to understand better the influence of signaling on morphological complexity and developmental instability. The ‘gradient autonomous ensemble’ uses the same developmental mechanisms as in the signaling ensemble but without extracellular signaling and keeping the initial gene gradient. With this ensemble we explored the morphogenesis that is possible from the gradient regulation of cell behaviors. In the ‘autonomous ensemble’, there is no extracellular signaling and one gene is homogeneously expressed through all cells whereas the rest of genes are not expressed at all. With this ensemble we explored the morphogenesis that is possible from the uniform regulation of cell behaviors. In the ‘autonomous biomechanics-only ensemble’ there are only three non-interacting genes and these are expressed everywhere in the initial developmental pattern. With this ensemble we explored the morphogenesis that is possible from the uniform regulation of cell behaviors with no changes in gene expression over time. See Fig. 3 and Tables S1 and S2 for an overview of the different ensembles.

In all ensembles we found complex morphologies, although at a low frequency (see Fig. 4, Fig. S3 for some example morphologies). In addition, we found that the developmental mechanisms producing complex morphologies tend to be more developmentally unstable than the developmental mechanisms producing simple morphologies (see Fig. 4I-L, Fig. S6).

Most remarkably, we found that the ensembles with signaling (the signaling and PCP ensembles) produced morphologies that were, on average, significantly less complex than the morphologies produced in the two autonomous ensembles (see Fig. 4). No significant differences were found between the signaling and PCP ensembles. The within-ensemble disparity was slightly smaller in the signaling ensemble than in the autonomous ensembles (Fig. S7). Disparity was measured as the sum of the morphological distances between each morphology in an ensemble and all other morphologies in the same ensemble divided by the number of distances measured.

The ensembles differed in developmental instability (see Fig. 4, Fig. S6). For morphologies of the same complexity, more developmentally stable morphologies were found in the signaling ensembles and in the gradient autonomous ensemble than in the other autonomous ensembles (Fig. S8). In other words, although extracellular signaling or gradients did not seem to be required for complex morphologies, they were required for complex morphologies to be developmentally stable.

Most complex morphologies in the autonomous and autonomous biomechanics-only ensembles consisted of highly folded epithelia, like crumpled paper balls, with lumens or mesenchymal cells inside. Often the position of epithelial folds was different in each run of the same developmental mechanism. In the signaling ensembles and in the gradient ensemble, many morphologies were also composed, totally or partially, of randomly folded epithelia but there were also complex morphologies made of folds that consistently appeared in the same location and, thus, were complex yet developmentally stable (see Fig. S9 for examples of both kinds of morphologies).

We found that the development of complex morphologies requires only that, over developmental time, a large proportion of the cells in an embryo change the regulation of cell behaviors or mechanical properties (Fig. S10). Cell contraction and cell division were the cell behaviors most often associated with the development of complex morphologies (Fig. S11). This is especially the case, but not exclusively, when cell contraction is asymmetric between the apical and basal side of each epithelial cell, as found in the formation of invagination and tubes in many animals (Martin and Goldstein, 2014). Cell division could lead to some buckling and wrinkling of the epithelium and then to some complexity, as in fact is also observed in the development of many animal organs (Bunn et al., 2011; Striedter et al., 2015). This was especially the case if the epithelium had regions with different values in nodes' mechanical properties or different rates of cell division. The regulation of other mechanical properties and cell behaviors also had an effect on complexity but to a lesser extent and only when contraction was also present (Figs S12 and S13 and Tables S4 and S5 for their statistics; see Table S6 for a summary of the mechanical node properties).

The emergence of complexity from homogeneous initial developmental patterns, especially the crumpled paper morphologies, was studied by following their development over time. We observed that the contraction of all cells in an epithelium leads to a homogeneous increase in its curvature. This increase, for the geometrical reasons described in Fig. 5, can resolve into the formation of several different invaginations. We found that the position of each such invagination is largely affected by noise. The size and shape of each invagination, instead, is affected by the degree of cell contraction: the larger the contraction, the larger the curvature and more and smaller invaginations form. As invaginations form, they affect each other's shape through mechanical interference and partial fusion. As a result of such interactions and the effect of noise, embryos' epithelia do not fold into a regular array of invaginations but, instead, into an intricate sea of partially fused invaginations (as in Fig. S9C). These embryos are complex precisely because it is hard to predict the position of one cell based on the position of its neighbors.

In order to understand better how extracellular signaling enhances robustness against noise, we first paid attention to the effect of noise on development when there is cell contraction because this is the cell behavior most often associated with complex morphologies in our simulations. Qualitatively, we found that if cells contract at the same time and with the same intensity over large regions of the embryo the resulting morphologies tend to be complex but unstable. If, by contrast, cell contraction occurs in different ways (at different moments or at different rates) over different regions of the embryo, development tends to lead to more complex but also more stable morphologies.

On purely geometric grounds it can be seen (see Fig. 5) that only a very specific number of cells can fit into an invagination. The larger the cell contraction, i.e. the larger the apical side is compared with the basal side or vice versa, the smaller the number of cells within the invagination is and the higher the curvature of the resulting invagination. This implies that, in our model, large fields of contracting cells will necessarily split into different invaginations and we found in such cases that the positioning of such multiple invaginations is quite sensitive to noise.

By following the development of many complex multi-invaginated morphologies in the ensembles, we found that the positioning of invaginations was highly dependent on noise and that the instability to noise was related to the non-linear nature of the invagination process. As an invagination starts to form in our in silico embryos, its cells start to rotate on its longest axis to align it with that of its neighbors that are also rotating. The rotation is, thus, in the direction in which neighbor epithelial cells have already contracted and rotated the most. In turn, the direction in which a cell rotates also affects the direction in which its neighbor cells rotate, further strengthening the rotation in one direction or another. When invaginations start to form, slight noise in the timing, rate or direction in which cells contract, gets easily amplified over time and affects where each invagination will form.

To explore the quantitative support of these qualitative observations, we ran several simulations in which epithelia of different sizes contracted all their cells in the same way and at the same time (as in the cell behaviors-only ensemble). As can be seen in Figs 6 and 7, the larger the epithelium, the larger is its developmental instability. Splitting the epithelium into regions contracting in slightly different moments or at slightly different rates, however, decreased developmental instability (Fig. S14). The same occurred if the epithelium was split into equally contracting regions separated by narrow non-contracting boundaries (Fig. S14). In other words, compartmentalizing the embryo in different regions largely reduces developmental instability without precluding the development of complex morphologies.

If developmental instability is related to the size of the epithelial regions contracting in the same way, then developmental instability and such size should correlate in the ensembles we simulated. Fig. 8 shows that this was indeed the case: the developmental stability of an embryo correlates with the size (in number of cells) of the largest regions in which cells are contracting in the same way (i.e. changing their apical or basal side at the same rate). These regions are larger in the autonomous ensemble and cell behaviors-only ensemble because, in such ensembles, all cells behave in exactly the same way. All the embryos in Fig. 8 (and most embryos in the ensembles) have the same number of cells so the relationship we found was not due to larger embryos being less stable.

The stabilizing effect of extracellular signaling was not found to be related to cell division. In contrast with contraction, compartmentalizing the embryo into regions where cell division would occur at different moments or in slightly different ways, as done for contraction in Fig. S14, had no clear effect on developmental instability (Fig. S15).

One main question of this article is whether there are some features of developmental mechanisms, e.g. at the level of gene network topology or cell behaviors, that are required for the development of complex robust morphologies. For the case of gene network topology, we found the answer to be ‘no’. The fact that complex morphologies were less common in the signaling ensembles than in the ensembles without signaling suggests that the development of morphological complexity as such does not require extracellular signaling, at least for the range of complexity observed in our study. Instead, we found that to produce complex morphologies the main requirement is that cell behaviors, especially cell division and contraction, should be activated over large regions of the embryo. One of our main results is also that gene networks and extracellular signaling facilitate that complex morphology to be robust to noise. In other words, even if the mere activation of cell behaviors can lead to complexity, for such complexity to be robust complex gene networks and extracellular signaling are also required.

In our simulations, complex morphology is often associated with cell division and contraction. That the other cell behaviors have only a modifying effect is not surprising. In epithelia, cell division and contraction are the cell behaviors that can generate forces (Bard, 1990) and, thus, lead to cell movement leading to morphogenesis. The secretion of extracellular matrix can also generate forces, but we rarely found such cell behavior associated with the development of complex morphologies in the ensembles. Other behaviors, e.g. cell adhesion, either do not generate forces, i.e. they only resist forces incoming from the environment or other cells or have a morphogenetic effect mostly when considering the mesenchyme. We only consider the morphology of the epithelium when measuring complexity. In our simulations, the morphology of the mesenchyme tended to mimic that of the epithelium or was very noisy. Simulating and measuring the complexity of mesenchyme in more detail, however, would be more computationally challenging than the present study. Thus, we cannot currently be sure that our results would apply when considering the mesenchyme too. However, as the morphology of epithelia is a large determinant of the body's morphology, we expect our results to provide insights for complex morphology and robustness. Cell polarity, such as in PCP, does not seem to have an effect on complexity, at least for the simple gradient of polarity we use; it simply leads to embryos that are more elongated in one direction than in others.

Our results also indicate that for complex morphologies to be developmentally stable, cell behaviors should not be activated in the same way over large regions of the embryo. This heterogeneity can be achieved by gradients of gene expression, as in the gradient autonomous ensemble, or by extracellular signaling compartmentalizing the embryo into different small compartments of gene expression, as in the signaling ensembles. In either case, cell behaviors can then be activated differently among contiguous compartments and, thus, lead to stable development. In addition, cell signaling allows for cells having specific patterns of gene expression (e.g. cell types) to be located in specific parts of a morphology.

The developmental mechanisms that can compartmentalize the embryo into regions of differential gene expression are those found in the signaling-only ensemble. These are essentially the same networks found in other studies using the ensemble approach only with extracellular signaling (Salazar-Ciudad et al., 2000). In broad terms, these can be classified into Turing-like reaction-diffusion mechanisms (Turing, 1952) or hierarchic-signaling developmental mechanisms (Salazar-Ciudad et al., 2000).

Our results on developmental stability suggest the existence of a strong constraint on the architecture of development. To ensure developmental stability, development should be structured in such a way that there are never large regions of cells regulating cell behaviors in the same way. Instead, the embryo should be undergoing constant compartmentalization so that, as it grows, the size of these compartments remains small enough and most of the embryo is subdivided into such compartments. This is consistent with the observed existence of compartments in organ development and their progressive sub-compartmentalization as they grow and develop further, most notably in wing imaginal discs (Garcia-Bellido et al., 1973) and brains (Kiecker and Lumsden, 2005). This progressive compartmentalization requires constant extracellular signaling as the embryo grows and deforms during morphogenesis. This simultaneity between signaling and cell movement has been previously suggested to be important in animal development, although on different grounds (Salazar-Ciudad et al., 2003) and especially for vertebrates (Salazar-Ciudad, 2010).

Our results are also consistent with the large, and mostly random, morphological variability of embryoids (Simunovic and Brivanlou, 2017), in which we observe how the simple activation of cell behaviors can lead to relatively complex morphologies. Embryoids lack the spatially confined patterns of gene expression observed in vivo and that in our study are required for developmental stability (Simunovic and Brivanlou, 2017). This, together with the inevitably noisier in vitro environment, may explain the large developmental instability of embryoids' morphology.

Our result that relatively complex morphologies can be attained by the spatially homogeneous activation of cell behaviors, most notably contraction, can be understood by considering that the mechanical interactions between cells are usually non-linear (Oster and Alberch, 1982; Forgacs and Newman, 2005; Taber, 2014). Thus, as in Turing-like reaction-diffusion systems, simple homogeneous developmental patterns may be unstable to small perturbations on the system variables, e.g. cell positions, and easily change to more stable but spatially non-homogeneous configurations, e.g. an invagination.

Ours is not the only report of complex morphologies developing from homogeneous fields of cells that mechanically interact without extracellular signaling. Certain aspects of organ morphology, such as gut folding (Savin et al., 2011; Thomason et al., 2012; Nerurkar et al., 2017) or brain cortical folding (Bayly et al., 2013), have been shown experimentally to develop from homogeneous fields of cells without extracellular signaling being involved in the process. In more general terms, it has been argued that the development of many morphologies can be explained as a simple consequence of the mechanical properties of cells and their extracellular matrix, even from homogeneous initial conditions (Newman and Comper, 1990).

The relevance of extracellular signal gradients for the robust specification of the different regions of gene expression in embryos, i.e. ‘patterning’, has received extensive attention (Hogeweg, 2000; Barkoulas et al., 2013). To our knowledge, however, their role in decreasing developmental instability at the level of morphogenesis has not been suggested before.

The questions we aimed to address in this study could not have been addressed before, at least directly. The closest study would be that by Hogeweg (2000) using a 2D Potts model with a Boolean gene network, extracellular signaling and cell adhesion to show that morphogenesis can evolve as a side effect of natural selection for diversity of cell types. Another study using a similar kind of model is that of Nissen et al. (2018). This is a 3D model of epithelia with planar cell polarity, cell division and adhesion. The model does not include other cell behaviors, extracellular signaling or gene networks. The initial conditions consist of balls of randomly positioned epithelial cells (e.g. cylinders as in our model). Each different initial condition gives rise, over development time, to a very distinct final morphology. Development is, thus, robust to noise but very dependent on the specific arrangement of epithelial cells in the initial condition.

There are several caveats related to our results. First, one may argue that complex animals are not just randomly folded epithelia and that, then, many of the morphologies we classify as complex do not necessarily resemble animal embryos. The way we measure complexity classifies seemingly random morphologies as more complex than morphologies that seem more animal-like (see Fig. S9). In fact, however, many existing measures of complexity, although not directly applicable to 3D morphology, acknowledge that what seems to be random may indeed be quite complex (Kolmogorov, 1998; Wolfram, 2002). For us to do otherwise, we would have to propose an unequivocal way to differentiate biological from non-biological complexity, a daunting task in itself. The only obvious difference we can detect is that, at least in most animals, there are three different axes of symmetry. In any case, both seemingly random and animal-like morphologies are classified as being rather complex according to our measures.

Second, our ensembles only consider small numbers of genes and cells and, thus, cannot lead to morphologies that look as complex as those of most actual animals. The only possible exceptions would be the morphologies of cnidaria and ctenophora, although these certainly involve developmental mechanisms that are more complex that the minimal ones we found (Abdol et al., 2017). In that respect, the inferences we make on metazoan evolution should be regarded as primarily applying to the early evolution of animals with true epithelia or eumetazoa (i.e. all animals but sponges and some other minor groups; Fidler et al., 2014).

Third, for computational reasons, we only simulated developmental mechanisms with up to ten genes and we stopped simulations at embryos of 5000 cells. We think, however, that our conclusion that large fields of cells activating the same cell behaviors would tend to be unstable and that partitioning such fields into smaller subfields increases their developmental instability should still hold. In order to shorten computational time, we also considered that each cell is made of a single cylinder or sphere. This precluded us from simulating PCP polarized cell contraction because this would require each cell to be made of many nodes, so that different parts of the cell can contract at different rates. In that sense we cannot conclude that polarized contraction does not have an effect on complexity and robustness.

The morphologies found all resemble each other in some fashion and are clearly distinct from those found in plants and fungi. This argues for the kinds of animal-tissue materials modeled by EmbryoMaker having inherent forms, a point that has been made in the literature for animal development (Newman and Comper, 1990). According to a hypothesis by Newman and Müller (Newman and Comper, 1990; Newman and Müller, 2000; Newman et al., 2006, Müller, 2006) early metazoans had relatively complex but very unstable morphologies. These authors argue that the behaviors and mechanical properties of animal cells allow for a relatively large repertoire of relatively complex morphologies. Strikingly, this is what we found in the autonomous and cell behaviors-only ensembles: complex unstable morphologies arising by the activation of cell behaviors and mechanical interactions without extracellular signaling or even without complex gene regulatory networks. In addition, these authors argue that, later on, these complex and unstable metazoans evolved stable morphologies through the recruitment of complex gene networks in development. Although these authors do not provide much detail on how this happens, our results are consistent with the view that the early function of developmental gene networks and extracellular signaling may have been in stabilizing development rather than in building complex morphology per se.

The above argument by Newman and Müller concerns early metazoan evolution. In many current eumetazoa, gene networks and extracellular signaling are pervasively important in the construction of morphology (Gilbert and Barresi, 2016). In addition, current complex eumetazoan morphologies consist of more than folded epithelia. These two facts suggest that beyond the earliest metazoan evolution the role of gene networks and extracellular signaling is not restricted to making complex morphologies stable but also extends to further increasing possible morphological complexity. This could be achieved by recombining existing developmental mechanisms in different stages and body parts. In other words, the use of gene networks and extracellular signaling allow a finer partitioning of the embryo, in each different stage, into territories and activate different developmental mechanisms in each of them. Although the basic morphologies are still those possible from the behaviors and mechanical properties of cells (such as rods, invaginations, cavities, etc.) (Newman and Comper, 1990; Newman and Müller, 2000) these are recombined through the different parts of the embryo to construct the complex and slightly modular anatomies observed in many current eumetazoa.

In conclusion, this study suggests that gene networks and signaling have a crucial role in stabilizing the development of complex morphologies against noise. Such stabilization is mediated through the compartmentalization of the embryo into small regions where cell behaviors are regulated differently. Without this compartmentalization, complex but non-robust morphologies are typically produced. We believe that the relationship between gene networks, signaling, compartmentalization and developmental instability identified in this article is general to animal development and, thus, should have wide consequences for the evolution of morphological complexity in animals.

A full description of EmbryoMaker can be found in its original publication (Marin-Riera et al., 2015). The cell behaviors considered in this article are: apoptosis, cell contraction (which can be asymmetric between the apical and basal side of epithelial cells and includes also cell expansion), cell division (and accompanying cell growth) and extracellular matrix secretion. In addition to cell behaviors, there are also a number of node mechanical properties such as their size, morphological plasticity (the plastic reduction of cell's size due to external pressure), cell adhesion affinities, resistance to compression and, for epithelial cells, resistance to epithelial bending. See supplementary Materials and Methods for a more detailed description of the properties considered in this article.

Random developmental mechanisms were built in four different ways. The set of morphologies that arise from each of these ways is termed an ensemble. Here, we describe the basics of these ensembles.

All simulations started from the same simple initial developmental pattern: a flat hexagonal sheet of 126 epithelial cells and an underlying layer of 126 mesenchymal cells (see Fig. 1A). Although EmbryoMaker allows for each cell to be made of several nodes, in this work, for simplicity, each epithelial cell was represented by a single cylindrical node and each mesenchymal cell by a spherical node. The initial values of the mechanical properties were the same in all cells in the initial developmental pattern. These initial values were such that, if unmodified over time, no changes in morphology would occur (see supplementary Materials and Methods, section 2.8 and Tables S7 and S8 for details on the parameter values).

Each developmental mechanism was built by making randomly chosen genes regulate other randomly chosen genes, mechanical properties and cell behaviors (Fig. 1B). Developmental mechanisms were randomly built but they were not random over time: once a developmental mechanism is built it does not change over time.

When building the gene network of a developmental mechanism each gene had a 50% chance of being either an extracellularly diffusible gene product (here we call these growth factors) or an intracellular gene product. The latter had a 0.25 probability of localizing at the apical side of the epithelial cells (0.25 probability for the basal side) and a 0.5 probability of localizing in both. We specified that gene 1 always directly activates a gene that can diffuse extracellularly. We built gene networks of ten genes by randomly wiring genes, each gene having a 0.2 probability of being a regulator of another gene in the network. Each regulation was, with equal chance, either positive or negative (transcriptional activation and repression) with a random regulative strength between 0 and t_max with uniform distribution. Thus, every gene had, on average, two positive and two negative connections (two efferent and two afferent). See supplementary Materials and Methods, section 2.8.1, for a description of how was t_max chosen. We considered only transcriptional regulation between gene products, although EmbryoMaker can implement regulation at other levels and molecules other than gene products.

When building a developmental mechanism, each gene is given a chance to regulate a mechanical property or a cell behavior (with each gene having a 0.5 probability of doing so). The value of such regulation was random with a uniform or logarithmic distribution that depends on the mechanical property and cell behaviors (see supplementary Materials and Methods). Each gene product had a randomly chosen degradation and diffusion rate.

In addition, all cells had a default activation of cell division and cell differentiation that was constant over time (in addition to the regulation they could receive from a developmental mechanism). This reflects the fact that cell divisions take place in basically every developing embryo. Cell differentiation causes cell behaviors to slow down over developmental time. Including this cell differentiation was motivated by the widespread slowing down of growth during embryonic development.

All simulations were numerically integrated using the order 4 Runge–Kutta method with a dynamic step size. Simulations were run for a fixed number of iterations that is roughly equivalent to three physical days (see supplementary Materials and Methods, section 1.7) unless largely aberrant morphologies were produced (e.g. consisting of broken epithelia). The number of iterations was chosen so that most embryos will finish development before this time. If largely aberrant morphologies were produced, e.g. broken epithelia (see supplementary Materials and Methods, section 2.1.7, for a full description), these were considered inviable and discarded.

The developmental mechanisms used in this ensemble are completely random, as described in the previous section. One gene is initially expressed in one cell in the center of the epithelium. This naive approach failed to find developmental mechanisms able to produce morphogenesis in most of the 100,000 random developmental mechanisms explored.

In this ensemble, cells were not allowed to move, grow or divide, but otherwise used the same strategy to make developmental mechanisms described before, although with no cell behavior being regulated. With this ensemble we identified which developmental mechanisms were able to lead to changes in gene expression over space in a temporally stable fashion (see Fig. S16), as in a previous publication (Salazar-Ciudad et al. 2000).

This ensemble was constructed using the developmental mechanisms identified in the signaling only ensemble. Then, as explained in the section 'Building random developmental mechanisms', some randomly chosen genes will regulate a cell behavior or node property, also chosen at random. The gene initially expressed, instead of being only in the central cell, is expressed in a gradient in the initial developmental pattern (Fig. 2I).

This ensemble is exactly as the signaling ensemble, but includes nine cells, at the margin of the epithelium, that constitutively secrete an extracellular signal. This signal diffuses over the embryo producing a concentration gradient (each cell's polarization vector points in the direction where the signal concentration decreases faster; see Fig. 2F). The polarization vector of each cell biases the direction of cell division and cell movement (Fig. 3E and supplementary Materials and Methods).

In this ensemble there is no extracellular signaling and one gene is homogeneously expressed through all cells in the initial developmental pattern whereas the rest are not expressed (Fig. 3C). Gene expression, thus, does not change in space but can change over time as a result of the model dynamics. Even if gene expression is homogeneous the biomechanical interactions between cells can lead to morphogenesis (a symmetry break is induced by noise and boundary conditions).

This ensemble is like the autonomous ensemble but without a regulatory gene network. Thus, gene expression does not change in space or in time but, as in the previous ensemble, morphogenesis can occur. Developmental mechanisms include only three genes: one is expressed in the apical part of cells, one in the basal side and one in both. These genes activate the cell behaviors that were regulated by the original developmental mechanism in the signaling ensemble and with the same intensity (see Fig. 3D, supplementary Materials and Methods, section 2.6).

This ensemble is just like the autonomous ensemble but only one gene is initially expressed, and it is expressed in a gradient over the epithelium (Fig. 3B).

Our measures of complexity are related to predictability, i.e. how likely it is to predict the position of an epithelial cell knowing the position of its neighbors' cells. In a flat epithelium, for example, one can easily predict the position of a cell from the position of its closest neighbors because it would have the same position in the z-axis. In a highly folded epithelium, this would be very difficult, unless the epithelium happens to fold regularly, for example following a sinusoid wave. Thus, very complex morphologies are folded irregularly. See Fig. 4 and Fig. S3 to get an intuitive idea of each complexity measure and to see the complexity of several example morphologies. We used two different measurements of complexity: angle-distance variance and orientation patch count.

This measure is based on the variation of angles between epithelial cells. The angle between two cells is calculated as the angle between the apical-basal vectors of cell i and the apical-basal vector between cells i and j (see Fig. S2A).

where D is the distance interval in which a cell has to fall in to be included in the category c, c defines the maximal and the minimal distance for each interval and is the average of all epithelial cells in the embryo. We started with c=3 to preclude noise from affecting the measurement. We continued up to c=9, to consider the macro-structure of the embryo.

where i is each of the epithelial cells, n is the total number of epithelial cells in the embryo, c is each of the category intervals and V_ic is the angle variation for cell i in the category c. Note that with this measurement a perfect sphere will have zero complexity.

Orientation patch count (OPC) is based on the number of differently oriented slope patches an epithelium has. This measure is a fully 3D version of a measure of tooth complexity that has been found to correlate with diet (Evans et al., 2007).

For this method, we first assigned each epithelial cell to one of eight categories. Each category corresponds to one octant (one of the eight divisions of a Euclidean 3D coordinate system defined by the signs of the coordinates, see Fig. S4B). To determine in which octant the basal node is, we simply checked the sign of each of the dimensions of the vector from the apical to the basal node of each epithelial cell.

Each cell was then further classified as belonging to a specific patch. A patch is a set of cells belonging to the same orientation category (of the eight possible ones) and globally connected to each other. This means that one can go from any cell in a patch to any other cell in the patch through a sequence of contiguous cells belonging to same orientation category (see Fig. S4B). By contiguous cells we mean cells that are in contact. Only patches with more than three cells were considered. Finally, we simply counted the number of patches in a morphology, which gave us the OPC value.

We define developmental instability as the morphological distance between the morphologies resulting from simulating the same developmental mechanisms with or without noise affecting how cells move. As with complexity measures, we only took into account epithelial cells. We used two different complementary measures of morphological distance: the Euclidean morphological distance and the orientation morphological distance (see supplementary Materials and Methods, section 3.0).

where n₁ and n₂ are the number of nodes in morphology 1 and 2, respectively, d_k,min(k,2) is the distance between node k in morphology 1 and its closest node in morphology 2, d_j,min(j,1) is the distance between node j in morphology 2 and its closest node in morphology 1. Note that one node in morphology 1 being the closest node to another node in morphology 2 does not imply that this latter node is the closest to the former (i.e. these minimal distance relationships are not symmetric).

We thank Miguel Brun-Usan, Lisandro Milocco, Renske Vroomans and Jukka Jernvall for useful comments, and Sophia Hagolani-Albov for correcting grammar and spelling in the text. The CSC (Finnish IT Center for Science) provided computer time. We thank the Konrad Lorenz Institute for its support.