## Abstract

The growth of a well-formed epithelial structure is governed by mechanical constraints, cellular apico-basal polarity, and spatially controlled cell division. Here we compared the predictions of a mathematical model of epithelial growth with the morphological analysis of 3D epithelial structures. In both in vitro cyst models and in developing epithelial structures in vivo, epithelial growth could take place close to or far from mechanical equilibrium, and was determined by the hierarchy of time-scales of cell division, cell–cell rearrangements, and lumen dynamics. Equilibrium properties could be inferred by the analysis of cell–cell contact topologies, and the nonequilibrium phenotype was altered by inhibiting ROCK activity. The occurrence of an aberrant multilumen phenotype was linked to fast nonequilibrium growth, even when geometric control of cell division was correctly enforced. We predicted and verified experimentally that slowing down cell division partially rescued a multilumen phenotype induced by altered polarity. These results improve our understanding of the development of epithelial organs and, ultimately, of carcinogenesis.

## Introduction

Epithelial morphogenesis is a complex process involving cell divisions, cell–cell and cell–ECM adhesion, cell migration, cell shape changes, and apoptosis, and represents a fundamental step in organogenesis. Indeed, these features are fundamental for the correct functioning of the tissue in terms of proliferation, survival, and differentiation. Aberrant epithelial architecture is most frequently found in the pathogenesis of epithelial tumors, and architectural patterns have been used for decades by pathologists to diagnose and classify carcinomas. The study of morphogenetic processes leading to the formation of epithelial tissues can thus be used to gain a better understanding of the development of epithelial organs and of carcinogenesis.

In vitro biological models have been successfully used to reproduce some of the key events involved in epithelial morphogenesis, and represent a fundamental tool to dissect the molecular cascade of events leading to the formation of tissues (O’Brien et al., 2002; Debnath and Brugge, 2005).

Cystogenesis is one of the best studied examples of epithelial morphogenesis in vitro (McAteer et al., 1986; O’Brien et al., 2002) and is considered to be a prototype for the development of several spherical structures encountered in vivo, such as acini, follicles, ampullae, and alveoli. Cysts are spherical monolayers of epithelial cells enclosing a central lumen (McAteer et al., 1986). Cells within cysts are connected by specialized junctions and cell–cell adhesion structures lying in the basolateral sides, whereas a strong apicobasal polarization characterizes the external surface, contacting the ECM, and the apical surface, facing the lumen. The correct architecture and the formation and maintenance of the lumen are crucial for normal cyst morphology and are altered in several common human diseases such as polycystic kidney disease (Boletta and Germino, 2003), hypertension (Iruela-Arispe and Davis, 2009), and many epithelial cancers, such as prostate carcinomas or preinvasive epithelial lesions (Debnath and Brugge, 2005).

Despite the specificity inherent to diverse types of tissues, recent findings support the idea that the formation of several spheroidal epithelial structures could be generated by common mechanisms, and that shared features underlie the appearance of aberrant phenotypes (Datta et al., 2011).

The first general key aspect involved in the process of cyst growth is the mechanics of cell contacts. Epithelial cells are physically connected to the ECM via integrin receptors (O’Brien et al., 2001), and neighboring cells are tightly linked by cell–cell junctions via adhesion receptors, such as cadherins and nectins (Harris and Tepass, 2010). Cell shape variations are caused by local deformations of the cortical actomyosin network. The cumulative effect of differential cell–matrix and cell–cell adhesion processes and of cortical elasticity can be described in terms of interfacial tensions, which have been shown to be the driving force behind tissue formation in several biological models (Käfer et al., 2007; Lecuit and Lenne, 2007; Manning et al., 2010).

A second aspect involves apico-basal polarization and the de novo generation of a luminal space. Luminogenesis proceeds through a coordinated series of molecular events starting with the exocytosis of apical membrane proteins (such as Crumbs3a [Crb3], podocalyxin [PCX], and Mucin 1 [Muc1]) to the cell surface, leading to the formation of the nascent lumen in a region termed the apical membrane initiation site (AMIS; Schlüter et al., 2009; Bryant et al., 2010). Similar structures have been observed during vascular lumen formation in developing mouse aorta (Strilić et al., 2009) and during neural rod formation in zebrafish (Tawk et al., 2007). After the formation of the AMIS, an asymmetric distribution of the phosphoinositides PIP2 and PIP3 is established (Shewan et al., 2011). In particular, the apical region is enriched in PIP2 and PTEN, whereas PIP3 is localized exclusively to the basolateral membrane. The AMIS matures to form a preapical patch (PAP), and eventually a lumen expands (Martín-Belmonte et al., 2007; Ferrari et al., 2008; Bryant et al., 2010; Datta et al., 2011).

A third aspect is the spatial control of cell division. The apico-basal polarization of specialized molecules such as PIP2, PTEN (Martín-Belmonte et al., 2007), Cdc42 (Jaffe et al., 2008), the Cdc42-specific exchange factors Tuba (Qin et al., 2010) and Intersectin-2 (Rodríguez-Fraticelli et al., 2010), Par3 (Hao et al., 2010), aPKC (Qin et al., 2010), and LGN (Zheng et al., 2010) restricts the formation of the mitotic spindle to the monolayer plane, ensuring the stability of the polarized architecture during subsequent cell divisions.

Remarkable progress in understanding the molecular cascade giving rise to the formation of a cyst has been made thanks to 3D cultures of MDCK cells and other similarly well-polarized epithelial cell lines, which provide an excellent model for epithelial morphogenesis in vitro and give valuable insights into the basic processes driving the formation of the correct tubular architecture and its disruption in cancer (Debnath and Brugge, 2005; Datta et al., 2011; Lo et al., 2012). Nonetheless, the precise connection between molecular and mechanical cues and the resulting dynamics is still largely unknown. In particular, it has not been demonstrated whether the simple aforementioned mechanistic principles are sufficient to reproduce cyst architectures. Furthermore, although it is clear that misorientation of the spindle is an ingredient leading to the aberrant multilumen phenotype (Zheng et al., 2010; Engelberg et al., 2011), it is not clear how and to what extent cyst architecture and orientation of cell division are related.

To shed light on these issues, we developed a mathematical model based on the phenomenology of cell–cell and cell–matrix interactions that reproduces the experimentally observed phenotypes in both the normal and the aberrant cases. Our model indicates that cystogenesis occurs through a sequence of states that are far away from mechanical equilibrium, and provides evidence that the appearance of the aberrant multilumen phenotype is determined by misorientation of the spindle only in conjunction with strongly out-of-equilibrium dynamics. We verify our predictions by performing computer-assisted, full 3D reconstruction at different stages of growth of MDCK cysts and of 3D epithelial structures derived from mouse embryos. We demonstrate that both cysts in vitro and developing epithelial structures in vivo share the same nonequilibrium phenotype. Our model predicts that equilibrium can be favored by slowing down cell division, by altering cell–cell and cell–ECM interactions and cortical contractility, or by increasing cell motility. Moreover, we demonstrate experimentally that slowing down cell division partially rescues a multilumen phenotype, thus confirming our predictions. We provide experimental evidence that inhibition of ROCK activity drives the appearance of equilibrium configurations in MDCK cysts. Finally, our model suggests that an increase in cell duplication rates, a typical hallmark of hyperproliferative diseases like cancer, could alone be responsible for the appearance of aberrant phenotypes, independent of the control of spindle positioning.

## Results

### Theoretical modelling of cyst morphogenesis predicts a crucial role for cell duplication and mechanical relaxation

The geometric structure of epithelia can be described in terms of the mechanical forces that result from cell–cell and cell–ECM interactions. The action of the actomyosin cortex can be described by an overall surface tension that acts so as to minimize the total cell surface (Lecuit and Lenne, 2007; Manning et al., 2010). Conversely, basal and lateral adhesion receptors favor the increase of contact surfaces (Käfer et al., 2007; Lecuit and Lenne, 2007; Manning et al., 2010). These ideas have been successfully used to quantitatively describe the geometrical configuration of ommatidia (Hilgenfeldt et al., 2008) and the growth of wings (Classen et al., 2005; Hufnagel et al., 2007; Lecuit and Lenne, 2007; Farhadifar et al., 2007) in *Drosophila melanogaster* and the spreading of murine sarcoma cell aggregates in vitro (Douezan et al., 2011).

We develop a simple theoretical model of a growing cyst by extending the reasoning presented in Manning et al. (2010) to a 3D context. Such a model, under the assumption of equilibrium dynamics, predicts that the observed configurations of a growing cyst are those that minimize mechanical energy (see Materials and methods and Eqs. 1 and 2) containing contributions from both elastic and adhesive processes. Indeed, minimizing this energy under the assumption of spherical symmetry, it is possible to find analytically the equilibrium configuration for a cyst-like cell aggregate of any number of cells (see Materials and methods). This equilibrium configuration has a central lumen and a soccer ball–like structure of the kind described in Cox and Flikkema (2010). Under more general conditions, an analytical solution is impossible, and it is necessary to resort to computer simulations. To this aim, we implemented a dynamical Potts model (Potts, 1952; Wu, 1982; Graner and Glazier, 1992; Glazier and Graner, 1993) with energy given by Eqs. 1 and 2. Within this model, random motility of cells plays a role analogous to temperature, and the associated thermal energy is always much smaller than the typical values of mechanical energy. (A more detailed description of the computational model is given in the Materials and methods section.)

The dynamics governed by Eqs. 1 and 2 depend crucially on the mechanical relaxation timescale *T*, i.e., on the time needed to reach global equilibrium, and on the characteristic duplication time τ. Whenever the ratio χ = *T*/τ of the relaxation time to the division time is much smaller than unity (χ ≪ 1), cystogenesis proceeds through a sequence of equilibrium states. Conversely, when χ is much larger than unity (χ ≫ 1), the next division round takes place before significant mechanical adjustments occur. In the latter case, the dynamics are constantly kept out of equilibrium by subsequent cell divisions. A deeper understanding of the problem comes from the observation that constrained systems often display a spectrum of relaxation times associated with a complex energy landscape with multiple local minima (Mezard et al., 1987). For the present purposes it will be sufficient to distinguish between the fast timescales that govern mechanical rearrangements without topology changes, such as the relaxation of the basal surface to a near-spherical shape, and the much slower timescales that involve a passage over energy barriers, such as the exchange of cell–cell neighbors. Such global topological rearrangements are associated with the longer timescale *T*.

Assessing the equilibrium or out-of-equilibrium nature of cystogenesis by means of a direct measurement of relaxation timescales is a daunting task. Here we propose a simpler method that allows unambiguous experimental discrimination between the two regimes. Assuming that the cyst is a monolayer and that spherical symmetry is preserved, the two opposite limits of equilibrium and strong deviation from equilibrium yield quite different predictions for the distribution of cell–cell contacts that characterizes the topology of the spherical monolayer.

If mechanical equilibrium holds, the problem of energy minimization with given cell volumes can be mapped on the problem of tiling a sphere with polygons of a given surface. This problem has been solved by a combination of mathematical arguments and extensive computer searches (Cox and Flikkema, 2010). For a cyst composed of several cells *N* > 13, the equilibrium topology is given by the “soccer-ball” configuration made of 12 pentagons and *N* − 12 hexagons.

In the opposite regimen mechanical equilibrium is never reached. In this case, for a 2D monolayer, our model predicts that the probability distribution of cell–cell contacts should reduce to that found in Gibson et al. (2006) and Dubertret et al. (1998), as mechanical relaxations become negligible compared with cell division. Indeed, the models described in Gibson et al. (2006) and Dubertret et al. (1998) consider the topological distribution of polygonal tilings emerging as a sole effect of random cell divisions, without taking into account mechanical relaxation. This nonequilibrium topological distribution is consistent with the observed topologies for a large class of planar metazoan epithelia (Lewis, 1928; Rivier et al., 1995; Dubertret and Rivier, 1997; Dubertret et al., 1998; Gibson et al., 2006; Li et al., 2012; Fig. S1 N). The same arguments hold unchanged for curved epithelial structures, provided that the monolayer structure is preserved and that the number of composing cells is large.

Our numerical simulations qualitatively and quantitatively reproduce the theoretically predicted cyst geometries. Indeed, for χ ≫ 1, the nonequilibrium distribution emerges, as shown in Fig. 1 (A and B), and the corresponding topological distribution also corresponds to the predictions of Gibson et al. (2006) in this limit case (Fig. 1 D). In the opposite regimen (χ ≪ 1), as shown in Fig. 1 C, the basal surface topology complies with the minimal energy tessellation of the sphere (Fig. 1 D, inset).

Thus, the theoretical model predicts two opposite scenarios: if cysts grow under mechanical equilibrium, their surface topology has to be that of a soccer ball (Cox and Flikkema, 2010). Conversely, if cysts grow out of mechanical equilibrium, their surface topology has to be the one predicted in Gibson et al. (2006).

### Cell–cell contact topology unveils the out-of-equilibrium nature of epithelial morphogenesis

To test whether cysts grow under mechanical equilibrium or not, we segmented confocal sections of experimental MDCK cysts and reconstructed the complete 3D structure of several cysts (see Fig. 1, E and F). We measured the topology of the cell aggregate by computing the neighbor map on the outer surface of the cysts. As shown in Fig. 1 D, the topological distribution of MDCK cysts is very close to the one given in Gibson et al. (2006) and definitely different from the equilibrium distribution (Fig. 1 D, inset). Therefore, we conclude that in real cysts, cell division is rapid with respect to mechanical relaxation, i.e., real cells are growing in the χ ≫ 1 regimen. Experimental cystogenesis occurs out of mechanical equilibrium. This theoretical conclusion is in agreement with existing data on real-time imaging of cyst growth showing that, both in cysts and in 2D epithelial layers, cell rearrangements (i.e., neighbor exchange) are rare (Gibson et al., 2006; Pearson and Hunter, 2007; Puliafito et al., 2012).

The direct experimental measurement of relaxation times is practically challenging, as it is difficult to isolate relaxation from other naturally occurring processes in the cyst, such as cell division or apoptosis. However, simulations allow us to easily perform such experiments numerically, and the details of out-of-equilibrium growth can be investigated. To this aim, we have simulated the growth of a cyst up to a size of a given number of cells, *n*_{cells}, blocked cell division, and let the cyst relax mechanically. Relaxation toward the equilibrium topology (Fig. 2, A and B) is a very slow process, as several energetically unfavorable random cell movements have to take place one after another to overcome the many local energy barriers that separate the original nonequilibrium distribution from the globally energy-minimizing configuration (Fig. 2 A, insets). The time needed to reach the predicted equilibrium configuration increases exponentially with the number of cells in the cyst, as shown in Fig. 2 C. This finding confirms that cysts evolve in a rugged energy landscape and are geometrically frustrated systems (Mezard et al., 1987). Geometrical frustration means that local moves that would, by themselves, tend to minimize energy are actually forbidden or strongly discouraged by the constraints acting on the system. As a result, global rearrangements are required, with ensuing long times for relaxation to equilibrium.

Quasi-equilibrium growth can be obtained if the cell division time is significantly slowed down (Fig. 2 D), as this allows mechanical relaxation to take place between subsequent cell divisions.

Although the typical nonequilibrium distribution described in Gibson et al. (2006) has been found in a variety of biological contexts, it has never been proved whether it is relevant in the context of developing epithelial structures in mammals. To test this hypothesis, we studied the topological distribution within tubular epithelia in developing mouse embryos. To this aim, we performed whole mount staining of lung and kidney samples and imaged them on a confocal microscope. After segmenting the data, the topology was calculated for rounded cyst-like alveoli from the lung (Fig. 3 A), the tubular ureter (Fig. 3 B), and the tips of the ureteric tree (Fig. 3 C) from the developing kidney. The analysis of the ureteric tree tips is a particularly relevant comparison to the in vitro cyst data, as MDCK cells were derived from an adult kidney. We found that the topological distributions from these tissues (Fig. 3 D) are indistinguishable from those found in an MDCK cyst (see Fig. 1 and Table S2). Therefore, the topological distribution is a universal fingerprint of nonequilibrium dynamics in cyst-like and tubular epithelial tissues both in vitro and in vivo.

### Relaxation of cortical tension in epithelial cysts favors equilibrium topology

Our theory predicts that cysts could be brought to quasi-equilibrium conditions by slowing down proliferation, thus allowing more time for topological relaxation, or by favoring the passage over topological barriers that separate distinct configurations.

The Rho-associated protein kinase (ROCK) inhibitor Y-27632 is known to relieve cortical tension (Matthews et al., 2006) through inhibition of actomyosin contractility, and we therefore expect it to ease topological rearrangements and induce the appearance of equilibrium phenotypes. Indeed, theoretically, ROCK inhibition would lower the overall mechanical relaxation time *T*, making χ smaller.

To verify the validity of our theoretical prediction, we compared the topology of growing cysts untreated and treated with ROCK inhibitor at different times from seeding, shown in Fig. 4 A. First, we evaluated the effect of ROCK inhibition on cysts by using GFP–myosin light chain (GFP-MLC) as an indicator for cortical tension (Yu et al., 2003). As shown in Fig. 4 B, MLC-GFP is concentrated at the basal surface in untreated cysts (a quantification is given in Fig. S5). When Y-27632 is applied from the first day of the experiment, MLC is no longer accumulated at the basal surface, and instead is diffuse in the cytoplasm, which is consistent with an alteration of cortical tension.

To compare topologies in treated and untreated cysts, we digitally reconstructed >50 cysts, and, for each cyst, we measured the topology of the outer surface and the number of composing cells. Representative images of experimental cysts are shown in Fig. 4 A and the complete set of reconstructed cysts is shown in Fig. S2. The comparison of theoretical predictions and experimental data for a set of cysts is summarized in Fig. 5 A. Cell–cell contact statistics of ROCK-inhibited cysts are significantly closer to the soccer-ball topology, which is the hallmark of equilibrium, as shown in Fig. 5 (B–D). Conversely, as expected, the contact statistics of control cysts are closer to the out-of-equilibrium topology.

To quantify the degree of equilibrium, we computed for each reconstructed cyst the normalized topological distribution and calculated the distance from equilibrium and nonequilibrium distributions (see Materials and methods). The distances from equilibrium of cysts under different treatments at 1–5 d from seeding are shown in Fig. 5 C.

According to our theoretical predictions, the time needed to reach mechanical equilibrium configuration depends on the number of cells in the cysts rather than on time (see Fig. 2 C). Thus, we expect the effect of ROCK inhibition to depend on the number of cells composing the cyst. Indeed, as shown in Fig. 5 (B and C), Y-27632 is very effective on cysts composed by *n* < 40 cells, and its effect becomes weaker for larger aggregates. Note that for both treated and untreated samples, cysts containing a smaller number of cells are generally closer to equilibrium than bigger cysts, as expected.

Treatment with Y-27632 has been reported to slow down cell divisions (Ishizaki et al., 2000). According to our theory, cell division slowdown can concur to relaxation toward the equilibrium configuration. To separate the contributions to this effect of cortical tension relief and cell division slowdown, we compared ROCK-inhibited cells with cells treated with the DNA synthesis inhibitor Aphidicolin. We performed experiments according to a variety of protocols, i.e., supplying the drugs in different combinations either from the first day of the experiment or from the third (Fig. 5 C). The observed slowing down of cell divisions induced by either Y-27632 or Aphidicolin is of a similar order (Fig. 5 E). However, topology is not significantly altered by Aphidicolin, as documented in Fig. 5 D. To explain this apparent discrepancy between our theoretical predictions and the experimental data, we recall the quantification of the effect of cell division slowdown on cell-contact topology (see Fig. 2 D). Indeed, the changes in the duplication time required to entail a significant effect on the relaxation toward equilibrium (on the order of 10-fold) are much larger than those experimentally accessible (for example, 30% to 40% in Fig. 5 E). According to our theoretical model, the duplication time is indeed independent of the number of cells, whereas the relaxation time grows exponentially with the number of cells (Fig. 2 C).

In summary, the observed equilibrium phenotype of Y-27632–treated cysts can most likely be attributed to the release of cortical tension induced by ROCK inhibition, rather than to an effect on proliferation.

### Multi-luminal cysts result from the combined effect of constrained out-of-equilibrium dynamics and disruption of mitotic spindle control

The nonequilibrium character of cyst growth manifests itself with particular evidence in the dynamics of luminal spaces. It is well known that cysts grown with cells that are defective in establishing a proper direction for the cleavage plane are characterized by several, interspersed luminal spaces and the lack of a clear monolayer structure, as shown in Fig. 6 (A and B; Martín-Belmonte et al., 2008). In vivo spherical epithelial structures equipped with such a phenotype would be severely impaired. A stringent control of the mitotic spindle orientation in the apico-basal direction therefore seems necessary. To test how sensitive the single-lumen phenotype is to perturbations in the alignment of the cleavage plane, we performed a large number of numerical simulations with decreasing control on the division geometry (Fig. 6, C–H). Indeed, our results indicate that larger errors in spindle positioning are associated with a dramatic increase of the emergence of aberrant phenotypes (Fig. 6 C, inset). To understand the dynamic origin of this phenomenon, we observed that for large proliferation rates, the process of lumen coalescence is bound to take place far from equilibrium, as new microlumens may be created at each cell division, and the number of new cell divisions per unit of time grows proportionally to the number of cell in the aggregate. In cysts with disrupted cell polarity, multiple lumens compete as aggregators of newly formed microlumens, as each new microlumen can merge with one of several different lumens in the cyst. If the typical time of lumen coalescence is small with respect to the cell division time, the monolumen phenotype is attained. Conversely, when coalescence is slow compared with cell division, the monolumen phenotype may never be reached. As normal MDCK cysts most often show a single lumen (Bryant et al., 2010), we expect the typical coalescence time to be faster than the cell division time.

To further corroborate these observations we performed the following in silico experiment: we simulated the out-of-equilibrium growth of a cyst up to a given size, *n*_{cells}, with disrupted orientation of mitotic spindle. The obtained multiluminal cyst was then relaxed, with cell divisions blocked. We measured the time needed to reach the monolumen phenotype as shown in Fig. 6 F. Such times are faster than those required to reach global equilibrium (and in particular a soccer ball–like configuration), as shown in Fig. 2, which suggests the absence of significant energy barriers. The total time *T*_{m} to reach the monolumen increases with the number of cells, as a larger number of microlumens are generated. In our simulations, smaller lumens were observed to coalesce toward larger lumens by sliding along cell–cell interfaces, as depicted in Fig. 6 G. The sliding movement can take place without appreciably changing the global surface energy (Eq. 2), as the total area of cell–cell contacts does not change. We verified this prediction by simulating the coalescence of a single microlumen into a central lumen and by measuring the characteristic sliding time as a function of surface tension. Our results indicate that, indeed, coalescence time does not depend on surface tension (Fig. 6 H, inset). This is compatible with the experimental observation that ROCK inhibition does not have a significant effect on the speed of lumen coalescence (Fig. S1 P). It is worth observing, however, that although the sliding of microlumens does not imply a significant energy expenditure, its final coalescence into a larger lumen provides an energy advantage, as a more symmetric and energetically favorable distribution of cellular interfaces is realized (Fig. 2 A, inset). The coalescence of lumens is therefore an energetically favored process.

Theoretically, as in the case of the topology of cell contacts, slow cell division compared with lumen coalescence yields equilibrium growth, i.e., monoluminal phenotype, whereas fast cell division favors the accumulation of multiple lumens. To better characterize the role of the proliferation rate, we simulated the growth of a polarity-disrupted cyst with increasingly long division times, starting from a single cell. In normal cysts, lumen coalescence times are typically faster than cell division times, thus allowing single lumen configurations to be reached after each cell duplication. This is also observed in experimental cysts (Fig. 1, E and F). For cysts that have faster cell proliferation, this ceases to be true, giving rise to multiple lumen phenotypes (Fig. 6 H).

To experimentally validate our theoretical predictions, we quantitatively characterized the maturation and dynamics of lumens in MDCK cysts (Fig. S1 P). In control cysts, the incidence of single lumens at day 2 is ∼50%, while it reaches 80% by day 5. As monoluminal cysts become more frequent, the number of immature and multiple lumens decreases.

To verify whether lumen dynamics are actually independent of cortical tension, as suggested by our simulations (Fig. 6 H, inset), we generated multiluminal cysts by using an aPKC pseudosubstrate inhibitor (aPKC-PSi). This inhibitor contains a pseudosubstrate for atypical PKC that has been myristoylated to be cell permeable and alters the correct orientation of cell division (Qin et al., 2010). Cells were then treated with aPKC-PSi or Y-27632, and the effect on lumen dynamics was quantified. As shown in Fig. S1 P, the treatment with Y-27632 has no significant effect on the fraction of monolumen versus multilumen cysts.

We then verified the effect of cell division on lumen dynamics. To this aim, we measured the frequency of the monoluminal phenotype in cysts with disrupted cell polarity where the proliferation rate was slowed down by Aphidicolin (Fig. 6 I). Interestingly, the slowing down of cell divisions induced a significant increase in the frequency of the monolumen phenotype (p-value: 10^{−5}). Although the control had 73.4 ± 3.8% monolumens, incubation with aPKC-PSi for 4 d decreased single lumens to 52.2 ± 6.1%. Slowing down of division by including Aphidicolin largely rescued the formation of single lumens, bringing it back to 66.1 ± 1.2%.

In conclusion, slowing down cell division partially rescues the multilumen phenotype induced by misoriented cell division.

## Discussion

Epithelial tissue growth is a complex process involving both biochemical signaling cascades and mechanical forces. The coordination between these two aspects is crucial for the correct development and maintenance of a well-ordered architecture, whereas aberrant phenotypes are the hallmark of life-threatening human diseases, such as adenocarcinomas and polycystic kidney disease. In this paper we focus on cystogenesis as an in vitro model of 3D epithelial growth and describe a novel theoretical framework that discriminates between the two radically different dynamical regimes observable when the cell aggregate grows either close to, or far away from, mechanical equilibrium.

We show that cell–cell contact topology can efficiently discriminate between these two regimes. By comparing theoretical predictions for cell arrangements and the 3D reconstruction of experimental MDCK cysts, we demonstrate that cystogenesis occurs out of mechanical equilibrium.

The out-of-equilibrium condition can be rephrased in the observation that topological rearrangements in cysts involving the exchange of next-neighbors are rare events, occurring on timescales much longer than the duplication time. This leads to a peculiar topological distribution of cell-to-cell contacts (Fig. 1 D) that is the true fingerprint of a strong departure from equilibrium.

Neighbor rearrangements require that the cysts transiently assume energetically unfavorable configurations, for instance during the shrinking of a basolateral cell–cell contact area that is required in the process of neighbor exchange. Random cell motility favors the crossing of such topological energy barriers, therefore playing the role that thermal agitation has in statistical physics systems. Indeed, induced activation of ERK1/2 in MCF10A acini results in increased cell motility and neighbor exchanges (Pearson and Hunter, 2007).

Our experiments on cyst-like and tubular structures show clear hallmarks of nonequilibrium (Fig. 1) both in vivo and in vitro. These features are not shared by all epithelial tissues. For instance, *Drosophila* wing epithelium at pupal stage shows a nearly perfect hexagonal packing, corresponding to the equilibrium configuration for planar tilings (Classen et al., 2005). Similarly highly ordered, close-to-equilibrium structures are observed in the fovea, lens, and corneal tissues of mammals (Bhat, 2001; Jalbert et al., 2003; Hofer et al., 2005). The achievement of this geometrical order has been related to an increased ability in rearranging bonds between cells (Amonlirdviman et al., 2005; Classen et al., 2005) induced by various mechanisms of active cell remodeling, such as increased trafficking of cell adhesion molecules (Classen et al., 2005) and anisotropic accumulation of myosin II (Rauzi et al., 2008). Cell division slowdown may also be at play in these tissues. Overall, these contributions act to favor barrier crossing and the approach to an equilibrium configuration.

In MDCK cysts it has been reported that ROCK inhibition rescues the orientation of polarity in polarity-disrupted cells (Yu et al., 2008). Here we have shown that the inhibition of ROCK activity gives access to the equilibrium topological distribution. ROCK inhibition is known to relieve cell cortical tension (Matthews et al., 2006), thus lowering the energy barriers that separate topologically different configurations and therefore favoring the approach to equilibrium. We have found that the multilumen phenotype induced by aPKC-PSi treatment cannot be rescued with a ROCK inhibitor, as shown in Fig. S1 P, which suggests that lumen coalescence and that topological rearrangements of cell–cell contacts are unrelated processes.

Another important feature that discriminates equilibrium versus nonequilibrium phenotypes is the presence of single or multiple lumens. Our results indicate that when proliferation rates are high with respect to the characteristic time for lumen coalescence, multiple lumens may accumulate. Conversely, when proliferation rates are low, multiple lumens tend to coalesce into a single lumen. Additionally, when proper alignment of the mitotic spindle is disrupted by perturbing the normal establishment of cell polarity, aberrant multiple lumen cysts are even more frequent.

In this case as well, nonequilibrium plays a determinant role: our results show that if cystogenesis were a pure equilibrium process, mechanical relaxation would ultimately always lead to the single-lumen phenotype independently of cell polarity, as shown in Fig. 6 E.

Our model explains why even in the case of normal cells, a significant fraction of cysts does not show a single lumen: indeed out-of-equilibrium growth can give rise to aberrant phenotypes even in presence of fairly accurate positioning of the spindle.

Cell-cycle signaling pathways are often altered in cancers, and accelerated proliferation represents one of the most typical features of carcinomas. In this context, the observation that control of cell duplication time alone can effectively have an impact on the properties of cyst growth is particularly relevant. Our simulations indicate that faster cell divisions would give rise to aberrant phenotypes even in the case where the geometry of cell division is perfectly controlled. These observations lead us to speculate that cancerous architectural patterns could also emerge due only to faster cell divisions, through a polarity-independent mechanism, thus calling for more detailed experimental investigations.

The present findings raise the following question. Since equilibrium dynamics would, by themselves, ensure a proper cyst geometry without the need for tight control over the mitotic spindle orientation, why does cystogenesis not take place in the equilibrium regimen? The answer probably lies in the fact that it is difficult to reconcile the contrasting requirements of fast growth and fast relaxation. Allowing sufficient time for mechanical relaxation to equilibrium inevitably makes the overall endeavor of cyst construction much slower. To make matters worse, the number of local energy barriers to be overcome before reaching equilibrium increases exponentially with the number of cells, with a consequent increase in the relaxation time. Eventually the cyst is trapped into an aberrant configuration. Control of the orientation of cell division circumvents this topological stalemate.

This leads us to speculate that evolution has warranted fast and orderly cyst growth through an accurate control of cleavage plane orientation, thereby remedying to the frustrating constraints imposed by mechanics.

## Materials and methods

### Experimental procedures

MDCK (subcloned from ATCC-CCL34 as described previously; Lipschutz et al., 2001) were cultured at 37°C, 5% CO_{2}, in Eagle’s minimal essential medium supplemented with 10% FBS, with the addition of l-Glutamine and antibiotics. To grow cysts, culture slides (Falcon 354108; BD) were precoated with a layer of basement membrane extracts (Matrigel; BD). Cells were seeded with a density of 50–100 cells/mm^{2} and fed every other day in complete media supplemented with 2% Matrigel. aPKC activity was perturbed by the addition of 50 µM C-Zeta pseudosubstrate, myristoylated (aPKC-PSi; Biosource) in growth media (Martín-Belmonte et al., 2007). ROCK activity was perturbed by the addition of 20 µM Y-27632 (EMD Millipore) to the medium. Aphidicolin was used at 1 µM. Cells were fixed and stained according to the protocol used in Yu et al. (2005). The primary antibodies used were: mouse anti-gp135/PCX (a gift from G. Ojakian, Department of Cell Biology, SUNY Downstate Medical Center, Brooklyn, NY), rabbit anti–B-catenin (Santa Cruz Biotechnology, Inc.), and rat anti–ZO-1 (R40.76; a gift from B. Stevenson, Nationwide Children’s Hospital, Columbus, Ohio). Alexa fluorophore–conjugated secondary antibodies (Molecular Probes), Hoescht and DAPI (to counterstain nuclei), and Alexa Fluor–conjugated Phalloidin were used. Whole-mount samples were dissected from embryonic day 13.5 (E13.5) C57/Bl6 mice and fixed in 4% paraformaldehyde for 10 min. Samples were washed in phosphate-buffered saline with 0.1% Triton X-100 (PBTX), blocked (PBTX + 10% FCS), and then incubated for 48 h in primary, then secondary antibody solutions, washing with PBTX in between. Mouse anti–E-cadherin (No. 610181; BD), and rabbit anti–ZO-1 (No. 40-2300; Life Technologies) primary antibodies were used and detected with Alexa flurophore–conjugated secondary antibodies. Samples were then imaged on a laser scanning microscope (LSM510 and Leica SP2 equipped with PMT detectors, 40× oil objective lens, 1.25 aperture; Carl Zeiss). The expression vector for MLC-GFP was transiently transfected into MDCK using Lipofectamine, according to the manufacturer’s recommended procedure, and imaged after 24 h.

### Image analysis

Image segmentation of cysts with formed lumen was performed with the help of ImageJ, SciPy, Mayavi, and MATLAB software. To measure topological distributions in vivo, sections from confocal stacks were semiautomatically segmented, and neighbors were counted by means of MATLAB custom-written algorithms. Cysts were visualized by 3D rendering with a custom written code in SciPy. Luminal surfaces were measured from segmented stacks by fitting the data with triangulated surfaces and calculating their areas. In Figs. 2 D and 5 B, the distance from equilibrium was measured as *d*_{eq}/(*d*_{eq} + *d*_{neq}), where *d*_{eq} and *d*_{neq} are the Euclidean distances of the vector

### Theoretical model

We assume that the mechanical energy associated with a given cyst geometry is*S*_{cell−cell} denotes the contact area between two cells (and similarly the contact area with ECM and lumen). The surface tensions *α*, *β*, and *γ* are hierarchically ordered: *γ* > *β* > *α* > 0. This ordering implies that cell–cell adhesion is energetically more favored than cell–ECM contact, leading to the formation of a spherical cell aggregate. The tendency to minimize mechanical surface energy is hindered by physical constraints, such as the limited variability of cell volume, and biochemical ones, such as the presence of signaling domains that control the size of the apical regions (Martín-Belmonte et al., 2007; Veglio et al., 2009; Bryant et al., 2010; see also Fig. S1). These two constraints are implemented through additional terms *C*_{V}(*V*_{cell}), where *V*_{cell} is the cell volume, and *C*_{A}(*S*_{cell−lumen}). These penalization terms have a single minimum at a given target cell volume and target cell-lumen area, respectively. The total energy therefore reads

When a division occurs, a cell is split in two along a cleavage plane picked randomly in the plane parallel to the luminal surface. For cells lacking proper spindle positioning, the orientation of the cleavage plane is chosen completely (or at different extents) at random in 3D space. The model can be easily extended to the case of shape-dependent cell division (Gibson et al., 2011), without substantially altering the predictions. After each division, the energy is redefined according to the new cell configuration. The constraint contribution to the energy has a much shorter relaxation time than the surface energy part, so that during the dynamical evolution, cell volume and apical surfaces are only transiently different from their target values.

### Numerical simulations

Cystogenesis was simulated using a dynamic Potts model (Potts, 1952; Wu, 1982; Graner and Glazier, 1992; Glazier and Graner, 1993). On a 3D cubic lattice, each site is occupied by a “color” variable that takes on integer values. To each value there corresponds an individual cell. At each division a new color is added. The exchange energies on neighboring sites provide a lattice approximation of the surface energies of cell–cell, cell–lumen, and cell–ECM contacts described in Eq. 1. Target cell volume and luminal surface are introduced on the basis of experimental observation and enforced by additional terms in the energy. The target cell volume and lumen areas are set to *V*_{0} = 1,200 µm^{3} and *S*_{0} = 80 µm^{2}, according to experimental values (Fig. S1). Simulations always started from a single cell. Duplication times are chosen from a uniform distribution with mean *τ* and coefficient of variation 0.1. To simulate normal cells, the division axis is chosen at random in the plane tangent to the lumen surface, whereas in the aberrant case it is allowed to depart from the tangent plane to various degrees. We verified that our results do not depend significantly on the precise values of the surface tension parameters in Eq. 1, as long as the hierarchy γ > β > α is respected. The results are similarly independent of the precise form of the constraints in Eq. 2, as long as they have single minima and a relaxation time much shorter than the typical relaxation time of the surface energy term.

### Stochastic algorithm

The dynamic Potts model is driven by a discretization of the energy assigned by Eqs. 1 and 2. We consider a 3D cubic lattice of *N* = *n*^{3} sites, with *n* = 100–300 and periodic boundary conditions. We assign to each site *i* a variable *s _{i}* whose value indicates if the site is occupied by ECM (

*s*= 0), lumen (

_{i}*s*= −1), or one of a collection of

_{i}*n*

_{c}cells (

*i*= 1, … ,

*n*

_{c}). Every color cluster corresponds to a cell. We refer to the

*s*values as “spins” or “colors,” following the traditional nomenclature of statistical physics. Discretized surface energies are computed by the formula

_{i}*i*to its 26 cubic neighbors. The symmetric tensor

*H*the penalization term

_{m}*λ*

_{A}satisfies a single requirement: it must be large enough to provide a relaxation time for the constraint contribution much faster than the characteristic time for the relaxation of the surface energy. Target volumes are enforced by the Monte Carlo dynamics by allowing only moves that tend to decrease the distance between the actual cell volume and the target volume

*V*

_{0}, as described in Creutz (1983). We start the simulations from a single spherical cell with “color” 1. Spins relax according to a Monte Carlo dynamics with inverse temperature β (Landau and Binder, 2000). The precise value of β is not relevant as long as β is large enough to smear lattice anisotropies associated with the cubic discretization and small enough to allow the formation of stable spin clusters. At each time step, a randomly chosen spin is “flipped” to another color chosen among those taken by its neighbors, corresponding to an elementary extension or retraction of one cell to the expenses of its neighbors. This random move is accepted if the change in energy Δ

*H*is negative. It is rejected with probability 1 −

*e*

^{−βΔH}if Δ

*H*> 0. When a move is accepted, another spin is flipped, and the procedure is repeated iteratively until the division time, chosen for each cell from a uniform probability distribution with mean τ and coefficient of variation 0.1, is reached. Then, the original cell is allowed to grow to the double of the target volume

*V*

_{0}and split in two along a randomly chosen cleavage plane. Two new colors are assigned to the newly formed cells. A germ of luminal surface is positioned at the interface between the two cells by assigning −1 spins in a localized circle around the center of the cleavage plane. The orientation of the cleavage plane is always randomly chosen, but the plane is constrained to be perpendicular to the cell apical surface for wild-type cells, and is instead completely unconstrained for cells lacking any control on mitotic spindle orientation. To determine the amount of error that is necessary to destroy the monolumen architecture, the intermediate case was also considered, where the positioning of the cleavage plane perpendicular to the apical surface could be affected by an angular error of increasing amplitude.

### Solvable model for cyst radii

The configuration of a cyst of *N* cells, subject to the constraints that cells have volume *V* and apical surface σ can be calculated analytically in the limiting case of sharp constraints. Indeed, the internal and external radii (*r* and *R*, respectively) obey: (4π/3)(*R*^{3} − *r*^{3}) = *NV*, 4π*r*^{2} = *N*σ. From these two relations we obtain *R* = [(*N* σ/4π)^{3/2} + (3/4π)*NV*]^{1/3} .

### Image analysis

DAPI- and Phalloidin-stained cysts stacks were acquired by confocal microscopy, with LSM510 (Carl Zeiss) and SP2 (Leica) microscopes (Fig. S3, A and D). The boundaries of cells in each slice were drawn semiautomatically by looking at the original images. Both nuclei position and actin staining were used to this aim. Each cell was tracked across the stack by assigning it a unique color in the stack. All cells in the cysts were then uniquely identified, as shown in Fig. S3 (B and E). Black traits corresponding to cell boundaries were then fused to neighboring cells by a standard erosion algorithm. A median filter was applied to remove isolated pixels and remaining segmentation errors (Figs. S3, C and F). Slices were then mounted together in a 3D image at the correct spatial resolution (i.e., accounting for different resolutions in the x-y plane and in the z direction). The result was then visualized by a rendering routine written in Python and MayaVi (Figs. S3, G and H). To extract the outer (or inner) topology, a 3-pixel-thick shell was extracted and contacts were then computed on the shell by calculating the neighbor map. Results were then carefully corrected by eye on the full 3D structure to identify possible errors coming from the finite spatial resolution. To quantitatively characterize MLC-GFP spatial distribution, we segmented midsections of digital images of stained cysts by using DAPI, Phalloidin, and GP135 staining. The three channels were used to reconstruct a binary mask identifying the section of the cyst. Within the mask, we identified cells by segmenting nuclei and then applying the watershed transform into the mask. Two circular layers spanning a quarter of the thickness of the cyst (one starting from the basal and the other starting from the apical side) were digitally constructed. The ratio of fluorescence intensity of the MLC-GFP channel between the basal and apical quarter was then used as a quantification (Fig. S5).

### Statistical analysis of topological distributions

The statistical likelihood that a given experimental distribution originates from nonequilibrium dynamics is determined as follows. The probability of observing *n _{i}* counts in a histogram for polygons of

*i*sides, with theoretically expected frequencies μ

*, is multinomial:*

_{i}*is zero. Thus, we have to restrict the comparison to nonzero entries of the expected distribution, i.e., to polygons with 4, … , 9 sides. We refer here to one of the distributions reported in Gibson et al. (2011; Fig. 6 B, gray and blue lines). In particular, we do not take into account refinements of that model (such as, for example, those described in Gibson et al. [2006] or Sandersius et al. [2011]). The original Gibson distribution has been derived by considering random unbiased cell division, which is known to be an oversimplification. However, the effect of all the corrections that have been made afterward is to change the final distribution by a few percentage points, and does not represent, from a quantitative point of view, a major difference. In our work we do not discuss any of these factors, such as size-dependent cell division rates, or shape-dependent control of mitotic spindle. The two opposite cases treated in the present work are the equilibrium and out-of-equilibrium cases, which do show significantly different distributions. The likelihood of observing a given realization, given the theoretical distribution, and the corresponding errors can be evaluated by Bayes’ formula:*

_{i}*P*(μ|

*n*) =

*P*(

*n*|μ)

*P*

_{0}(μ)/

*P*(

*n*). Following a maximum likelihood approach, we assume a uniform prior, i.e.,

*P*

_{0}= const, and maximize

*P*(μ|

*n*) by maximizing

*P*(

*n*|μ). Corresponding values of μ* such that the likelihood function is maximized are those satisfying:

*λ*is a Lagrange multiplier used to impose ∑μ

_{i}= 1. The solution is μ* =

*n*

_{i}/

*N*. We estimate the Euclidean distance between the estimated μ* and the expected μ as

^{2}log

*P*/∂μ

*∂μ*

_{i}*= (*

_{j}*N*

^{2}/

*n*)δ

_{i}*. Here, the curvature of the Hessian represents the error. For the*

_{ij}*k*−th bin, this is

*d*must be compatible with zero with a standardized test, i.e.,

*d*± 3Δ

*d*must include 0 at the 99% confidence level. To implement these quantifications, we generated several synthetic distributions by extracting random numbers according to the probabilities given in Table S1. The corresponding distances and probability values are given in Fig. S4, C and D. These statistical methods, applied to the data extracted from Figs. 1 and 3, give the results reported in Table S2. The same methods, applied to the data extracted from Fig. 4 C, give the results reported in Fig. S4 (A and B).

### Online supplemental material

Fig. S1 provides quantitative data about measured epithelial geometry and topology. Fig. S2 shows digital 3D reconstructions of fixed and stained MDCK cysts. Fig. S3 details the image segmentation process. Fig. S4 contains data about the statistical analysis of cyst topology. Fig. S5 shows a quantification of GFP-MLC distribution in untreated and ROCK-inhibited cysts. Table S1 contains predicted nonequilibrium probabilities for the number of cell–cell contacts. Table S2 shows distances of experimental data from the nonequilibrium distribution. The following custom written routines for 3D stochastic simulations and image visualization are provided in a zip file: potts3d.f is a Fortran code for dynamical simulations of 3D Potts model; image_analysis.py is a Python code for image analysis using SciPy and MayaVi libraries; invivo_topology.py is a MATLAB code for the analysis of in vivo data; and palette is the list of colors used by the routines in “image_analysis.py.” Online supplemental material is available at http://www.jcb.org/cgi/content/full/jcb.201305044/DC1.

## Acknowledgments

This study was supported by Fondazione Piemontese per la Ricerca sul Cancro-Onlus, intramural grant 2010-2012. Additional financial support from Telethon GGP09175, Associazione Italiana per la Ricerca sul Cancro (AIRC) IG9211, FPRC/PAMI/2010, Istituto Nazionale Di Fisica Nucleare (INFN) TO61, and HuGeF is gratefully acknowledged. This work was supported by National Institutes of Health grants R01 DK074398, 5R37 AI25144, P01 AI53194, and R01 DK091530 to K.E. Mostov, and National Health and Medical Research Council (HHMRC) C.J. Martin Fellowship to A.M. Shewan. M.H. Little is an NHMRC Principal Research Fellow.

Author contributions: B. Cerruti, A. Puliafito, A. Celani, and A. Gamba developed the theoretical model in collaboration with K.E. Mostov and G. Serini. A.M. Shewan, W. Yu, F. Chianale, and A. Puliafito performed the experiments. B. Cerruti performed computer simulations and computer-assisted analysis of the data on cysts. A. Puliafito performed computer-assisted analysis of the in vivo data. L. Primo provided reagents and materials. A.N. Combes and M.H. Little provided images of mouse tissues. B. Cerruti, A. Puliafito, A. Celani, and A. Gamba wrote the manuscript. A. Celani, K.E. Mostov, and A. Gamba share authorship of this work.

## Footnotes

- Abbreviations used in this paper:
- AMIS
- apical membrane initiation site
- E
- embryonic day
- MLC
- myosin light chain
- PCX
- podocalyxin
- PSi
- pseudosubstrate inhibitor

- Submitted: 9 May 2013
- Accepted: 22 September 2013

This article is distributed under the terms of an Attribution–Noncommercial–Share Alike–No Mirror Sites license for the first six months after the publication date (see http://www.rupress.org/terms). After six months it is available under a Creative Commons License (Attribution–Noncommercial–Share Alike 3.0 Unported license, as described at http://creativecommons.org/licenses/by-nc-sa/3.0/).