Theory bridging cell polarities with development of robust complex morphologies

Despite continual renewal and damages, a multicellular organism is able to maintain its complex morphology. How is this stability compatible with the complexity and diversity of living forms? Looking for answers at protein level may be limiting as diverging protein sequences can result in similar morphologies. Inspired by the progressive role of apical­basal and planar cell polarity in development, we propose that stability, complexity, and diversity are emergent properties in populations of proliferating polarized cells. We support our hypothesis by a theoretical approach, developed to effectively capture both types of polar cell adhesions. When applied to specific cases of development – gastrulation and the origins of folds and tubes – our theory makes experimentally testable predictions pointing to the strength of polar adhesion, initial and boundary orientation of cell polarities, and the rate of cell proliferation to be major determinants of morphological diversity and stability.


INTRODUCTION
Multicellular organisms are amazing in their ability to maintain complex morphology in face of continuous cell renewal and damages. Adult salamander can regenerate entire limbs (Eguchi et al. 2011) , and, during development, some regions can maintain patterning when moved to different parts of an embryo or if the size is varied (Lyons, Kaltenbach, and McClay 2011) . Given the vast complexity and diversity of living shapes, how can we reconcile the robustness to perturbations with flexibility to diversify? While undoubtedly the end result is encoded in the DNA and protein networks, looking for an answer at this level is challenging. Examples of phenotypic plasticity (Libby and Rainey 2011) , convergent evolution, and contrasting rates of morphological and protein evolution (Cherry et al. 1979) show that morphological similarity may not couple to the protein sequence similarity (Stephen Jay Gould 1970) . Inspired by the unfolding of morphological complexity in development, we propose that cellular polarity may be the key for reconciling complexity, robustness, and diversity of organismal morphologies.
During early development the increase in morphological complexity coincides with the progressive polarization of cells -first apicalbasal (AB) polarity and then planar cell polarity (PCP) (Müller and Bossinger 2003; Roignot, Peng, and Mostov 2013; Andrew and Ewald 2010; R. Li and Bowerman 2010 . This theme is ubiquitous across vertebrates and invertebrates: starting from the single fertilized egg cell, first the morula formed by nonpolarized cells turns into the blastula -a hollow sphere of cells with AB polarity. Then, as cells acquire additional PCP, primary headtail axis forms and elongates during gastrulation and neurulation (Loh, van Amerongen, and Nusse 2016) (Figure 1). Because of the optical transparency these stages are particularly prominent in sea urchin. At the morula stage, a lumen in the center is formed and is gradually expanding as cells proliferate and rearrange into the hollow sphere. Next, during gastrulation, a group of cells invaginate and rearrange into a tube that narrows and elongates primarily by cell rearrangement and convergent extension movements (Martik and McClay 2017) . The tube then merges with the sphere at the side opposite to invagination, and as a result, the sphere transforms into a torus. Emerging data suggest that PCP drives both invagination and tube elongation (Nishimura, Honda, and Takeichi 2012; Croce et al. 2006; Long et al. 2015) -a recurring theme in gastrulation across species. Figure 1. Two symmetrybreaking events, gain of apicalbasal polarity and planar cell polarity, on cellular level coincide with the appearance of a rich set of morphologies. Starting from an aggregate of nonpolarized cells (globular symmetry), individual cells can gain apicalbasal polarity and form one or multiple lumens (spherical symmetry). Additional, gain of planar cell polarity allows for tube formation (axial symmetry). Complex morphologies can be formed by combining cells with none, one, or two polarities.
Mutations in PCP pathways produce shorter and wider tubes (OchoaEspinosa, Baer, and Affolter 2012; Saburi et al. 2008; Kunimoto et al. 2017 , somites (Song et al. 2010) , embryos (Gong, Mo, and Fraser 2004) , and also result in neural closure defects (Nishimura, Honda, and Takeichi 2012) . Formation and elongation of the tubes can proceed without cell division and cell death by cells rearranging along the tube's axis, termed convergent extension (CE) (Andrew and Ewald 2010; Tanimizu, Miyajima, and Mostov 2009; Martik and McClay 2017 . While the importance of PCP in gastrulation and tubulogenesis is well established (Andrew and Ewald 2010; Tanimizu, Miyajima, and Mostov 2009; Martik and McClay 2017; Kunimoto et al. 2017 , it is unclear how polarity may control tube stability and dimensions. The bulklumensfolds/tubes transition seen across animal species in early embryogenesis, is also a key feature of the later organ formation. The early stages of organogenesis in liver, kidney, brain, gut, and pancreas are apparently so robust, that they can be recapitulated in vitro , allowing for advanced quantification and manipulation (Little 2017) . The case of pancreatic organoids is interesting as it illustrates an increase of morphological complexity from spheres to folds. Cells in in vitro pancreatic organoids first grow as a bulk and later acquire AB polarity and develop lumens. Depending on the growth conditions organoids develop into a hollow sphere or acquire a complex folded shape (Greggio et al. 2013) . It is currently unknown what drives the transition from sphere to folded state; the two main hypothesis are rapid proliferation or physical pressure by growing into a stiff matrigel.
Is the apparent link between cellular polarity and morphological complexity accidental? Or, could it be that morphological transitions, stability, and diversity are emergent features in a population of proliferating polarized cells? If true, can we identify what drives the transition from lumens to folds and tubes? Why are these stable? Can we predict what controls fold depth, and tube length and width? To answer these questions, we lack a unified theory that could bridge polar interactions between single cells to the global features emerging on the scale of thousands of cells in 3D.
Starting with D'Arcy Thompson's seminal contribution (Thompson and Others 1942) , quantitative models aided in understanding specific morphogenetic events. Among these are invagination (Odell et al. 1981; Rauzi et al. 2015; Polyakov et al. 2014; Hočevar Brezavšček et al. 2012 , primitive streak formation in gastrulation (Newman 2008) , convergent extension (Collinet et al. 2015; Belmonte, Swat, andGlazier 2016) , epithelial folding (Buske et al. 2012; Osterfield et al. 2013; Monier et al. 2015; Murisic et al. 2015 , emergence of global PCP alignment from local cell-cell coupling (Amonlirdviman et al. 2005; Le Garrec, Lopez, and Kerszberg 2006; Burak and Shraiman 2009 , origins of tubulogenesis (Engelberg et al. 2008) , and recently statistical properties of branching morphogenesis (Hannezo et al. 2017) ( Figure S1). However, they are often on either of the two ends of the spectra: those modeling single cells explicitly, often rely on vertexbased approaches and are limited to dozens of cells (Alt, Ganguly, and Salbreux 2017; Misra et al. 2016; Aigouy et al. 2010; Le Garrec, Lopez, and Kerszberg 2006 . To capture the large features spanning thousands of cells, one typically turns to elastic models where AB polarity is implicit and epithelia is presented as a 2D elastic sheet (Hannezo, Prost, and Joanny 2014; Etournay et al. 2015; Hufnagel et al. 2007; Nagai and Honda 2009; Aliee et al. 2012; Nagai and Honda 2001 . We developed a theoretical approach that, with only a few parameters, bridges cellular and organ scales by integrating both types of polarity. A main difference to earlier approaches is that a cell's movement is coupled to how its AB polarity and PCP are oriented relative to each other and relative to neighbor cell polarities. In other words, the adhesion strength between neighbor cells is modulated by the orientation of their polarities. We find that polarity enables complex shapes robust to noise but sensitive to changes in initial and boundary constraints, thus supporting that morphological stability and diversity are emergent properties of polarized cell populations. Lumens, folds, and stable tubes emerge as a result of energy minimization constraints. We make testable predictions on morphological transitions in pancreatic organoids, tubulogenesis, and sea urchin gastrulation. Our approach illustrates the evolutionary flexibility in the regulatory proteins and networks, and suggests that despite differences in proteins between organisms, the same core principles may apply. Cell i and j do not interact if the ij 's midpoint is inside of the Voronoi diagram for cell k (shaded in grey). (C) Cell i and j interact because cell k is further away than the distance ij /2 and ij 's midpoint therefore lie outside of cell k 's Voronoi diagram. See also Figure S2 where we test the sensitivity of our model to the details of the potential and neighborhood assignments.

MODEL
There are three key points that allow us to bridge the scale from cellular level to macroscopic stable morphologies.

(1) Cells are approximated by point particles
Cell-cell adhesion is modeled by repulsive and attractive forces acting between cell centers. This allows a substantially gain in computation time compared to vertex based models where adhesion of cell surface adhesion are explicitly considered (Alt, Ganguly, and Salbreux 2017) . The potential for pairwise interaction between two interacting neighbors, i and j , separated by distance r ij is where the the first term corresponds to repulsion and the second term to attraction (see Figure 2A). For a pair of nonpolar cells the strength of attraction S = 1. β > 1 is the parameter that sets how much longer the attraction range is compared to repulsion. We set β = 5 throughout the paper, but our results and conclusions are consistent for smaller β . The main results are also not sensitive to the exact choice of the potential, thus for example the higher power in the exponential, give qualitatively similar results (see Figure S2C). The potential energy of a cell is the sum of pairwise neighbor interactions (3)

(2) Cells interact with (a subset of Voronoi) neighbors
Interacting neighbors of cell i are selected from a subset of cells sharing a Voronoi surface. The subset is limited to the nearest neighbors j which are closest to the midpoint between i and j ( Figure 2B-C). This constrain effectively corrects for the finite volume associated with point particles and assures that two cells will not interact if the line of sight between their centers is separated by a surface of a third cell. Without this constraint, the macroscopic morphologies collapse. However, our results are robust to replacing the line of sight constraint with full Voronoi and a cutoff distance for attraction force ( Figure  S2D).

(3) Cells -cell adhesion depends on the orientation of polarity
To capture directional adhesion, we set the strength of attraction, S , to be dependent on the relative orientation of the polarities in each of the cells. We assume that AB polarity and PCP are orthogonal, and that the polarities of one cell align with the polarities of it's neighbor cells. Mathematically, we introduce unit vectors, and representing AB polarity, and and representing PCP for cell i and j , p respectively. We set S = λ 1 S 1 + λ 2 S 2 + λ 3 S 3 , and require that λ 1 + λ 2 + λ 3 = 1 , to satisfy the constraint that perfectly aligned cells always have a steady state distance of 2 cell radii (for β = 5). To capture that in an epithelial sheet, AB polarities align parallel to each other and tight adherens junctions form in the plane perpendicular to AB polarity, we introduce the quadruple product This makes two interacting cells with AB polarity maximally attracted ( S 1 = 1) if the two apical sides are next to each other. On the other hand, if apical side of one cell is next to the basal side of another cell, the two cells will be maximally repulsing ( S 1 = 1), see Figure 2A. In case of planar polarization, we define which makes the attraction maximal if the PCP of two cells are parallel to each other and perpendicular to their AB polarities. In addition, we assume that similarly to AB polarity, two cells with PCP are maximally attracted if their PCP are parallel and cells have the same kind of pole (e.g. Vanglenriched) next to each other, We later show that this assumption makes neighbor exchange on a sheet possible and results in CE. However, unlike with tight junctions, there is no biological confirmation for preferred directional adhesion with PCP.
The motion of the cells and their polarities are calculated assuming overdamped dynamics + η , dt dr i = dr i dV i + η , and + η where the and differentiation takes into account the rotation of dt dp i = dp i dV i dt polarity vectors, and η is a random uncorrelated Gaussian noise. In practice, we use the Euler method, and perform the differentiation along the polarity by differentiating along all three cartesian coordinates (see Model in the supplementary material). After each time step, we normalize the updated polarity vectors. The above differentiation does not include the change in partners when neighborhood changes. This is treated as a nonequilibrium step where potential energy can increase (Equation 3). Biologically this is similar to cells spending biochemical energy as they rearrange their neighborhood.
The point particle approximation has been utilized earlier for modeling nonpolar cell adhesion in early blastocyst (Krupinski, Chickarmane, and Peterson 2011) , slug formation in amebae (Dallon and Othmer 2004) , and PCP organization in primitive streak formation (Newman 2008) . The main novelty of our approach is in the dynamical coupling of cell positions and polarity orientations (Equation 6-8).

RESULTS
We have recently introduced effective representation of AB polarity, and showed that it is sufficient for capturing spherical trophectoderm in the early blastocyst (Nissen et al. 2017) . Expanding on that work, we here explore how AB polarity supports diverse yet stable and complex morphologies.

Stable complex shapes emerge from randomly polarized cell aggregates
To test if cellular polarity could enable stable morphologies with lumens and folds as in developing organisms, we first performed a series of tests with AB polarized cells ( Figure 3 and Movie 1).
When starting a bulk of cells with AB polarities pointing randomly, an initial rapid expansion ( Figure  3A-C) stabilizes into a complex topology of interconnected tunnels ( Figure 3C-E). The shape remains unchanged for at least 10 times longer than the initial expansion ( Figure 3C-E). The stability of the shape is illustrated by the time evolution of the average energy per cell ( Figure 3F) that after an initial fast drop converges to a constant value. As expected, this value is higher than the energy of a hollow sphere (yellow dot in Figure 3F) -a configuration obtained if we start with radially, instead of random, polarized cells and preserve radial polarization at all times. This behaviour is not sensitive to the shape of the potential ( Figure S2C) but is sensitive to how the neighborhood is defined ( Figure S2D -F). Rerunning the simulation in Figure 3 with different initial conditions results in a different shape ( Figures S2A and S3). At the same time the shapes are robust to dynamical noise ( Figures S2A-B and S3).
The obtained diversity of metastable topologies suggests a metaphor for evolution of body plans. Diversity and metastability goes hand in hand to create a world of different, yet qualitatively similar shapes. Each could switch into a new distinct geometry by localized changes of a subset of the polarities at the initial state. Thereby, our results support a metascale view of evolution where small modifications can lead to large morphological changes, reminiscent of punctuated equilibrium (S. J. Gould and Eldredge 1993) . , if the polarities are allowed to change dynamically and the noise is high, the resulting shape consists of three nested "Russian doll"like hollow spheres that will never merge due to opposing polarities. See also Figures S2A-B and S3 where we quantitatively study the importance of noise and initial polarities.

The final shapes are robust to noise but sensitive to initial and boundary conditions
To investigate sensitivity to boundary conditions, we consider three cases where polarities are fixed at all times and point either radially out from center of mass ( Figure 4A), radially out from a central axis ( Figure 4B), or pointing away from a central plane ( Figure 4C). As anticipated, the difference in symmetries of boundary conditions result in a sphere, a cylinder, or two parallel planes. At the same time, in these symmetric cases the differences in initial conditions (no boundary conditions) are not sufficient to generate different structures; they all converge to the nested "Russian doll"like hollow spheres ( Figure  4D). In development, this highlights the importance of the neighboring tissues for defining boundary conditions.
Our results thus support the idea that polar adhesion enables stable and robust macroscopic shapes. The closest biological parallels would be the complex luminal morphologies emerging in reaggregation experiments on e.g. Hydra (Seybold, Salvenmoser, and Hobmayer 2016) or in vitro culture of purjunkie brain cells (Muguruma et al. 2015) . Together with our simulations, these experiments highlight how stable and complex topologies can develop in nonproliferating populations from cell rearrangements alone.  Figure S4 for additional measurements on the differences between rapid growth and pressure. In these simulations, two cells do not attract to each other if the angle between their apicalbasal polarity is larger than π/2.

Folding by pressure or rapid proliferation result in different foldmorphologies
Growth is inherent for developing organs (and organisms) and morphological complexity increase with increasing number of cells. In vitro pancreatic organoids, started with a few nonpolarized cells expand about 50 fold in 7 days. The complexity of the shape correlates with organoid size: small tend to be more spherical but the largest 20 % show folded and branched morphologies reminiscent of pancreas (Greggio et al. 2013) .
Transitions from spheres to folded shapes are ubiquitous in 3D organoid systems. Folds are an important part of in vivo organ development, and the composition of cell types in the folded organoids is closer to that in real organs (Greggio et al. 2013) . To date, it is unclear what drives the transition from spheres to folded lumens. One possibility is that it is driven by the mechanical properties of the matrigel that effectively may place the growing organoid under pressure. Alternatively, data from 3D brain organoids suggests that the rapid cell proliferation leads to the emergence of surface folding (Y. Li et al. 2017) .
The simplicity of our tool allows to explore both of these scenarios. To model dividing cells, we pick a random cell from the entire population and introduce a new daughter cell with inherited polarity direction placed in a random location a half cell radius away from the mother cell. This event introduces dynamic perturbation by locally increasing cell density and requires some time to relax back to equilibrium. If proliferation is slow and the time between cell doubling is longer than the relaxation time, the system approaches equilibrium and will expand as a sphere. However, if proliferation is increased, the system will be pushed out of equilibrium and folds will emerge ( Figure 5A, Methods).
As cells divide faster, our simulations predict a transition from a smooth spherical shell to an increasingly folded structure with multiple folds and pronounced local maxima, in line with the observation of brain organoids proliferating at different rates (Y. Li et al. 2017) . In comparison with the model for cortical convolutions by Tallinen et al. (2016) in which folding is a result of expanding cortical sheet adhered to the nonexpanding white matter core, our mechanism does not require a bulk core. Instead folds emerge in a fast expanding sheet when the growth is faster than the global relaxation to dynamical equilibrium.
While we find that the external pressure is not necessary for folding, pressure alone can also drive folding ( Figure 5B, Methods). However, this scenario contradicts the observation that pancreatic organoids can grow as spheres or folded morphologies in gels with same stiffness but different media composition (Greggio et al. 2013) .
In principle, both scenarios may contribute to folding, but visually the fold morphologies are very distinct.
To differentiate between the two, we have quantified the final folded structures in terms of its local minima ( Figure 5, Methods). Our simulations predict that in the pressure driven case, the number of local minima does not exceed a upper limit of about 10-12. In the case of outofequilibrium proliferation, new f olds can continue forming as organoids grow. Increased proliferation causes more and shallower folds. These folds are different than obtained with pressure which causes a shape with fewer but deeper minima. Quantitatively, both the depth and the horizontal extension of the folds are about double as large with pressure than with growth induced folding ( Figure S4). For each value of λ 3 , we initialize 1000 cells on a hollow sphere with PCP whirling around an internal axis. Semimajor axis (dark blue) and semiminor axis (light blue) are measured at the final stage (Methods). Images show the final, frozen state. Throughout the figure, λ 2 = 0.5 and λ 1 = 1 λ 2 λ 3 . The animated evolution from sphere to tube is shown in supplementary Movie 2.

PCP enables convergent extension and robust tubulogenesis
Despite the numerous evidence supporting the role of PCP in tubulogenesis, it remains unresolved whether oriented cell division or the extent of convergent extension controls tube length and width (Karner et al. 2009; Carroll andYu 2012) . It is also debated if it is important for the stability of the tubes, or if it is only important for tube initiation and growth (Kunimoto et al. 2017) .
The simplicity of our approach allows us to address these questions by introducing cell-cell interactions through PCP. This term favours frontrear cell alignment in the interaction potential with o nly two additional parameters: t he strength of the orientational constraint of AB polarity and PCP, λ 2 , and the strength of PCP, λ 3 (see Equations 4-8). For simplicity, we focus on the stability and tube morphogenesis in systems without cell division.
Inducing PCP in a spherical lumen leads to two significant events. First, independent of initial orientation, after some transient time PCP becomes globally ordered and point around a sphere. This arrangement has the lowest energy state. Second, cells start intercalating along the axis perpendicular to PCP orientation, gradually elongating the lumen (Movie 2). The intercalations along the axis continue until the force balance between AB polarity and PCP is restored at a new equilibrium. Thus, our model predicts that the strength of PCP (λ 2 together with λ 3 ) relative to AB polarity (λ 1 ) determines the width and the length of the tube. Note that this result is very different in nature from the tube presented in Figure 4B as both AB polarity and PCP can now reorient in each cell at any time.
These results support the observations that stable tubes can emerge without cell proliferation. In addition, when first the tube is formed, loss of PCP does not lead to cyst formation as recently shown by Kunimoto et al. (2017). However, localized cysts could result if the lumen is initialized with varying strength of PCP along the axis.

Two polarities are sufficient to explain major features of sea urchin gastrulation
Currently invagination in neurulation and gastrulation is understood and quantitatively modeled as a process driven by changes in cell shapes or the mechanical properties of cells with AB polarity (Rauzi et al. 2013; Tamulonis et al. 2011; Misra et al. 2016; Hočevar Brezavšček et al. 2012) . This process is often assumed to be driven by apical constriction and decoupled from the eventual tube formation and elongation. However, emerging data suggests that first, apical constriction may not be necessary (Chung, Kim, and Andrew 2017) , and second, that PCP drives both invagination and tube elongation (Nishimura, Honda, and Takeichi 2012; Croce et al. 2006; Long et al. 2015 . To probe the limits of our approach, we investigated if AB polarity, PCP, and boundary conditions reminiscent of posterior organizer (Loh, van Amerongen, and Nusse 2016) are sufficient to recapitulate the main stages of sea urchin gastrulation: invagination, tube formation, elongation (by CE), and finally merging of the gastrula tube with the pole opposite to invagination site.
In sea urchin, PCP is likely to be induced in a subpopulation of cells placed as a concentric ring at the site of invagination. PCP is necessary for both invagination and tube extension (Croce et al. 2006) . Motivated by this, we assign PCP to all cells in the lower third of the sphere by setting their λ 2 to 0.1 keeping λ 3 = 0 ( Figure 7A, Methods) . Further, we assume that PCP initially points toward the center of mass due to signals from the primary mesenchymal cells, a directional constraint that then forces the AB polarity to reorient. As a result of local optimization, cells start rearranging. A subpopulation of cells flatten the bottom ( Figure 7B ) and bends inward ( Figure 7C ). At this time point, we allow PCP to relax. When PCP has reached its low energy state ( Figure 7D ), we increase λ 2 to 0.4, and set λ 3 = 0.1. The subsequent relaxation to these stronger PCP interactions drives the tube formation ( Figure 7E ) that elongates and eventually merges with the top of the sphere ( Figure 7F and Movie 3 ). In line with experimental observations, the tube elongates due to cells moving into the tube (Martik and McClay 2017) .
The above model leads to a testable prediction: During primary axis formation, local organizing signals (e.g. WNT) (Loh, van Amerongen, and Nusse 2016; J. Wu et al. 2013; Chu and Sokol 2016 orient PCP in cells along the future axis, away from the organizer. As a consequence, we predict that the "wrapping" in neurulation and gastrulation in Drosophila ( Figure S5) vs. "budding" in sea urchin and organogenesis Ewald 2010; Zegers 2014) are outcomes of different boundary conditions. In both cases, we initiated the process without the use of CE (λ 3 = 0), but with different initially imposed symmetries. The sea urchin gastrulation was driven by PCP directed toward a point, whereas neurulation (and similarly configured gastrulation in Drosophila ) was driven by PCP initially directed relative to a line. Thus, we propose that different manifestations of gastrulation reflects different outside constraints on the PCP.

DISCUSSION
Despite the stunning diversity and complexity of morphologies, the same concepts seem to emerge across organismal development. One of them is the link between local, cellular, and global, organs/whole organism, symmetry breaking. We know, from experimental and, to a lesser degree, theoretical work that cellular polarity is essential for forming axis, complex folded sheets, and interconnected tubes ( Figure  S1). What we do not know is why are these shapes so stable, and where do the differences between species and organs come from. To understand the differences, we typically compare genes or gene regulatory networks, thus limiting ourselves to related processes in a few related species.

Main results
To address the origins of morphological diversity and stability across species and organs, we focused on a phenomenological description of polarized cell-cell interactions. This allowed us to bridge local single cell symmetry breaking events to global changes in morphologies spanning tens of thousands of interacting cells. With this tool at hand, we find that with only a few parameters, we can recapitulate the two global symmetry breaking events: formation of epithelial sheets and folds by cells with AB polarity, and emergence of global axial symmetry (tubes) among cells with PCP.
Remarkably, our results show that interactions among AB polarized cells lead to metastable topologies, that after initial relaxation remain indefinitely in their final configuration. The topologies are robust to noise, growth, and local damage. These results may explain how organs and embryos preserve their architecture while growing. Polar cell-cell interactions not only provide clue to the morphological stability, but also point to a simple explanation to the origin of the diversity. We find that the exact morphological details are defined by initial conditions, e.g. initial positions and orientation of polarities, and boundary conditions, e.g. fixed polarities for a fraction of the cells in some time. It is thus tempting to speculate that diverse shapes do not require multiple interacting morphogen gradients, but simply can be a result of differences in initial and/or boundary conditions: as for example presence of yolk cells at start and boundary constraints by vitelline membrane (Schierenberg andJunkersdorf 1992; T. Wu et al. 2010) .
The diversity of shapes and forms is further enriched by a second symmetry breaking event, PCP, oriented perpendicular to AB axis. Within our phenomenological framework, addition of PCP component is simple, and requires only two additional parameters: one favouring perpendicular orientation of AB polarity and PCP within a cell, and another, favouring parallel PCP alignments between neighbor cells. These constraints are the coarsegrained representation of the wellestablished experimental and computational results on intracellular symmetry breaking events and global ordering of planar polarities mediated by cell-cell coupling (Le Garrec, Lopez, and Kerszberg 2006; Amonlirdviman et al. 2005; Wang, Badea, and Nathans 2006 . The first constraint allowed formation of axial symmetry and in combination with AB polarity, stable tubes. The second constraint resulted in cell rearrangements and intercalations consistent with CE typically associated with PCP. The mechanism of the CE in our model is in line with the results by "filopodia tension model" where elongated structures of many cells emerge from local cell-cell interactions in a direction defined by PCP (Belmonte, Swat, and Glazier 2016) .
Combining AB polarization and a local induction of PCP in a subpopulation of cells was sufficient to obtain main stages of sea urchin gastrulation: invagination, tube formation, and elongation through CE as well as merging of the tube with the animal pole at the top of the blastula. By altering the initial and boundary conditions to resemble those of neural plate folding, we, in line with the experimental observations, could obtain transition from an neuroepithelial sheet to a sheet and a tube.

Testable predictions
In addition to our conceptual findings, we arrived to three testable predictions. First, we predict that the two potential mechanisms behind the emergence of folds in pancreatic organoids -matrigel resistance and rapid, outofequilibrium, cell proliferation -will result in distinct morphologies. Our results suggest that in case of rapid proliferation, the growing structure will develop many folds close to the surface which later tend to deepen. In contrast, external pressure causes fewer and much deeper folds that form early during growth. Visual inspection of published morphologies seems to support the outofequilibrium growth (Greggio et al. 2013; Y. Li et al. 2017 . Our model predicts that quantitative measurements of the fold depth and length relative to the size of the growing structure will discriminate between the alternative hypothesis. This can be done in in vitro organoids by either phase or confocal fluorescence microscopy (Greggio et al. 2013; Y. Li et al. 2017) . In addition, distributions of cell shapes can be used to asses the outofequilibrium growth in 3D organoids (Cerruti et al. 2013) .
Our second prediction is that in case of tubes formed by nonproliferating cells, the length and width of the tubes are controlled by the relative strength of AB polarity and PCP. This result calls for quantification of adhesion proteins along the AB and PCP axes. In PCP mutants with shorter and wider tubes one would expect higher fraction of "ABadhesion" relative to "PCPadhesion".
Our third testable prediction is on the conditions differentiating between tubes forming perpendicular (e.g. sea urchin gastrulation) or parallel (as in Drosophila gastrulation or neurulation) to the plane of epithelium. We predict that the outcome will be defined by the orientation of PCP in the invaginating region and at the boundaries set by e.g. WNT organizing signals. Recent development in imaging localization of PCP complexes in single cells (J. Wu et al. 2013; Chu and Sokol 2016; Minegishi et al. 2017; Habib et al. 2013) allows monitoring localization of PCP complexes, and thus PCP orientation, in individual cells. By placing WNTsoaked (Habib et al. 2013) beads or WNTsecreting cells (Chu and Sokol 2016) one can vary PCP orientation in the cells at the epithelial boundary facing WNT and test for the direction of the tube formation.
Our results open for a series of biological generalizations both in development and diseases. On one hand, we now may be able to explain and unify the apparently very distinct morphological transitions during gastrulation in flies, frogs, fish, mice, and humans by accounting for different initial and boundary conditions. Our model predicts that a moderate change in expression of polarities during some critical evolutionary stages could lead to widely different final morphologies. Thereby, development driven by cell-cell polarity interactions could provide major morphological transitions from local and transient modulations in polarity. This suggests how evolution may bypass larger scale protein evolution to obtain dramatic changes in body plans.
On the other hand, it becomes possible to think of gastrulation, neurulation, tubulogenesis, and organogenesis as the same class of phenomena, where the orientation of the tube is guided by local organizers, and lengths/widths of the tubes are determined by the relative strength of AB polarity and PCP. At the same time, there is an emerging view that wound healing and cancer are local perturbationse.g. local loss of cells, dysregulation of cell polarities (MartinBelmonte and PerezMoreno 2011) , proliferation, or autonomously induced organizing signals -of otherwise conceptually the same developmental processes (Humphries and Mlodzik 2017) . The power of our here presented tool is that it allows to quantitatively address these hypothesis through predictive models for the dynamics of many cells that interact through combinations of AB polarity and PCP.

METHODS
Throughout the paper, we use the Euler method to integrate the ordinary differential equations stated in the Model section with dt set to 0.1 or 0.2. Lower dt values will give qualitatively similar results but with increased simulation time. Higher dt values, will result in a collapse of the presented morphologies.
The noise parameter η = 10 4 where nothing else is stated. Lower noise values will give more smooth simulations, while η on the order of 10 1 will result in collapsing shapes (see also Figure S2A-B).

Proliferation rate in growing organoids
In Figure 5, the number of cells, N , at a given time, t , is given by N = 200 exp( Kt ) where K is the growth rate. For reference, the timescale in Figure 5 is a factor thousand slower than the timescale in Figure 3.

Modeling resistance from matrigel
In Figure 5, we model resistance from the matrigel by imposing a surface force pointing towards the center of mass. The potential of the pressure in the growth medium is given by V M = Pr 2 / r max where P is the stiffness of the medium, r is distance from the center of mass, and r max is the distance to the cell that is the furthest away from the center of mass. The resulting force will be constant in time at the periphery. Thus, all cells on a growing sphere will be exposed to a force of equal size. However, cells that end up deep inside a folded morphology will experience weaker resistance.

Quantification of the local minima
In Figure 5, the number of local minima is defined as the number of cells that do not have any neighbor cells that are closer to the center of mass than themself, and at the same time have an average angle between their AB polarity and their neighbor cells displacement vector that is less than π/2.

Measuring the tube length and width
In Figure 6, the semiminor and semimajor axes correspond to the halfwidth and halflength of the tubes, respectively. As cells on opposite sides of a tube have AB polarity pointing in opposite directions, we approximate the semimajor and semiminor axes, by finding the half of the maximum and minimum distance between two cells with AB polarity pointing in opposite directions.

Modeling cells with different polarities
In Figure 7, each cell is assigned a specific value of polarity strengths ( λ 1, i , λ 2, i , and λ 3, i ). We define the mutual interaction strength between a pair, i and j , of cells in Equation 4 with different polarity strengths by setting λ 1 = max( λ 1, i , λ 1, j ), and λ 2 = min( λ 2, i , λ 2, j ) as well as λ 3 = min( λ 3, i , λ 3, j ). This choice makes sure that two neighbor cells interact with a force with equal magnitude but opposite sign. The maximum for λ 1 is chosen to keep a steadystate distance of 2 cell radii between perfectly aligned cells. The pairwise average between polarity strengths will give qualitatively similar results.

ACKNOWLEDGEMENT
We would like to thank Anne GrapinBotton for insightful discussions on organoids.

DECLARATION OF INTERESTS
The authors declare no competing interests. Figure 1) Figure S1. Overview of the existing literature on models addressing specific developmental events discussed in our work. For more references on vertex models see Alt et al. (2017). Figure S2. Dependence on noise, the shape of the physical potential, and the interaction partners.

Figure S2 (Related to Figures 2 and 3)
(A) By applying the neighborhood function as shown in Figure 2, the system reaches a metastable state with high noise ( η = 10 2 ) that is comparable to the state obtained under (B) low noise ( η = 10 5 ). (C) Changing the shape of the potential to the short range potential written in Equation 2, the system unfolds and reaches a metastable state ( η = 10 5 ). (D) Full Voronoi interactions with a cutoff does also lead to a metastable state, although a few cells might lose interaction with the majority of cells (cutoff at 3 shown). (E) Simple cutoff (and no Voronoi) does not result in metastable topologies but broken sheets on top of each other (cutoff at 2.5 shown). (F) When all cells always interact with their six nearest neighbors, the system breaks up and no metastable state is obtained. The initial positions and polarities are identical for all six simulations. Figure 3) Figure S3. The final shapes are more sensitive to initial polarities than to noise. (A) The pairwise distance between cells for three systems with identical initial polarities but different noise and three systems with identical noise but different initial polarities. (B) For the same set of aggregates the angle between the pairwise polarities is calculated. The initial positions are the same for all systems. Each system has 8000 cells. Cells are pairs if they were initiated with identical position. Here, the noise level is η = 10 5 . Twosample KolmogorovSmirnov tests showed p < 0.001 statistical significance (marked by *). Comparing noise levels give similar results as comparing noise seed. The initial polarities are random like in Figure 3. Figure 5) Figure S4. Organoids grown in a matrigel have deeper and longer folds compared to organoids grown with rapid cell proliferation. We fill the surface of the organoids with water until halfway between the maximum and minimum radius of the system. Then we measure the relative depth and circumference of these 'lakes'. (A-C) Deepest point of the 'lakes' (folds) relative to the water level. The probability of having a lake at a given depth is normalized to the number of 'lakes'. (D-F) Length of the 'lakes' relative to the entire circumference at this same level. Length of a lake is defined from the angle between the two cells at lake shore that are the furthest away from each other. Pressure and proliferation rate increase from upper to lower panels. Twosample KolmogorovSmirnov tests showed p < 0.001 statistical significance (marked by *). The shown histograms are for the 16.000 cell stage, which compares to the dark blue line in Figure 5. Comparing the initial stage to the final stage, the overall direction of PCP in the plate is conserved while in the tube PCP goes around an internal axis. For this simulation, we set λ 2 = 0.5 and λ 3 = 0. Turning on convergent extension (λ 3 ) at the final stage would allow for elongating the system along the axis going through the tube and narrowing it in another direction. The concept is similar to gastrulation in Drosophila and sea urchin. The latter is seen in Figure 7.

Movie 1 (Related to Figure 3)
An aggregate of 8000 cells all with apicalbasal polarities initially pointing in random directions develops into a final stable complex morphology. During the simulation the polarities and positions are updated dynamically with equal speed. λ 1 = 1, and there is no planar cell polarity. Panel A shows the entire system, panel B shows half of the system viewed from the inside, while panel C shows half of the system viewed from an angle slightly above and slightly from a side. The color scheme is as described in Figure 3. Figure 6) Model of tubulogenesis. A spherical lumen consisting of 1000 cells with apicalbasal polarity pointing radially out gets planar cell polarity (PCP). In equilibrium, PCP will curl around an internal axis. Depending on the relative strength of the polarities the sphere will elongate and transform to a tube with a given length and width. This movie illustrates the most extreme scenario in Figure 6 (λ 1 = 0.41, λ 2 = 0.5, and λ 3 = 0.09). Panels are similar to Movie 1. The color scheme is as described in Figure 3.

Movie 3 (Related to Figure 7)
Model of sea urchin gastrulation. Starting from a hollow sphere of 1000 cells with apicalbasal polarity pointing radially out. By giving a fraction of cells planar cell polarity (PCP) pointing towards the center, the bottom flattens and invaginates. The tube forms, elongates, and connects with the upper part of the blastula when PCP has reached equilibrium curling around the axis. For more details see the caption to Figure 7. Panels are similar to Movie 1. The color scheme is as described in Figures 3 and 7.