Exploring the effects of mechanical feedback on epithelial topology

During development, growth must be tightly regulated in order to ensure the formation of organs that have a specific size and to maintain their structural integrity. Epithelia, for example, generally need to consist of tightly packed cells and be capable of resisting applied mechanical stress. The wing imaginal disc of Drosophila serves as a model system to study growth and epithelial packing. This disc is an epithelial sac-like structure that grows from about 30 to roughly 50,000 cells during the larval stages and develops into the adult wing during metamorphosis (Milan et al., 1996). During the larval stages, cells have different numbers of neighbors, which we here refer to as having different polygon numbers. The disc shows a reproducible distribution of polygon numbers, which is conserved among Metazoa (Gibson et al., 2006). This topology most probably emerges from characteristics of single cells, such as their mechanical properties, interactions with neighboring cells and proliferation rate (Farhadifar et al., 2007; Nagpal et al., 2008). Insight into this process could reveal aspects of how epithelial integrity is obtained in a growing tissue or about certain features of growth regulation itself.

Studies in the wing imaginal disc have revealed several signaling pathways that are involved in growth regulation, including insulin, Wnt and Hippo signaling (reviewed in Johnston and Gallant, 2002). We have hypothesized that growth is also regulated by mechanical stress and have included this idea into a model for size regulation (Aegerter-Wilmsen et al., 2007). Hufnagel et al. have proposed an alternative model for size control in which mechanical forces also play a central role (Hufnagel et al., 2007). Even though both models differ fundamentally in how stresses build up, they both postulate that compression increases in the center of the disc and that this compression terminates growth. In contrast to mammals, a growth regulatory effect of mechanical forces is purely hypothetical in Drosophila. For the wing disc, such an effect is difficult to assess experimentally because the disc cannot readily be accessed mechanically in vivo and there is no satisfactory in vitro culture system available.

The wing imaginal disc comprises two cell layers, the columnar layer and the peripodial membrane. In the more apical part of the columnar layer, cell shapes approximate polygons. The more basal parts of the columnar layer do not form simple columns of polygonal shape but have a more irregular shape (see Fig. S1 in the supplementary material). Most of the cell body is concentrated around the nucleus, and cell bodies of neighboring cells are packed on top of each other (Hipfner et al., 2004). Therefore, it seems unlikely that the apical cell area is strictly determined by the volume of the cell. Instead, it seems as if the apical junctions are under tension and that cell area at this level is shaped by forces that result from interactions with neighboring cells in combination with forces that are generated by the cytoskeleton. Indeed, cutting a single-cell boundary with a laser beam resulted in cell shape changes in theneighboring field of cells at the level of apical junctions up to about 8-times the average cell side length away from the cut. This tension at cell boundaries was shown to depend on actin-myosin contractility (Farhadifar et al., 2007). Thus, apical cell shapes are at least partly determined by mechanical forces that are exerted by neighboring cells.

In order to gain insight into mechanisms underlying tissue topology, Gibson et al. have used mathematical modeling to study the theoretical effect of cell division (Gibson et al., 2006). As can be seen in Fig. 1A, division of a cell can alter the polygon number of daughter cells as well as neighboring cells. Gibson and colleagues argued that, if all cells have asynchronous, but roughly uniform cell cycle times, division leads to an average polygon number of six. Furthermore, they concluded that such uniform cell division can theoretically lead to the convergence of epithelial topology to a fixed equilibrium distribution of cellular polygons and that this distribution is in agreement with the experimentally observed distribution. Moreover, their model predicts that each cell in the population must gain one side on average per cell cycle and they indeed found that the average mitotic cell in the wing disc has seven sides instead of six (Gibson et al., 2006). The authors have not directly compared the polygon distribution of mitotic cells with a predicted distribution.

Farhadifar et al. also studied mechanisms underlying epithelial packing in the wing imaginal disc (Farhadifar et al., 2007). They used a two-dimensional vertex model that takes into account physical properties of the cells. In addition to division, the model also allowed for cell rearrangements (also referred to as T1 transitions), which affect the polygon numbers of the cells involved (Fig. 1B). They found a region of the parameter space where cells showed the experimentally observed polygon distribution. Part of this region also reproduced quantitative responses of the wing disc to perturbations by laser ablation, as well as the positive correlation between average cell area and polygon number that they observed experimentally. The simulation reproduced this correlation in a cell-volume-independent way (Farhadifar et al., 2007). It has not been tested whether this model can also reproduce the polygon distribution of mitotic cells.

Thus far, it has not been studied how growth regulation by mechanical stress could have an effect on topology. If such a regulation exists, it would be expected to affect topology based on the following. The polygon distribution can be altered if division frequencies are shifted in a polygon-number-dependent way. Cells with different polygon numbers could sense different mechanical forces from their neighbors (Fig. 1C). This is because a cell with a low polygon number is surrounded by neighboring cells that tend to be elongated along its circumference as a result of geometric constraints, whereas a cell with a high polygon number is surrounded by neighboring cells that are elongated perpendicular to its circumference (Fig. 1C). Because the contractility of circumferential actin bundles at the apical level of an epithelium stimulates cells to round up (Owaribe et al., 1981) and thus oppose elongation, a cell with a low polygon number tends to be more compressed by its neighbors, whereas a cell with a high polygon number tends to be more stretched. This is supported by the experimental finding that apical cell area positively correlates with polygon number (Farhadifar et al., 2007). Because of these mechanical effects of neighboring cells, growth regulation by mechanical stress at the apical junctions of cells could lead to polygon-number-dependent growth rates (where a growth rate is here referred to as the rate at which cell volume increases). If it is assumed that cells divide once they reach a certain threshold volume, these polygon-dependent growth rates would lead to polygon-dependent division rates. Note that, even though we use the definition given above for growth rate, including its relationship with volume in this paper, a similar argument can be given if growth is seen as cell cycle progression. The precise relationships between volume increase and cell cycle progression in wing discs is not yet clear (Neufeld et al., 1998; Coelho and Leevers, 2000). Thus, in summary, growth regulation by mechanical stress could affect the polygon distribution, as cells with increasing polygon numbers could sense increasing mechanical stretching exerted by their neighbors and this would lead to increasing growth rates for cells with increasing polygon numbers. These polygon-dependent growth rates would in turn cause polygon-dependent shifts in division rate, which could then affect the polygon distribution.

Here, we extend the work of Gibson et al. and Farhadifar et al., by investigating whether division alone, and division in combination with rearrangements, respectively, can also theoretically account for the polygon distribution of mitotic cells. In addition, we study the theoretical effect of mechanical-force-dependent growth rates on topology. We formulated four scenarios (Table 1) and compared simulation outcomes of these scenarios with experimental data we generated. In the first scenario, the polygon distribution emerges from roughly uniform cell divisions only. In the second scenario, cells are also allowed to rearrange. In the third and fourth scenarios, cell cycle length depends on deviations of average cell areas, which we use as a measure of mechanical stress. The fourth scenario also allows for rearrangements. For all scenarios, we use a modeling method that is very similar to the one employed by Farhadifar et al., and thus includes mechanical properties of cells (Farhadifar et al., 2007). The polygon distribution of mitotic cells is sensitive to whether neighboring cells are in mitosis at the same time. As such clustering of mitotic cells has previously been found in the wing disc (Milan et al., 1996), we also compared the mitotic clustering that was produced by the simulations with the clustering that we measured experimentally.

Surprisingly, we found considerably less clustering of mitotic cells experimentally than what has been reported previously. All scenarios reproduce our data on clustering. When comparing polygon distributions of the entire population and of the mitotic cells, we found that only scenarios 3 and 4, which include area-dependent growth, are in agreement with the data, thus arguing in favor of a role for mechanical stress in growth regulation. These scenarios also show decreased apoptosis and cell rearrangements, respectively, suggesting that, if growth regulation by mechanical stress would exist, it could be favorable for maintaining a structurally integral tissue during growth, while largely preventing the extrusion of healthy cells at the same time.

For immunocytochemistry and Phalloidin staining, standard protocols were used. Images for polygon count analysis comprised immunostained, y,w genotype, third instar wing discs that were collected with a Leica TCS SP5 confocal microscope and a 63× 1.4 NA objective. Alexa Fluor 647 Phalloidin (Invitrogen) was used at 1:40, rabbit anti-pH-H3 (Upstate) was used at 1:400, and mouse anti-Dlg (Developmental Studies Hybridoma Bank) was used at 1:300. The secondary antibodies were Alexa Fluor 488 goat anti-rabbit and Alexa Fluor 594 goat anti-mouse (both 1:500, Molecular Probes). Images for apical surface area and volume measurement were collected using a Zeiss LSM710 confocal microscope with a 63× 1.4 NA objective. Clones were induced by a 10-minute, 37°C heatshock in mid-third instar larvae carrying an Act-FRT-CD2-FRT-GAL4 flip-out cassette (Pignoni and Zipursky, 1997) and a UAS-CD8-GFP transgene (Bloomington Stock Center line 5137). To ensure small clone size, discs were dissected and stained 4 hours after induction. These discs were stained with rabbit anti-GFP (1:400) to improve GFP brightness and mouse anti-Dlg (Developmental Studies Hybridoma Bank, 1:300). Secondary antibodies used were Alexa Fluor 488 goat anti-rabbit and, Alexa Fluor 594 goat anti-mouse (both at 1:500, Molecular Probes). These image stacks were collected using Nyquist-optimized sampling rates, and deconvolved using the Zeiss ZEN software package. Measurement of volume and apical surface area was carried out using the ‘contour surface’ and ‘oblique slice’ functions of the Bitplane Imaris software package, respectively, where the ‘oblique slice’ function was used to compensate for folds in the tissue.

In order to obtain the polygon distributions of the entire cell population, the most apical Dlg-stained cell shapes were analyzed using a packing analyzer developed by B. Aigouy (Farhadifar et al., 2007). The error was considered to comprise three parts. First, there was a variation from disc to disc, which was obtained by calculating the standard deviation from measurements of seven images (2928 cells in total). Furthermore, when the length of a cell side is only 1 pixel, the program cannot determine which cell it belongs to and does not assign it to any cell. The frequency by which this occurs was calculated using the fact that one would theoretically expect six sides per cell on average (Gibson et al., 2006) and was found to be about 10% of the counts per polygon class. Lastly, sometimes it is very difficult to assign sides to cells even if they are a few pixels wide. This is because connections may also seem different in different confocal planes (Fig. 5, not only the connection between the mitotic cells, but also a connection at the top right corner of the mitotic cells). We estimate this to constitute an absolute error of 1%.

In order to obtain the distribution of mitotic clusters, we manually analyzed 28 images (304 clusters). Phospho-histone-H3-stained cells were only considered if we could unambiguously determine which apical cell shape belonged to the mitotic cell and when all the neighboring cells were clearly apparent. The low mitotic frequency, in addition to a variation in mitotic frequency among different discs, did not allow for a direct comparison of counts in different images in order to obtain the variation. Instead, clustering was determined by counting the number of neighboring mitotic cells at the apical side. We estimated the error by taking the theoretical error due to counting (square root of the number of occurrences) and adding an estimated absolute disc and interpretation error of 0.1%.

The polygon distribution of mitotic cells was obtained manually from the same 28 images (334 mitotic cells) and was determined for the apical sides of the Phopho-histone-H3-stained cells. Phalloidin and Dlg antibody double-staining allowed for improved polygon number analysis, as these stainings do not completely overlap: Phalloidin stains adherens junctions most intensely, whereas Dlg is localized to the septate junctions, which are positioned basally from the adherens junctions (Condic et al., 1991; Woods et al., 1997). Here we also used the theoretical error and estimated the additional absolute disc and interpretation error to be 1%.

All errors were added quadratically. Distributions were compared using a χ²-test.

For modeling growth, a similar method was used as was employed previously (Farhadifar et al., 2007). In short, apical cell shapes are considered as polygons that are defined by their vertex positions, where three cell sides (edges) meet. Periodic boundary conditions are employed. Polygon shapes are obtained by minimizing the following energy function using a conjugate gradient method (Farhadifar et al., 2007):

The meanings of the terms and parameters are as previously described (Farhadifar et al., 2007). The parameter values correspond to case I in Farhadifar (Farhadifar et al., 2007), where all cells are identical and Γ′Γ/K A⁽⁰⁾=0.12 (−) and Λ′=Λ/K A^(0)3/2=0.04 (−).

Growth is implemented by letting cell divisions occur along a plane with random orientation. In order to allow the tissue to increase in size, the x- and y-axis are rescaled by an energy relaxation step. Dividing cells are selected based on a volume function. For scenarios 1 and 2 the following equation was used for the increase of volume for cell α per time step (ΔV): ΔV_α=V_d·0.018·r_α, where V_d is the threshold volume at which a cell divides and r_α is a cell-specific growth rate. This rate is a random dimensionless number between 0.25 and 1.75 that is assigned to daughter cells upon division. Without loss of generality, V_d was set to 1.

For scenarios 3 and 4, the following equation was used: ΔV_α=V_d·[0.015 + μ·(A_α-<A>)/A⁽⁰⁾], ΔV_α≥0, where A_α denotes the area of a cell, <A> the mean area of all simulated cells, A⁽⁰⁾ the target area, and μ is a dimensionless number characterizing the strength of cell area dependent growth regulation.

For scenarios 1 and 3, cells were removed and counted as apoptotic when they had an area smaller than 0.1× the mean cell area. For scenarios 2 and 4, rearrangements took place when the distance between two vertices was smaller than 0.2× the mean edge length (Fig. 1). When rearrangements resulted in the formation of two-sided cells, these were removed and counted as apoptosis.

In the simulations it is assumed that the preferred (apical) area of a cell (A_α⁽⁰⁾) is independent of its volume. As a control, simulations were run where the preferred cell area depends on cell volume via: A_α⁽⁰⁾/<A_α⁽⁰⁾>=1 + ν·(V_α-<V_α>)/<V_α>, where the average preferred area, <A⁽⁰⁾>, was equal to A⁽⁰⁾ in the other simulations. Increasing ν leads to a stronger correlation between polygon number and cell area. In order to compensate for this, K_α in the energy equation was increased by a factor of 1.3 for ν=0.7 (−) and by a factor of 1.5 for ν=1 (−).

All simulations started with 36 cells and were run until 10,000 cell divisions had occurred and polygon distributions were stable (see Fig. S2 in the supplementary material). For the implementation of the simulations, MATLAB (The MathWorks) was used.

Because we aimed to study whether the different scenarios can reproduce the polygon distribution of the entire cell population, as well as that of the mitotic cells, we first measured these distributions experimentally. Anti-disc large (Dlg) antibodies and Phalloidin, which stain septate and adherens junctions, respectively (Condic et al., 1991; Woods et al., 1997), served to define apical cell shapes. The most apical region of Dlg stainings have been shown to overlap with markers of adherens junctions (Woods et al., 1997), and we used the more apical confocal sections. Mitotic cells were identified using Phospho-histone H3 antibodies. The polygon distribution we measured for the entire cell population (Fig. 2) is in good agreement with data from Gibson et al. (Gibson et al., 2006), as well as from Farhadifar et al. when they performed the analysis in an identical manner (see Materials and methods) (Farhadifar et al., 2007). For the polygon distribution of the mitotic cells,we also measured a distribution that is similar to the one found by Gibson et al. (Gibson et al., 2006) (Fig. 2).

The distribution of the polygons of mitotic cells is influenced by clustering of mitotic cells. For example, because a cell can gain additional sides when neighboring cells divide, the polygon number of a cell at the beginning of mitosis is sensitive to whether the neighboring cells divided earlier, later or simultaneously. Therefore, in order to compare simulated and observed mitotic polygon distributions, clustering must be measured experimentally and then be reproduced by the simulations. As we use Phospho-histone H3 staining in order to identify mitotic cells for the polygon measurements, we use the same staining to analyze clustering. We found that 2.2±0.3% of the cells are mitotic, which is in good agreement with earlier data (Milan et al., 1996). When considering all clusters (including ‘clusters’ of single cells), we found that 91% contain one, 8% contain two and 1% contain three or more cells. This amounts to considerably less clustering than has been reported previously, where about 35% contained one, 25% contained two, 20% contained three and the rest contained four or more cells (Milan et al., 1996). Possible explanations for this difference are treated in the discussion.

In order to model the apical cell shapes in a growing wing disc, we used a simulation method where cells are described as polygons (Farhadifar et al., 2007; Hufnagel et al., 2007; Rauzi et al., 2008) (see Movies 1-4 in the supplementary material). Cell contacts change upon cell division, rearrangements and apoptosis. Cell shapes are obtained by minimizing an energy function that was developed by Farhadifar et al. (Farhadifar et al., 2007) (see Materials and methods). The function contains terms for area elasticity, line tension and cell perimeter contractility respectively. Even if the function would not precisely describe the cellular mechanics, it can still quantitatively reproduce different aspects of the biological behavior, as it has been fitted to different experimental data on the wing disc (Farhadifar et al., 2007).

In contrast to Farhadifar et al., dividing cells are not chosen at random. Instead, each cell is assigned a volume, which increases over time, dependent on the growth rate (see Materials and methods). When this volume reaches a certain threshold, the cell divides and the daughter cells each obtain half of the volume of the mother cell. The model represents cell shapes as they are found in the more apical parts of the cells, where cell shapes resemble polygons. In order to get an indication of whether cell volume plays an important role in determining cell area at this level, we measured the correlation between cell volume and cell area (n=60; see Fig. S1 in the supplementary material). As we did not find a significant correlation, we assume that the volume does not directly determine cell area. A growth-regulating effect by mechanical stress was included by making the growth rate dependent on the deviation of the area of a cell from the mean area of all cells, such that the growth rate increased with increasing area (see Materials and methods).

In order to mimic the biological situation as well as possible, we aimed to choose the growth step size in such a way that the same fraction of cells is in mitosis as is in vivo (2.2±0.3%). When all cells are assigned with exactly the same growth rate and consequently have the same division rate, two daughter cells will divide again at the same time, and the resulting four cells will again divide at the same time. Therefore, large groups of cells will divide simultaneously after several rounds of division. This is biologically unrealistic. Furthermore, it has been shown that the lengths of phases in the cell cycles vary largely among different cells. For example, some cells progress from S phase into mitosis in only 30 minutes, whereas others require 3.5 hours (Milan et al., 1996). In the scenarios where mechanical stress is involved in growth regulation, growth rates will naturally vary, as different cells will be exposed to different forces from their neighbors. For the other scenarios, a growth rate was randomly assigned to each daughter cell upon division, and this rate varied over a factor of 7.

In order to compare the simulation outcomes of scenarios 1 and 2 with the experimental data, no parameter fitting was required so that the experimental data could directly be compared with the simulation outcomes. This is because we used the energy equation, including its parameter values, that has been published by Farhadifar et al. (Farhadifar et al., 2007) and the simulations do not include any additional free parameters. Scenario 1, where the distribution results from division only, clearly did not reproduce the observed polygon distribution of the whole cell population (Fig. 2A, χ²=80, P<0.0001), nor that of the mitotic cells (Fig. 2C, χ²=366, P<0.0001). As this scenario was able to reproduce the observed mitotic clustering, the latter discrepancy cannot be explained by alterations in this clustering. By contrast, scenario 2, where cell rearrangements are allowed, could reproduce the polygon distribution of the entire cell population (Fig. 2D, χ²=9, P=0.2). However, it could not reproduce that of the mitotic cells either (Fig. 2F, χ²=142, P<0.0001), even though it was in agreement with the data on mitotic clustering (Fig. 2E, χ²=0.51, P=0.8).

In contrast to scenarios 1 and 2, scenarios 3 and 4 could not directly be compared with the data, as growth regulation by mechanical force is hypothetical and the values of the parameters involved are therefore not known. We assumed a linear relationship between cell area and growth rate (see Materials and methods) and varied the gradient μ systematically from 0.01 to 0.12. The obtained simulation outcomes are shown in Fig. 3. Both scenarios 3 and 4 can reproduce the polygon distribution of the entire cell population for the whole range of values of μ tested (Fig. 3A,A′,D,D′). They can no longer reproduce this distribution if μ is higher than 0.45 and 0.25 respectively. The fraction of hexagons becomes too large then and the fraction of cells with lower polygon numbers too small. Scenarios 3 and 4 are in agreement with the data on mitotic clustering for μ=0.02 to 0.12 (Fig. 3B,B′,E,E′). When μ=0.01 was used, more clustering was obtained. Finally, the scenarios can reproduce the polygon distribution of mitotic cells for μ=0.03 to 0.08, and for μ=0.01 to 0.1, respectively (Fig. 3C,C′,F,F′).

Thus, only scenarios 3 and 4 can reproduce the observed polygon distribution of mitotic cells. When considering all of the data, scenario 3 is in agreement with them when μ is between 0.03 and 0.08. For scenario 4, this is the case when μ is between 0.02 and 0.1. This indicates that μ does not have to be regulated very tightly for both scenarios to be in agreement with the data and that the system is thus relatively robust against changes in the value for μ. Scenario 3 as well as scenario 4 show the lowest average Chi square when μ=0.04. For comparison, the results for this best-fit have been used for Fig. 2.

In order to assess whether our conclusions are dependent on the precise details of the scenarios, we also tested a number of modifications. The growth functions were changed in several ways. For scenarios 1 and 2, the random variation in growth rates was changed from a factor of 7 to factors of 1.5 and 15, respectively. For scenarios 3 and 4, cell perimeter was used instead of cell area. Lastly, for each scenario, growth was made exponential instead of linear. None of these changes affected our conclusions (data not shown).

In the scenarios, we used the energy function that was developed by Farhadifar and colleagues (Farhadifar et al., 2007). The parameter values corresponded to the ones these authors used for their case I. These values lie in the center of a region of the parameter space where the energy function can reproduce their data. When testing parameter value combinations at the border of this region [(Γ′,Λ′) corresponds to (0.10,0.04), (0.16,0.04), (0.16,0.03), and (0.20,0.03), respectively], we reached the same main conclusions, ie. only scenarios 3 and 4 can reproduce all of the data (data not shown).

As we did not find a significant correlation between cell volume and apical cell area, we assumed that cell volume does not affect cell area and therefore no volume term was included in the energy function. However, there still exists a possibility that we missed a (weak) correlation between cell volume and cell area because of the large variation in cell volumes. Therefore, we also performed a control, where cell volume affects the preferred area of a cell. The variable ν (see Materials and methods) was varied from 0 (no effect of volume on preferred area) to 1 (the preferred area depends linearly on the volume as columns of polygonal shape). As a measure of mechanical stress, the deviation of the cell area from the average area of cells with the same volume was taken. This average area was determined based on a linear fit between the volumes and areas of all cells in the simulation. Scenario 1 did not reproduce the observed polygon distributions for any value of ν. Scenario 2 reproduced the polygon distribution of the entire cell population for values of ν up to 0.2 and did not reproduce that of the mitotic cells for any value of ν. Scenarios 3 and 4 reproduced all of the data for values of ν up to 0.7 and 1.0, respectively (data not shown). When ν had a value of 0.7, a strong correlation between cell volume and apical area could be observed in the simulations (ρ=0.87; see Fig. S3 in the supplementary material). Based on our experimental data (see Fig. S1 in the supplementary material), such a strong effect of volume on apical area is highly unlikely, indicating that scenarios 3 and 4 can reproduce the experimental data in the biologically likely range of values for the effect of cell volume on apical area.

When cells with large cell areas divide with a higher frequency than those with small cell areas, this could lead to mechanical stress release and thus might have an effect on the structural integrity of a growing tissue. In order to identify such possible effects, we compared the simulation results of the four scenarios. When visually checking the cell shapes that were obtained by the simulations, scenario 1 appeared to yield much less regular shapes than the other scenarios and scenario 4 seemed to yield slightly more regular shapes than scenarios 2 and 3 (Fig. 4; see Movies 1-4 in the supplementary material). Indeed, scenarios 2, 3 and 4 show lower variations in cell side lengths and cell areas (Table 2). Furthermore, the number of short edges (less than 0.3× the average length) was highest in scenario 1 and lowest in scenario 4. Rearrangements and area-dependent growth can thus lead to a more regular tissue structure. This could be favorable for maintaining tissue integrity, as a highly irregular structure may be associated with the presence of local stresses that could increase the risk of buckling.

During the simulations, cell shapes are found by minimizing an energy function. In this process, some cells lose their apical area. They are subsequently removed, which we refer to as apoptosis. Furthermore, some vertices (points where three cells meet) come very close to each other. If the scenario allows for rearrangements, cell contacts change (Fig. 1). As both apoptosis and rearrangements occur as a result of minimizing the system's energy, they represent a form of stress release. When comparing rearrangement and apoptosis rates among the different scenarios, we found striking differences (Table 2). Area-dependent growth regulation decreased the rearrangement frequency by approximately a factor of 3 (compare scenarios 2 and 4) and the apoptosis frequency by approximately a factor 70 (compare scenarios 1 and 3). Allowing for rearrangements also decreased the apoptosis frequency, but to a lesser extent (a factor of 18; compare scenarios 1 and 2). When both mechanical feedback was included and rearrangements were allowed (scenario 4), apoptosis disappeared almost entirely.

Here we compared four different scenarios, which vary with respect to the growth regulation mechanism and to the occurrence of cell rearrangements. We show that only the scenarios that include growth regulation by mechanical stress are in agreement with the experimental data. These scenarios contain specific equations and parameter values. In order to check whether our results are very sensitive to specific details of the scenarios, we also modified different aspects of the scenarios. This included the precise growth equation, the parameter values applied in the energy function, and the relationship between cell volume and cell area. For these modified scenarios, we reached the same main conclusions, showing that they are not very sensitive to the details of the scenarios. This does not exclude the possibility that scenarios could still be formulated that have roughly uniform growth and that reproduce the data, nor does it provide a direct proof for the involvement of mechanical stress in growth regulation. However, as scenarios with mechanical-stress-dependent growth are currently the only ones that can account for the data and do this rather robustly, we conclude that together our studies argue in favor of such a role for mechanical forces.

Gibson et al. have previously concluded that uniform cell division can yield the polygon distribution of the entire cell population in the absence of rearrangements (Gibson et al., 2006). In order to show this, they used a global method that does not take into account single cells. Their method implicitly assumes that each cell gains one additional vertex from neighboring dividing cells, independent of its polygon class. However, a cell can gain an additional side from each neighboring cell with a certain probability. Therefore, the more neighbors a cell has, the higher the average number of additional sides it will gain from dividing neighbors (see Fig. S4 in the supplementary material). In a very recent paper, the same group has found that a very different polygon distribution is obtained when a modeling method is used that takes single cells into account (Patel et al., 2009), indicating that this neighbor effect cannot be ignored. Subsequently, the authors have developed an alternative hypothetical algorithm for the positioning of division planes (Patel et al., 2009). When we exactly reproduce their simulation method, we find that the polygon distribution of the whole cell population is indeed in agreement with the data, but the distribution of the mitotic cells is not. When implementing the same division plane selection algorithm in scenario I, we find that it cannot reproduce the polygon distribution of the mitotic cells either (see Fig. S4 in the supplementary material). Thus, as neither scenario I, nor Gibson et al.'s original biological model, nor their recent modified version, can reproduce all experimental data, we consider it unlikely that it will be possible to find any model that is in agreement with the data and where the polygon distribution is solely determined by roughly uniform division rates.

Farhadifar et al. have previously developed a model that can account for the observed packing geometry (Farhadifar et al., 2007). This model has not been compared with topology data on the mitotic cells. Scenario 2, which is similar to this model, can also reproduce the data on a whole cell population, but it cannot reproduce those on the mitotic cells. As scenario 2 implicitly contains assumptions, such as how dividing cells are selected, we cannot exclude the possibility that there would be similar scenarios that would be able to account for these data on mitotic cells. However, when altering the selection algorithm, we did not find simulations that could reproduce the polygon distribution of mitotic cells.

As clustering of mitotic cells influences the polygon distribution of these cells, we measured this clustering experimentally. Surprisingly, we found much less clustering than was reported previously by Milan et al. (Milan et al., 1996). This may partly be owing to differences in experimental methods, as we used Phospho-histone H3 staining instead of Hoechst 33258 and Propidium iodide in order to identify mitotic cells. Differences in clustering are, however, not owing to the percentage of cells considered, as both methods identified about 2% of the cells as being mitotic. The difference may also be partly owing to a difference in definition used. For a good comparison with the modeling results, we only considered clustering of mitotic cells at the level where cell shapes resemble polygons. Milan et al. did not use such a strict definition. This might lead to differences as we found, for example, that mitotic cells sometimes appeared to be direct neighbors more basally, whereas this was not the case more apically (Fig. 5). For leg discs it has also been reported that cells are not always surrounded by the same neighbors apically and basolaterally (Condic et al., 1991). We do not exclude the possibility of more global clustering where mitotic cells are concentrated in small spots covering the disc.

The mitotic clustering we observed was reproduced by each scenario. As polygonal packing is not completely regular, and because the polygon distribution of mitotic cells is shifted compared with that of the entire cell population, it is difficult to determine the clustering that would occur based on chance. However, when 2% of the cells are randomly selected in such a way that the selected cells have an average polygon number of 7, it can be estimated that roughly 7% of the cells would be part of a cluster of two cells. This implies that about 96% of the clusters would comprise one cell and 4% would comprise two. When modeling, the percentage of mitotic pairs is about twice as high. This could be explained by taking into account that, upon division, daughter cells of one division are nearly always neighbors and these daughter cells have a higher chance to be in mitosis simultaneously than other cells as they have the same age. Because modeling each of the scenarios reproduces the observed apical clustering, there is no need for additional specific clustering mechanisms in order to explain this mitotic clustering.

As apical area increases with increasing polygon number, scenarios 3 and 4 predict a polygon-number-dependent variation in cell cycle lengths. Cell cycle times vary over more than a factor of 4 in these scenarios (see Fig. S5 in the supplementary material). It has indeed been shown that there is considerable variation in cell cycle time (Milan et al., 1996; Dubatolova and Omelyanchuk, 2004). It would be interesting to measure cell cycle lengths more directly in the future and compare them with the modeling outcomes.

In mammalian cells, growth regulation by mechanical stress often seems to have a biological function. For example, growth regulation by shearing forces in veins allows the veins to adjust their proportions to the actual blood flow, and regulation of growth and differentiation in bones by mechanical compression yields bones that are mechanically robust and light at the same time (Murray and Oster, 1984; Hudlicka, 1994; Resnick et al., 2003; Nowlan et al., 2007; Judex et al., 2009). Our modeling results suggest that regulation by mechanical stresses could also have a biological function in uniformly growing tissues. Counterintuitively, uniform growth may lead to the build-up of local stresses. These could, for example, be released by potentially harmful buckling of the tissue, apoptosis of healthy cells or rearrangements. If rearrangements would readily occur, this could decrease the resistance of a tissue against externally applied mechanical forces. Including growth regulation by mechanical stress decreased the frequency of apoptosis and rearrangements in the simulations and led to a more regular tissue structure. Therefore, such regulation could play a role in releasing local mechanical stresses, thus allowing for a strong, non-buckling tissue, in which little apoptosis occurs.

Thus, growth regulation by mechanical stress may help to ensure the structural integrity of a uniformly growing two-dimensional homogeneous tissue. However, the wing imaginal disc is not completely flat and contains regions with different protein activity profiles, which might correspond to different mechanical properties. Furthermore, we found that cells have an irregular cell shape, with volumes that vary by almost a factor of 4 within a single disc (see Fig. S1 in the supplementary material). These 3D shapes must change dynamically over time, not only because of growth, but also during mitosis or apoptosis, when cell bodies are found apically or cells are extruded basally, respectively. The wing disc thus needs to maintain its structural integrity under more difficult circumstances than have been treated here, and it will be an interesting future challenge to combine experimentation and 3D modeling to elucidate the underlying mechanisms.

We thank Gerald Schwank and Maria Willecke for comments on the manuscript, Sabine Schilling and Simon Restrepo for assistance and useful discussions and Christian Lehner for advice. We also thank the Developmental Studies Hybridoma Bank at the University of Iowa for supplying the Dlg antibody (originally contributed by C. Goodman). Funding for this work was partially provided by the Roche Research Foundation, the SystemsX.ch initiative within the framework of the WingX Project and the Forschungskredit of the UZH.