Abstract
Robustness of biological systems is crucial for their survival, however, for many systems its origin is an open question. Here, we analyze one subcellular level system, the microtubule cytoskeleton. Microtubules selforganize into a network, along which cellular components are delivered to their biologically relevant locations. While the dynamics of individual microtubules is sensitive to the organism’s environment and genetics, a similar sensitivity of the overall network would result in pathologies. Our largescale stochastic simulations show that the selforganization of microtubule networks is robust in a wide parameter range in individual cells. We confirm this robustness in vivo on the tissuescale using genetic manipulations of Drosophila epithelial cells. Finally, our minimal mathematical model shows that the origin of robustness is the separation of timescales in microtubule dynamics rates. Altogether, we demonstrate that the tissuescale selforganization of a microtubule network depends only on cell geometry and the distribution of the microtubule minusends.
Introduction
The correct positioning of intracellular components such as proteins and organelles is critical for correct cellular function (St Johnston and Ahringer, 2010; Ryder and Lerit, 2018). These components are transported to their biologically relevant locations by motor proteins moving along the cytoskeleton (Gagnon and Mowry, 2011; Kapitein and Hoogenraad, 2011; Franker and Hoogenraad, 2013), or through active diffusion often dependent on these motors (Drechsler et al., 2017; Colin et al., 2020). Therefore, the direction of cytoskeleton filaments guides the direction and efficiency of intracellular transport. One common type of cytoskeleton used for transport is microtubules (Franker and Hoogenraad, 2013). These are highly dynamic unstable polymers that switch between periods of growth and shrinkage. During the growth phase, GTPtubulin dimers are added to the microtubule plusend forming a GTPcap. GTPtubulin stochastically hydrolyzes into GDPtubulin, and the loss of the GTPcap results in a catastrophe: a switch to depolymerization (BowneAnderson et al., 2013). Microtubule dynamical properties are influenced by a wide range of both internal and environmental factors. For example, the dynamics and number of individual microtubules in a cell depend on the expression of particular plus and minusend binding proteins (Akhmanova and Steinmetz, 2015), the interaction between microtubules is affected by the presence of crosslinking and motor proteins (Kapitein and Hoogenraad, 2015), and the stability of the microtubule network is affected by external factors, e.g. temperature, which changes microtubule rigidity (Kawaguchi and Yamaguchi, 2010).
Animal cells take multiple shapes and forms depending on their function, ranging from neuronal cells with meterlong projections, to epithelial cells, to migrating ameboidal leukocytes. Therefore, it is not surprising that microtubule systems similarly acquire a multitude of organizations. In undifferentiated cells, microtubules form a radial array with minusends at a microtubule organizing center at a centrosome. Upon cell differentiation, microtubules are reorganized into noncentrosomal arrays of varying geometry, ranging from unidirectional bundles (e.g. axons), to bidirectional microtubule systems (e.g. subapical microtubules or microtubules in dendrites; Muroyama and Lechler, 2017). Microtubule selforganization in unidirectional or radial microtubule systems has been extensively studied (e.g. Surrey et al., 2001; Nédélec et al., 2003; Dehmelt, 2014; Kapitein and Hoogenraad, 2015). Directionality and alignment of these networks depend on several interconnected factors. These include the localization of microtubule minusends, which could be concentrated at a single location (microtubule organizing center), distributed uniformly at the cell surface, or targeted to specific locations, for example, to the sites of cellcell contacts (Muroyama and Lechler, 2017). The other two factors affecting microtubule network organization are the geometrical constraints of a cell, for example, in a long and thin axon or dendrite, microtubules can only grow in specific directions; and the presence of crosslinking and motor proteins, which promote assembly and orient microtubule bundles (Zemel and Mogilner, 2008; Zemel and Mogilner, 2009). Bidirectional microtubule networks have an additional level of complexity, as several other factors independently contribute to their organization: the dynamics of individual microtubules, their interactions, and the often dispersed distribution of minusends.
In this paper, we explore the selforganization of bidirectional microtubule networks, which are particularly common in differentiated epithelial cells – one of the four fundamental tissue types found in all animals (Gilbert et al., 1991; Bulgakova et al., 2013; Muroyama and Lechler, 2017; Tateishi et al., 2017). Here, our focus is the mixedorientation microtubules just under the apical surface of epithelial cells that are seeded from sites of cellcell adhesion at the cell periphery (Toya and Takeichi, 2016), and not unidirectional apicalbasal microtubules. The advantage of this system is that it is physically constrained in space, which allows us to model it as quasi2d (Gomez et al., 2016). We define microtubule selforganization as the degree of alignment of individual microtubule filaments with each other, which we quantify using a lengthweighted microtubule angle distribution. Hence, more peaked (or flat) microtubule angle distributions correspond to more aligned (or disordered) microtubule networks. In this work, the microtubule angle distribution describes longterm steady state of the system. In particular, in both the stochastic simulations and analytical model, it was computed as a long time average, and in vivo it was computed using averaging cells of the same eccentricity within a tissue. Using other metrics, for example, the 2d nematic order parameter $S}_{2$, which is a quantitative measure of the orientational ordering (de Gennes and Prost, 1993), was prohibitive, since the microtubule density was too low in both our stochastic simulations and experimental images. In contrast, using microtubule angle distribution allowed for accurate comparison between simulations, analytical, and experimental results. Furthermore, it was not possible to perform a systematic analysis of microtubule selforganization dependence on many parameters of microtubule dynamics in vivo. This was due to, first, extremely large number of combinations of individual dependencies, which goes beyond the microtubule network itself; and second, altering the microtubule network has profound consequences on processes relying on microtubules, and thus, on cellular functions. We, therefore, analyzed our system via mathematical modeling and validated the model predictions in vivo.
Various modeling approaches have been used for describing microtubule selforganization. However, many of them are specific to a particular tissue. In plants, it was shown that microtubule zipping strongly affects their selforganization (Tindemans et al., 2010), and that tension can have a nonnegligible effect on stabilizing microtubules (Hamant et al., 2019). In larger cells, such as Drosophila oocytes, microtubule nucleation at the cortex was shown to be important (Khuc Trong et al., 2015). Models that include the hydrodynamic effect of the cytoplasm and molecular motors’ effect on microtubule selforganization are summarized in Shelley, 2016 and Belmonte et al., 2017. Our published stochastic model successfully recapitulates the organization of microtubule networks in various epithelial cells (Gomez et al., 2016). It is a minimal in silico 2dmodel, in which the microtubules are seeded on the cell periphery, grow stochastically to capture the dynamic instability (as in Peskin, 1998), and follow geometric interaction rules.
Here, we use this stochastic model for simulations exploring the typical parameter space of microtubule dynamics, discovering that the average microtubule selforganization is robust. We confirm the robustness in vivo using genetic manipulations of epithelial cells in the model organism Drosophila. Finally, we build a minimal probabilistic model, which accurately predicts the experimental results and reveals that the reason for robustness is the separation of timescales in microtubule dynamics. This model shows that the details of microtubule dynamic instability are irrelevant for microtubule selforganization within their biologically relevant ranges and that the only biological quantity beyond cell shape that affects microtubule alignment is the minusend distribution. Therefore, we demonstrate the extreme robustness of bidirectional quasi2d microtubule selforganization, which can be explained by simple mathematical rules. This suggests the general applicability of our findings to quasi2d microtubule networks, and provides a foundation for future studies.
Results
A conceptual geometric model accurately predicts in vivo alignment
Cells of the Drosophila epidermis elongate during stages 12–15 of embryonic development, changing their eccentricity from 0.7 to 0.98 (Figure 1A and Gomez et al., 2016). As cells elongate, initially randomly oriented microtubules become gradually aligned (Gomez et al., 2016). The simplest thought experiment to visualize how cell elongation translates into microtubule alignment is the following. Imagine a ‘hairy’ unitcircle on the $(x,y)$ plane, where ‘hairs’ are microtubules. Turn it inside out (Figure 1B). The microtubules are randomly pointing inside the ball, representing the absence of microtubule alignment in nonelongated cells; at each microtubule minusend on the cell boundary, the mean microtubule direction is normal to the cell boundary. We now deform both the cell and the filament directions by stretching the cell uniformly from its center in the $y$direction by a factor $b$. This results in an ellipse with eccentricity $e(b)=\sqrt{1{b}^{2}}$, where the minusend positions move proportionally to deformation, the filament directions point towards the cell center and their lengths remain unaltered. Mathematically, the distribution $\rho (\theta )$ of microtubule angles $\theta \in [0,180]$ changes from uniform, $\rho (\theta )=1/180$, to the angle distribution we call hairyball distribution, ${\rho}_{HB}(\varphi )$, which is the inverse Jacobian of the stretching mapping $F:\theta \to \varphi $, where $\theta ,\varphi \in [0,180]$, given by $\mathrm{tan}(\varphi (\theta ))=b\mathrm{tan}(\theta )$,
where $M=180/b$ is the normalization constant. This result gives a surprisingly good agreement with the experiment (Figure 1C), especially considering that this model does not take into account the underlying biological processes, for example, microtubule dynamics. Therefore, while a detailed mathematical model is required to understand how various biological processes control microtubule alignment, the hairyball angle distribution formula in Equation 1 provides a valuable shortcut for the analysis of biological data and parameterizations of microtubule angle distribution, where it can be used to fit data with one parameter – the ‘effective aspect ratio’ $b$, and, therefore, eccentricity.
Stochastic simulations demonstrate robustness of microtubule selforganization for a wide range of parameter values
We now explore how microtubule selforganization depends on dynamics and interactions of individual microtubules in silico. To this end, we use the same model setup as previously published (Gomez et al., 2016), as this stochastic model recapitulated microtubule selforganization observed in vivo. To focus the study on the role of microtubule dynamic instability and interactions and reduce the number of free parameters, we kept the density of the microtubule minusends on the cell boundary uniform. We ran the simulations on cells with fixed shapes, since in vivo the cell shape evolves on a longer timescale (hours) as compared to the time required for the microtubule network to stabilize (several minutes) (Gomez et al., 2016). We chose the cell shape to be an ellipse since we demonstrate below that the averaged experimental cell shape is an ellipse as well. Finally, the cell eccentricity range in the simulations, 0.7 – 0.98, mimicked the experimental one.
To capture the dynamic instability, we model microtubules as follows (Figure 2A). Since the microtubule width (24 nm) is much smaller than the typical cell size (2–10 µm) (Bulgakova et al., 2013), we model microtubules as 1d filaments. They are composed of equal length segments, representing microtubule dimers, whose dynamics is governed by a continuous time Markov chain (Figure 2A; Peskin, 1998; Gomez et al., 2016). The microtubule grows (polymerizes) at the rate $\alpha$, and shrinks (depolymerizes) at the rate $\beta >\alpha $; it switches from the polymerizing to depolymerizing state at the catastrophe rate ${\beta}^{\prime}$; and undergoes the reverse switch at the rescue rate ${\alpha}^{\prime}$. These dynamic instability rates (apart from $\beta$) depend on the concentration of free tubulin dimers in cytoplasm (Walker et al., 1988), which is reported to vary between 30 and 75% of the total tubulin in cells in vivo (Pipeleers et al., 1977; Reaven et al., 1977; Zhai and Borisy, 1994). However, after the microtubule network stabilizes, the total amount of tubulin in microtubules, and, therefore, in the cytoplasm, remains approximately constant, leading to approximately constant dynamic instability rates. Since we investigate the statistics of the microtubule network in steady state (see Materials and methods), we use constant microtubule dynamic instability rates throughout the simulation, and set the same rescue rate for a completely depolymerized microtubule as for any nonzero length microtubule. To account for potentially different microtubule network organizations due to the initial tubulin concentration, we investigate a broad range of parameters. As discussed below, we find that the microtubule angle distribution of stabilized microtubule networks is not sensitive to the parameters of the dynamic instability.
To include microtubule interactions with other microtubules and cell boundaries in the model, we parameterize them using the known parameterization in plants (Dixit and Cyr, 2004; Tindemans et al., 2010) and Drosophila cells (Figure 2A; Gomez et al., 2016), which relies on two parameters – the critical angle, ${\theta}_{c}$, and the probability of catastrophe, ${p}_{cat}$ – as follows. Upon encountering another microtubule at an angle θ, if $\theta \ge {\theta}_{c}$, the growing microtubule either undergoes a catastrophe with probability ${p}_{cat}$ or crosses the microtubule otherwise. If $\theta \le {\theta}_{c}$, it collapses with probability $p(\theta )=\frac{\theta}{{\theta}_{c}}{p}_{cat}$, crosses with probability $p(\theta )=\frac{\theta}{{\theta}_{c}}(1{p}_{cat})$, and bends to change its direction and continues to grow parallel to the existing microtubule otherwise (the microtubule is said to zip). Upon reaching a cell boundary at an angle $\theta \le {\theta}_{c}$, the microtubule zips along it with probability $p(\theta )=1\frac{\theta}{{\theta}_{c}}$, and depolymerizes otherwise (Gomez et al., 2016). These rules for microtubule interactions were originally inspired by the wellestablished induction of catastrophe when a microtubule grows against a barrier (Janson et al., 2003). The role of cell borders as barriers for microtubule growth in epithelial cells is supported by the observations that microtubules buckle at the cell cortex (Singh et al., 2018) and that the microtubule catastrophes at the cell boundaries are angledependent in vivo (Gomez et al., 2016). We envision that the exact critical angle, ${\theta}_{c}$, and the probability of catastrophe, ${p}_{cat}$, depend on the presence of specific crosslinking and motor proteins that promote microtubule bundling and stability (Yan et al., 2013; Takács et al., 2017). Finally, we did not include in the model the effect of microtubules sliding along each other promoted by crosslinking proteins. Sliding has a profound effect on microtubule selforganization in diverse other systems ranging from oocytes to neurons (Zemel and Mogilner, 2008; Lu et al., 2016; Winding et al., 2016). However, we assume this effect to be secondary in our system, since microtubule minusends being anchored at the cell boundary (as opposed to being free) prevents sliding. Instead, the leading cause of microtubule selforganization in our system is that microtubule plusends are dynamic (as opposed to stabilized). This allows microtubules to 'sense’ the scale of the cell: when a microtubule grows toward a cell boundary, upon reaching it, the microtubule either continues growing parallel to it or undergoes a catastrophe.
This setup allows us to investigate a broad range of biologically relevant microtubule dynamics scenarios: in an organism, the dynamic instability parameters are linked to the expression of plus and minusend binding proteins and severing factors, while the interaction parameters ${\theta}_{c}$ and ${p}_{cat}$ are linked to the presence of crosslinking and motor proteins, as described above, and temperaturedependent microtubule rigidity (Kawaguchi and Yamaguchi, 2010). The relation of the model parameter values to their dimensional equivalents is as follows. The typical observed microtubule growth speed is 0.15 µm/s (Gomez et al., 2016). Expressing it as $\alpha \times d\times R$, where $\alpha =1000$ is the nondimensional base growth rate, $d=8.2$ nm is the height of one dimer, and $R$ is the dimensionality coefficient, we find $R$ to be 0.0183 s^{1}. Therefore, the dimensional rates are: the microtubule growth speed ${\alpha}_{dim,speed}=0.15$ µm/s; shrinking speed ${\beta}_{dim,speed}=0.52$ µm/s; rescue rate ${\alpha}_{dim}^{\mathrm{\prime}}=0.07316$ s^{1}; and catastrophe rate ${\beta}_{dim}^{\mathrm{\prime}}=0.01829$ s^{1}.
We varied $(\alpha ,\beta ,{\alpha}^{\prime},{\beta}^{\prime})$ independently, each from 0.5 to 1.5 times the baseline value $(\alpha ,\beta ,{\alpha}^{\prime},{\beta}^{\prime})=(1000,3500,4,1)$. We, therefore, tested the following range of dimensional rates: the microtubule growth speed ${\alpha}_{dim,speed}=0.0750.225$ µm/s; shrinking speed ${\beta}_{dim,speed}=0.260.78$ µm/s; rescue rate ${\alpha}_{dim}^{\prime}=0.03660.1097$ s^{1}; and catastrophe rate ${\beta}_{dim}^{\prime}=0.00910.0274$ s^{1}. These are in the biologically relevant range, since in vivo the parameters of microtubule dynamics depend widely on the cell type, with the reported ranges of growth being 0.05–0.5 µm/s; shrinking 0.13–0.6 µm/s; rescue 0.01–0.17 s^{1}; and catastrophe 0.003–0.08 s^{1 }(Rogers et al., 2002; Komarova et al., 2005; Bulgakova et al., 2013). We further varied $\theta}_{c$ between 0° and 40° and $p}_{cat$ between 0.01 and 0.3, where ${\theta}_{c}=10,\phantom{\rule{thinmathspace}{0ex}}40$ and ${p}_{cat}=0.1,\phantom{\rule{thinmathspace}{0ex}}0.3$ are reported in Figure 2B. As the strength of microtubule interaction increases with $\theta}_{c$, the case ${\theta}_{c}=0$ corresponds to noninteracting microtubules, which would occur, for example, in the absence of crosslinkers.
The simulations gave an unexpected result (Figure 2B) that the microtubule angle distribution varied only slightly, suggesting that in this model microtubule selforganization is robust in individual cells. In particular, it does not depend on the details of dynamic instability and the strength of microtubule interaction, which could range from strong (at ${\theta}_{c}=40$) to nonexistent (at ${\theta}_{c}=0$). Despite the robustness of microtubule selforganization measured via the microtubule angle distribution, we found that some features of the microtubule network varied spatially inside cells (Figure 3). In the absence of zipping (at ${\theta}_{c}=0$), the microtubule density per unit area was the highest near the longer cell sides, and these peaks shifted toward the cell center with increasing zipping (higher ${\theta}_{c}$). Altogether, our model suggests that the microtubule angle distribution for any given cell is robust, while other microtubule network characteristics could vary spatially inside it.
In vivo manipulations of microtubule dynamics and stability alter microtubule density but not alignment
To test in vivo the robustness of microtubule selforganization predicted by the stochastic simulations, we examined how changes in microtubule dynamics and stability affect the organization of subapical microtubules in cells of the Drosophila embryonic epidermis, where microtubules are constrained to the thin 1 µm apical layer of the cell and grow in the plane of the adhesion belt (Gomez et al., 2016). Using genetic manipulations, we could either increase the catastrophe rate ${\beta}^{\prime}$, simultaneously increase the catastrophe rate ${\beta}^{\prime}$ and shrinkage rate β, or reduce the number of minusends, therefore reducing the density of the network and thus encounters and zipping between microtubules. In particular, we increased ${\beta}^{\prime}$ by overexpressing a dominantnegative variant of End Binding protein 1 (EB1DN), a form which increases the number of catastrophes without changing other parameters of microtubule dynamics (Bulgakova et al., 2013; Komarova et al., 2005) or increased both ${\beta}^{\prime}$ and β by overexpressing the protein Spastin, which severs and disassembles microtubules (Bulgakova et al., 2013; RollMecak and Vale, 2005; Komarova et al., 2005). These proteins were overexpressed using the UASGal4 system, in which the Gal4 protein expressed from a tissuespecific promoter induces overexpression of the protein of interest by binding the Upstream Activating Sequence (UAS) (Brand and Perrimon, 1993). Specifically, we used engrailed::Gal4, which drives expression in stripes along the dorsoventral axis of embryos, which correspond to posterior halves of each segment (Figure 4A). In this instance, we avoided abolishing all microtubules by using mild Spastin overexpression (Figure 4B). Overexpression of a CD8Cherry protein, which does not alter microtubules, was used as a control. We also reduced the number of minusends using a null mutation of the minusend capping protein Patronin (Nashchekin et al., 2016), one of the best characterized proteins that protects microtubule minusends (Goodwin and Vale, 2010). We used zygotic mutants  embryos carrying two mutant alleles of Patronin but produced by heterozygous mothers with one wild type allele. Therefore, some Patronin protein is present in homozygous mutant embryos due to maternal contribution  the protein supplied by mothers into eggs, leading to subapical microtubules being reduced but not abolished. Finally, we altered microtubule interactions by manipulating the motor protein Kinesin1 and the microtubuleactin crosslinker Shortstop (Shot). Kinesin1 crosslinks microtubules (Yan et al., 2013), thus its downregulation simulates reduced microtubule zipping. In contrast, Shot stabilizes microtubules and crosslinks them with actin cytoskeleton (Takács et al., 2017). We used the same engrailed::Gal4 to knockdown the Kinesin1 using the expression of interfering RNA (RNAi) against its heavy chain (KhcRNAi), and complemented this knockdown by using zygotic mutants for Khc and shot.
We quantified how the above manipulations altered the organization of the microtubule network in cells by obtaining images of microtubules stained with an antibody recognizing αTubulin, and analyzing them on a cellbycell basis. From these images, we obtained two types of information about microtubule organization (see Materials and methods). First, we determined how the manipulations altered the amount of microtubules in cells by quantifying the percent of the cell area covered by αTubulin signal. This measure is a good proxy for the number of microtubules in cells, but it is also sensitive to local arrangements of microtubules – for example, less clustered microtubules due to reduced crosslinking produce a larger area of αTubulin signal, as shown below. Therefore, the interpretation of changes in αTubulin signal area is done in relation to the known effects of each modification. Second, we determined if the genetic manipulations altered the microtubule alignment. To do this, we determined the direction and magnitude of change of the αTubulin signal at each position within the cell (see Materials and methods and Gomez et al., 2016), and produced the microtubule direction distributions.
We focused on a late stage of Drosophila embryo development (stage 15), as the amounts of protein expressed using UASGal4 system increases, whereas amounts of protein supplied from mothers into zygotic mutants (the maternal contribution) decrease with the progress of embryonic development. The overexpression of both Spastin and EB1DN reduced the area of αTubulin signal in cells (pvalues $p<0.0001$ and $p=0.003$, respectively, Figure 4B–C), consistent with their functions. Similarly, the area of αTubulin signal was reduced in heterozygous Patronin^{+/} embryos in comparison to wildtype controls (pvalue $p=0.02$, Figure 4A–B), and even further reduced in homozygous Patronin^{ /} embryos (pvalues $p=0.0003$ and $p=0.002$ in comparison to wildtype control and heterozygous siblings, respectively, Figure 4B–C). This result is consistent with the dosedependent protection of the microtubule minusends by Patronin. In contrast, zygotic loss of Kinesin1 did not alter the area of αTubulin signal as compared to heterozygous siblings (Figure 4B–C). In cells both hetero and homozygous for the Kinesin1 mutation, this area was the same and significantly smaller than that in controls presented above. The same signal area was measured in two independent experimental repeats (Appendix 1—figure 1), which highlights both the reproducibility of our approach and variability of this measure in response to genetic background, unlike the microtubule angle distribution. In contrast, the KhcRNAi led to a slightly elevated signal area ($p=0.01$, Figure 4B–C), which is in agreement with both the presence of less clustered microtubules and the role of this motor protein in crosslinking. We therefore suggest that the absence of detectable difference in αTubulin signal area between cells hetero and homozygous for the Kinesin1 mutation is either due to rescue by maternally supplied protein or the use of internal controls allowing for more sensitive detection of differences. Finally, although the average αTubulin signal area was not affected in shot mutant cells (Figure 4B–C), we observed that microtubules seemed more likely to be in close proximity of cell boundaries, similar to the organization reported before (Takács et al., 2017). Such organization is predicted to result from reducing the zipping strength in our stochastic simulations (Figure 3) and is consistent with the role of the Shot protein in microtubule crosslinking with actin. Despite the observed changes in area and subcellular distribution of αTubulin signal, and differences between controls, the microtubule angle distributions did not differ between all genotypes and relative to controls in all cases (Figure 4D). Altogether, these results support the robustness of microtubule selforganization despite variable microtubule dynamics and amounts.
To capture a wide range of eccentricities, we used the embryos at different stages of development during which the epidermal cells progressively elongate from eccentricities around 0.7 to 0.98 (stages 12 through 15). To this end, we focused on the genotypes with the strongest effects at stage 15 (Spastin, EB1DN, and CD8Cherry, as a control) and used paired::Gal4, which although leading to milder overexpression than engrailed::Gal4, is expressed in broader stripes along the dorsoventral axis of embryos. Overexpression of Spastin reduced αTubulin signal area in comparison to both neighboring cells, which did not express paired::Gal4, and control cells expressing CD8Cherry (both pvalues $p<0.0001$, Appendix 1—figure 2AB). Similarly, the αTubulin signal area was reduced in embryos homozygous for Patronin^{ /} in comparison to heterozygous siblings (pvalue, $p=0.04$, Appendix 1—figure 2AB). We suggest that the lack of a difference between heterozygous Patronin^{ +/} and the wildtype control observed here might be due to measurements being taken over a broader developmental time, which makes it more difficult to detect changes. Overexpression of EB1DN did not change the area covered by the αTubulin signal per cell (pvalue, $p=0.98$, Appendix 1—figure 2AB). This can be explained by either weaker expression of paired::Gal4 in comparison to engrailed::Gal4, or lower amounts of Gal4 at earlier developmental stages. The microtubule angle distributions for each eccentricity (binned at a particular eccentricity ±0.025) did not significantly differ between all genotypes and relative to controls (Appendix 1—figure 2C). Although only two of the above manipulations affected microtubule density, these results further support that microtubule selforganization is indeed robust in vivo.
Having found robustness in the embryonic epidermis, we sought to test that the same is true in another epithelial tissue – pupal wing. These cells are elongated to a lesser degree than average epidermis cells in the apical plane, but display similar subapical microtubule networks that align along cells’ major axes (Gomez et al., 2016). We examined the effects of overexpressing EB1DN and Spastin in pupal wings, as these had the largest effects on αTubulin signal area in embryos, and used CD8Cherry as a control. We have compared cells in engrailed::Gal4expressing posterior compartments of pupal wings with anterior compartments, which do not express engrailed::Gal4 but are otherwise genetically identical. Expression of EB1DN reduced the area of αTubulin signal, although to a lesser extent than in embryos ($p=0.07$, Figure 5A–B). This reduction was not reflected in any change in the microtubule angle distribution (Figure 5C). The expression of Spastin in pupal wings did not produce a consistent phenotype – some wings looked as in control, whereas in others microtubules were clearly severed (Figure 5A), which lead to no apparent change in the bulk αTubulin signal area ($p=0.36$, Figure 5B). To determine if the microtubule angle distribution was altered in cells with severed microtubules in comparison to nonsevered ones, we compared cells with signal areas of αTubulin lower than 0.45 with cells where it was greater than 0.45. We found no differences between these two groups of cells (Figure 5C). The distributions did not differ from those in cells from anterior compartments, which did not express Spastin, and from cells expressing CD8Cherry as a control. Altogether, our findings demonstrate that microtubule angle distribution is robust to perturbations in dynamics of individual microtubules, microtubule interaction, and the number of microtubules, in two independent epithelial tissues in vivo. As the microtubule angle distributions were produced by averaging cells of the same eccentricity within a tissue, we term this robustness on the tissue scale.
An analytical model shows that microtubule selforganization depends on the cell geometry and minusend distribution
Given that the microtubule angle distribution both in silico and in vivo only weakly depends on microtubule interactions, we propose a minimal mathematical model with noninteracting microtubules, which is analytically tractable. Here, independent microtubules cross upon reaching one another and fully depolymerize upon reaching a cell boundary. Their averaged behavior is the average of 1d behaviors of individual microtubules growing from different positions on the cell boundary. While we expect the microtubules to have an overall alignment along the major cell axis due to the high buildin catastrophe rate at cell boundaries, our goal is to obtain the full microtubule angle distribution.
Our model setup is as follows (Figure 6A). Consider a convex 2d cell with the boundary parameterized by the arclengthcoordinate $\zeta$ increasing in a counterclockwise direction. From microtubule minusends distributed with density $\rho (\zeta )$ on the cell boundary, the microtubules grow into the cell interior at an angle $\theta$ to the cell boundary and at an angle $\phi$ with respect to the horizontal. Note that the cell shape is fully determined by the function $a(\zeta ,\theta )$ – the length of a crosssection that starts at $\zeta$ at an angle $\theta$ with respect to the boundary. When $a(\zeta ,\theta )$ is considered as a function of $(\zeta ,\phi )$, we denote it by $\stackrel{~}{a}(\zeta ,\phi )$ to avoid confusion. Microtubules undergo dynamic instability by switching between the states of growth, shrinking, catastrophe, and rescue at the rates $\alpha$, $\beta$, ${\beta}^{\prime}$, and ${\alpha}^{\prime}$, respectively. After fully depolymerizing, they regenerate at the rate ${\alpha}^{\prime}$ from the same minusend but in a new direction at an angle $\theta$ taken from a uniform distribution on $[0,180]$. Thus, over a large fixed time interval $t\in [0,T]$, a microtubule undergoes a large number $N$ of growth and shrinkage lifetimes, which are separated by periods of average duration $1/{\alpha}^{\prime}$ when the microtubule has zero length.
The first quantity of interest is the microtubule mean survival time. Since both in vivo and in silico, the microtubule angle distribution is lengthweighted, we include the general case of weighting the mean survival time by a function $\gamma (x)$ of microtubule length $x$. Then the mean survival times, $f(x)$ and $g(x)$, of polymerizing and depolymerizing microtubules of length $x$ satisfy
where the terms on the righthand side are the contributions from growing (and shrinking in the $g$ case), switching, and the time increment weighted by $\gamma (x)$, which we specify below. Here, $dt$ is a small timeincrement. Expanding Equation 23 in Taylor series and neglecting terms of the second and higher order in $dt$, we obtain that $f(x)$ and $g(x)$ are governed by
where $p=\frac{{\beta}^{\prime}}{\alpha}\frac{{\alpha}^{\prime}}{\beta}$ and $q=\frac{1}{\alpha}+\frac{1}{\beta}$. We assume that $g(a)=0$, that is, once the microtubule reaches the cell boundary at the length $a(\zeta ,\theta )$, it quickly depolymerizes. Finally, for a zerolength microtubule $g(0)=0$, and hence $h(0)=f(0)$. Note that the only quantity of interest is $f(0)=h(0)$, since it is the lifetime of a microtubule when it starts in a growing state with zero length.
The two choices $\gamma (x)=1$ and $\gamma (x)=x$ give the solutions for the notweighted and the lengthweighted mean survival times, denoted ${f}_{nw}(x)$ and ${f}_{w}(x)$, respectively
For each microtubule minusend location, $\zeta$, the average time between any two regrowths of a microtubule is the sum of the averaged waiting time $1/{\alpha}^{\prime}$ and the average of the nonweighted lifetime over all the growth angles $\theta$
Then the average number of lifetimes with direction $(\phi ,\phi +\delta \phi )$ with respect to the horizontal is $N\delta \phi /\pi$, and their contributions to the lengthintegral is ${f}_{w}(0)N\delta \phi /\pi =({f}_{w}(0)/{T}_{ave})T/\pi \delta \phi$. Integrating it over the cell boundary weighted by the density of minusends ${\rho}_{m}(\zeta )$ and using Equation 8, we obtain the lengthweighted microtubule angle distribution
where $M$ is the normalization constant. The cell crosssection is denoted by either $a(\zeta ,\theta )$ or $\stackrel{~}{a}(\zeta ,\phi )$ depending on its arguments. This analytical prediction matches the stochastic simulations (Figure 6B).
The dependency of the resulting microtubule angle distribution on the biological rates reduces to two parameters
If both of them are small, $p\ll 1$ and $\frac{{\alpha}^{\prime}}{\pi}\left(\frac{1}{\alpha}+\frac{1}{\beta}\right)\ll 1$, as is always observed in biological systems (see above), the microtubule angle distribution formula can be significantly simplified. In particular, for small $p$, the distribution in Equation 9 reduces to
to leading order, and to
when, in addition, $\frac{q{\alpha}^{\prime}}{\pi}\ll 1/{\int}_{0}^{\pi}a(\zeta ,\theta )\mathit{d}\zeta $. This becomes exact in the limit of deterministic microtubules ($\beta =\mathrm{\infty}$, ${\beta}^{\prime}=0$). Note that while $p$ is required to be nonnegative in models of microtubules on an infinite line (Peskin, 1998), our setup does not have this restriction, as the microtubule lifetimes ${f}_{w}(0)$ and ${f}_{nw}(0)$ are positive even for negative $p$. Furthermore, the analytical microtubule angle distribution, Equation 9, is independent of the multiplicative change in the minusend density, ${\rho}_{m}(\zeta )$, which would be absorbed into the normalization constant $M$. Only nontrivial changes to the density of minusends that vary along the cell boundary affect the microtubule angle distribution.
It is required that ${\alpha}^{\prime}\ll \alpha $ and ${\alpha}^{\prime}\ll \beta $ for the second parameter in Equation 10 to be small, and ${\beta}^{\prime}\ll \alpha $ for the first parameter to be small as well. This separation of timescales in microtubule dynamics is observed in vivo, as described above, where the rates of polymerization and depolymerization are much higher than those of catastrophe and rescue. Therefore, the microtubule angle distribution depends only on the cell geometry and the minusend distribution, and the underlying reason for that is the separation of timescales in microtubule dynamics.
The analytical model accurately predicts microtubule selforganization given both the experimental cell shape and distribution of microtubule minusends
To validate our analytical model’s predictions, we sought to compute the analytical microtubule angle distributions in Equation 9 using experimental cell shape and minusend density. Note that Equation 9 predicts that the microtubule angle distributions is robust in each individual cell, where each distribution is determined by the cell geometry and the distribution of minus ends. For example, a cell corner at multicellular junctions (where three or more cells contact) leads to ‘corners’ in microtubule angle distributions. However, cell shapes in tissue are highly variable, even when they have the same eccentricities. We, therefore, sought to perform comparison between the model and experiment using the averaged values of the microtubule minus end locations and cell shape.
First, we determined the localization of microtubule minusends. Since Patronin localizes at the microtubule minusends, we analyzed its distribution in epithelial cells in the Drosophila embryonic epidermis using PatroninYFP (Nashchekin et al., 2016). As expected, PatroninYFP is mostly localized at the cell boundaries with few speckles inside cells (Figure 7A). We quantified the distribution of PatroninYFP at the cell boundaries by measuring its asymmetry, namely the ratio of PatroninYFP average intensity at dorsoventral borders to that of anteriorposterior borders (see Materials and methods, Figure 7B, and Bulgakova and Brown, 2016). The asymmetry of Patronin distribution was a linear function of the cell eccentricity (Figure 7C), suggesting that Patronin becomes enriched at the dorsoventral boundaries as the embryo develops and cells elongate. Additionally, when comparing boundaries in cells with similar eccentricities (Stage 15 embryos only), the intensity of PatroninYFP was decreasing with the border angle relative to the anteriorposterior axis of the embryo (Figure 7D). Several lines of evidence support that the observed enrichment of PatroninYFP at the dorsoventral boundaries is due to the asymmetry of cellcell adhesion in these cells. Indeed, microtubule minusends were shown to be tethered by cellcell adhesion in some epithelial cells (Meng et al., 2008). Concurrently, Ecadherin, the key component of cellcell adhesion in Drosophila embryonic epidermis, is asymmetrically distributed in stage 15 embryos (Bulgakova et al., 2013), with enrichment at the dorsalventral borders similar to that of Patronin. Finally, asymmetries of both Patronin and Ecadherin progressively increase from stage 12 to stage 15 of Drosophila embryonic development (Figure 7C and N.A.B. personal communication).
To use this in our analytical model, we used leastsquares to simultaneously fit the asymmetry data with a linear function of eccentricity, and the normalized intensity of PatroninYFP with an exponential function of the cell border angle. We imposed a constraint that the asymmetry value at eccentricity 0.98 (Stage 15) is the same as the normalized intensity of PatroninYFP at the dorsoventral border (border angle ${0}^{o}$). The resulting formula used in the analytical model is
where $\psi$ is the cell border angle with respect to the horizontal, ${C}_{1}=0.4144$ and ${C}_{2}=0.0231$.
Next, we determined the average cell shape. In tissue, each cell has a unique shape, and cells with the same eccentricities may differ significantly in their geometry. Therefore, to test and validate the analytical solution of microtubule selforganization we have generated masks of epithelial cells in the Drosophila embryonic epidermis, which provided us with coordinates of cell boundaries (see Materials and methods). Dividing all cell shapes into groups by eccentricity, we computed the average cell shape for each group as follows. First, the cells were recentered to have their centers of mass at the origin. We then rotated them so that they are elongated along the vertical axis (the direction of elongation is the first singular vector, see Materials and methods). Finally, we rescaled all the cells to have unit area. The average of the distance from the center of mass in a particular direction to the cell boundary traced the boundary of the averaged cell. Surprisingly, we found that the average cell shape for a given cell eccentricity is an ellipse (Figure 8A).
The analytical microtubule angle distribution, Equation 9, computed on an averaged cell shape using the experimental minusend density, Equation 13, gives surprising agreement with experimental in vivo data (Figure 8B). This agreement is better for the case of asymmetric minusend distribution, comparing to the uniform one, which supports our prediction that the minusend distribution does indeed influence the microtubule selforganization.
Discussion
Here, we present several novel findings that describe the fundamental rules underlying selforganization of microtubule networks in epithelial cells. Firstly, we have shown robustness of microtubule selforganization both in silico and in vivo; secondly, in addition to the known importance of cell shape for microtubule organization (Gomez et al., 2016), our minimal analytical model predicted the importance of the asymmetric minusend distribution as the only other impactful parameter, which we then confirmed in vivo in Drosophila epidermal cells; and finally, we have shown that this robustness originates from the intrinsic separation of timescales of microtubule dynamic instability.
Robustness on the tissue scale
Biologically, the discovered robustness of microtubule organization makes perfect sense, given the fundamental importance of microtubule functions in a cell. Most intracellular trafficking events require microtubules for the delivery of various cellular components to their relevant biological locations by motor proteins (HammAlvarez, 1998; Apodaca, 2001). This process must be reliable, as mislocalization of cellular components leads to cell death or disease (Levy et al., 2006; Lopes et al., 2010; Baek et al., 2018). However, this delivery mechanism is highly stochastic, given the microtubule dynamic instability and the dynamics of molecular motors (Kolomeisky and Fisher, 2007; Brouhard, 2015; Goodson and Jonasson, 2018). We suggest that it is the average microtubule angle distribution, which is likely to guide the net outcome of intracellular trafficking and focus on this distribution in this work.
We have shown that going from the subcellular to the tissue scale, the microtubule organization becomes more and more robust. On the subcellular level, some features of microtubule networks are variable, for example, the timeaveraged microtubule density per unit area in our stochastic simulations depends on the strength of microtubule interactions. In particular, in systems with strong microtubule interactions, the microtubule density plateaued at the cell center, while in the absence of microtubule interactions, the microtubule density was higher near the cell boundaries, which we also observed experimentally in vivo in shot mutant cells. Moving to the scale of individual cells, our stochastic simulations and the analytical model showed that the averaged microtubule angle distribution is already robust. This is because the average microtubule selforganization depends to leading order only on the slowly evolving parameters, such as the cell shape and the density of the minusends on the cell boundary. However, microtubule selforganization is distinct in individual cells, where individual cellshape features such as cellcorners at the multicellular junctions are reflected in the microtubule angle distribution, since the latter is proportional to the integral over the cell boundary of the squared crosssection of the cell. Finally, on the multicellular level, the microtubule selforganization is robust on the tissuescale, since the effect of individual cell shape variability averages out. In particular, while most of the cells of a given eccentricity in biological tissues are polygons, we found that the average cell shape of a particular eccentricity is an ellipse.
One of the findings of our analytical model is that the robustness of microtubule selforganization exists only as long as the microtubule dynamics exhibits a separation of timescales, ${\alpha}^{\prime},{\beta}^{\prime}\ll \alpha ,\beta $, a rule which is observed in all published data about microtubule dynamic parameters (Shelden and Wadsworth, 1993; Rogers et al., 2002; Komarova et al., 2002; Komarova et al., 2009; Dhonukshe and Gadella, 2003; Zilberman et al., 2009; Zhang et al., 2011; Bulgakova et al., 2013; Shaebani et al., 2016). Our mathematical model shows that if this rule is not observed, the microtubule organization becomes sensitive to changes in these rates. As these rates depend on multiple internal and external factors, such as changes in gene expression and temperature, the microtubule organization would be unpredictable in a cell in a biological tissue. Therefore, breaking this rule will impair cellular function over time, which suggests that any mutations that led to such change were likely to cause cell lethality and did not fix in evolution.
Hairyball distribution
We demonstrated that the microtubule angle distribution is accurately predicted by the hairyball distribution. The reason for this excellent agreement with the experimental data remains an open question, as we were unable to show that hairyball, which is a result of a conceptual 0th order model, has any relation to our analytical distribution, neither as an approximation nor as a limiting case. As presented, the hairyball distribution does not include the effect of nonuniform microtubule minusend distribution. We found that including it (as a multiplicative factor in Equation 1) does not significantly change the agreement with the experimental data. We suggest that the best use of the hairyball distribution is as a simple adhoc formula to parameterize the microtubule angle distribution in cells up to eccentricity $0.95$ using the single ‘effective aspect ratio’ parameter. This could prove useful in investigating, for example, correlations and interdependencies between the microtubule network organization (e.g. their overall direction and spread) and dynamic intracellular processes (e.g. signaling and transport).
The importance of microtubule interactions
While we show that microtubule selforganization does not depend on microtubule interaction, we admit that this is true for the measure of selforganization being the timeaveraged microtubule angle distribution. Here, such effects of microtubule interactions as zipping and bundling disappear due to timeaveraging of the dynamics. However, from the biological point of view, what matters for an organism is the longterm behavior, because most of the processes such as microtubulebased transport occur on much longer timescales than microtubule network rearrangements (Jankovics and Brunner, 2006; Bulgakova et al., 2013; Iyer et al., 2019). Therefore, we hypothesize that it is the microtubule angle distribution that affects tissue behavior on long timescales. We further hypothesize that such effects as bundling and zipping will affect shortterm intracellular transport. For example, the presence of Spastin, the microtubule severing protein, leads to a change in the delivery of the Ecadherin, the protein responsible for the cellcell adhesion delivered along the microtubule network (Bulgakova et al., 2013). A more detailed modeling approach that includes the effect of microtubulemicrotubule interaction on intracellular transport is outside the scope of this article and will be considered in a separate publication.
Our in vivo experiments show that microtubulemicrotubule and microtubuleactin crosslinkers such as Kinesin1 and Shot, which are known to alter microtubule dynamics and local organization (Figure 2; Jolly et al., 2010; Takács et al., 2017; Drechsler et al., 2020), do not alter robustness of average microtubule selforganization. We have recently discovered that molecules localized inhomogeneously in the plane of microtubules (e.g. actin stress fibers) can alter microtubule angle distribution by reorienting the microtubule angle distribution away from the cell’s major axis (Delon and Brown, 2009; Płochocka et al., 2019). However, the average microtubule selforganization in homogeneous environments remains robust in all tested cases.
Our system is a particular but generalizable scenario
We suggest that our findings are applicable beyond apical microtubules in the Drosophila embryonic epidermis, the dynamics of which is quasi2d. Previously, we demonstrated that a similar relationship between microtubule organization and cell shape is observed in other Drosophila epithelia, including cells in pupal wings and ovaries (Gomez et al., 2016). Here, we conducted in vivo experiments in multiple different genotypes in two tissues: embryonic epidermis and pupal wings. Since the tissue and genetics were varied in these scenarios, the microtubule growth and interaction parameters were altered as well. While this led to noticeable differences in such network properties as microtubule density in different controls and local microtubule organization, the average microtubule selforganization reflected by the microtubule angle distribution remained robust on the tissue scale. There are multiple other instances in which maturing and differentiating epithelial cells develop an apical microtubule meshwork, including cells of mammalian airways and even cells in culture (Gilbert et al., 1991; Herawati et al., 2016; Takeda et al., 2018). Our findings are likely to hold true in these systems as long as in them the microtubules are anchored, the plusends are dynamic, and the separation of timescales holds. Furthermore, these rules are likely to apply to squamous cells, which, despite having specialized apical microtubules, have such a small cell depth that microtubules are constrained within a thin plane similar to that in our experimental model (Gomez et al., 2012; Pope and Harris, 2008). The validation of our findings in other cell types and other evolutionary divergent organisms, as well as how the discovered robustness of microtubule selforganization ensures the reliability of intracellular transport, are important questions for future research.
Materials and methods
Computing the eccentricity of a cell
Request a detailed protocolWe uniformly distributed points inside the experimental cell boundary data and used the singular value decomposition on the resulting dataset. The eccentricity then is $ecc={(1{(a/b)}^{2})}^{1/2}$, where $a<b$ are the singular values.
Derivation of the hairyball distribution
Request a detailed protocolUpon stretching a circular cell by a factor of $b>1$ in the vertical into an ellipse with eccentricity $ecc={(11/{b}^{2})}^{1/2}$, the angles $\theta$ with uniform distribution $R(\theta )=1/180$ are mapped into the angles $\psi$ with the distribution $\rho (\psi )$ such that $R(\theta )=J\rho (\psi )$, where $J$ is the Jacobian of the stretching map. Since $R(\theta )=const$, $\rho (\psi )\propto {J}^{1}={\left({\mathrm{sin}}^{2}(\psi )+{b}^{2}{\mathrm{cos}}^{2}(\psi )\right)}^{1}$.
Parameters of stochastic simulations
Request a detailed protocolWe discretized the cell boundary with 180 points, and uniformly placed 200 microtubule minusends along the boundary. The short axis of the cell was kept constant at 60 numerical dimer lengths. The critical angle, ${\theta}_{c}$, was varied in the range (0, 10, 20, 30, 40) degrees, and the probability of catastrophe, ${p}_{cat}$, was varied in the range (0.01, 0.1, 0.2, 0.3). The microtubules zipping along cell boundaries were imposed to undergo catastrophe upon reaching the tip of the ellipse or or a sharp cell corner. This mirrors the scenario in vivo that due to the high stiffness (Gittes et al., 1993), microtubules undergo catastrophe upon reaching an acute cell corner (Figure 1A). All the microtubule angle distributions were computed as a timeaveraged microtubule dimer angle histograms. For each parameter set, 500 simulations were run to the nondimensional time $T=10$, which corresponds to 550 s. At the start of each experiment, all the microtubule lengths were set to zero. During the simulation, the angles with respect to the horizontal of all the dimers of nonzero length microtubules were saved at regular time intervals. We binned these angles into 180 onedegree bins, and then averaged the resulting histograms over the last 2.5 nondimensional timeunits of the simulation, which corresponds to 250 timepoints. The in silico model is available from the authors upon request.
Fly stocks
Request a detailed protocolpaired::Gal4, engrailed::Gal4, UAS::CD8Cherry, UAS::KhcRNAi, Khc^{8}, (Bloomington stocks 1947, 30564, 27392, 35770, and 1607, respectively), UAS::EB1DN (Bulgakova et al., 2013), UAS::Spastin (Sherwood et al., 2004), shot^{3} (gift from K.Röper), Patronin^{05252}, PatroninYFP (Nashchekin et al., 2016), and Ubip63E::EcadGFP (Kyoto DGGR stock 109007). patronin, Khc, and shot mutants are depicted by ^{} in the text. The flies, embryos, and pupae were kept at 18°C.
Embryo fixation, pupal wing dissection, and antibody staining
Request a detailed protocolThe embryos were fixed as described in Gomez et al., 2016. In brief, the staged embryos were dechorionated in 50% bleach for 4 min, and then fixed in 1:1 10% formaldehyde (methanol free, $\mathrm{\#}18814$, Polysciences Inc) in PBS:heptane for 20 min at room temperature (RT) and postfixed/devitellinized for 45 s in 1:1 icecold methanol:heptane. Finally, embryos were washed three times in icecold methanol, kept in methanol between 6 and 24 hr at −20°C, rehydrated in 1:1 PBS with 0.3% Triton X100 (T9284, Sigma):methanol, and washed one time in PBS with 0.3% Triton X100. Rehydrated embryos were blocked for 2 hr in 5% Normal Goat Serum (ab7481, Abcam) in PBS with 0.3% Triton X100.
Prepupa individuals were collected and aged for 24 hr at 25°C. Pupae were fixed in 10% paraformaldehyde (PFA, #R1026, Agar Scientific) after external cuticle removal for 50 min at room temperature. Then, wings were removed from the carcass, released for their cuticles and postfixed in fresh 10% PFA for 10 min. Wings were then blocked for 1 hr at room temperature in 10% Native Goat Serum in PBS with 0.3% Triton X100.
Primary antibody incubations were carried out overnight at 4°C. Primary antibodies used were mouse antiαTubulin 1:1000 (T6199, Sigma), and rat antiEcadherin 1:200 (DCAD2, Developmental Studies Hybridoma Bank). In embryos expressing KhcRNAi and all pupal wings, cell boundaries were labelled by native fluorescence of EcadGFP instead of antibody. Incubation with secondary antibody was performed for 2 hr at 25°C. Alexa Fluor fluorophore Alexa Fluor 488, 594 and 647coupled secondary antibodies (Jackson ImmunoResearch) were used in 1:300. Finally, embryos and wings were mounted in Vectashield (Vector Laboratories).
Image acquisition
Request a detailed protocolAll images were acquired at RT (2022°C). For quantification of microtubule selforganization, we acquired 16bit depth images on the Zeiss AiryScan microscope, using a 60x objective lens. For embryos, zstacks of seven sections with 23.5 px/μm in XY resolution and 0.38 µm distance between zsections were taken. For pupal wings, zstacks consisted of five sections with 23.5 px/μm in XY resolution and 0.185 µm distance between steps. All processing was done at 6.5 power in the ZEN software. For an analysis of PatroninYFP distribution, an upright confocal microscope (FV1000; Olympus) using 60 x 1.42 NA oil PlanApoN objective lens was used. 16bit depth images were taken at a magnification of 12.8 px/μm with 0.38 µm between zsections using Olympus Fluoview FVASW.
The subapical domain of the cell was determined using the Ecadherin junctional signal as a reference for its basal limit, and the absence of αTubulin signal as a reference for its apical limit. All imaging was done on dorsolateral epidermal cells, which excluded the leading edge (first row) cells, given its different identity to the rest of the dorsolateral epidermal cells. Pupal wings were imaged in their anterior and posterior compartments, the latter detected by expression of mCherry protein. Average projections were made using Fiji (http://www.fiji.sc) for measuring signal area and PatroninYFP distribution. Figures were assembled using Adobe CS3 Photoshop and Illustrator (http://www.adobe.com). The processing of images shown in figures involved adjusting gamma settings.
Image analysis to quantify microtubule organization
Request a detailed protocolWe used the same workflow as described in the Gomez et al., 2016. In short, the Ecadherin signal from a max intensity projection was used to obtain cell outlines using Packing Analyser V 2.0 software (Aigouy et al., 2010). These cell outlines were used to identify each cell as an individual object and fit it to an ellipse to calculate eccentricity and the direction of the ellipse major axis using a script in MATLAB R2017b (Mathworks, http://www.mathworks.co.uk). The αTubulin signal within each cell was filtered using the cell outline as a mask, and the magnitude of the signal according to its direction analyzed by convolving the filtered αTubulin signal with two $5\times 5$ Sobel operators (Gomez et al., 2016). The resulting magnitude gradient and direction gradient were integrated into a matrix to assign each pixel a direction and magnitude of change in the intensity of αTubulin signal. To reduce the noise, the pixels with a magnitude less than 22% of the maximum change were discarded. The remaining pixels were binned with respect to their direction gradient (bin size $=4$ degrees). The resulting histogram was normalized. The script is available at https://github.com/nbul/Cytoskeleton; Płochocka, 2021; copy archived at swh:1:rev:f94200246cc2389e3c165b5c55bb1ed87ebf8d25.
To calculate the area of αTubulin signal within each cell, a custommade MATLAB script was used. The average projected images were adjusted such that 0.5% of the data with the lowest intensities were set to black and the 0.5% of the data with the highest intensities were set to white in order to compare datasets across genotypes and compensate for potential variability of antibody staining and laser power. The threshold was then calculated using Otsu’s method, and a binary image was created using the calculated threshold multiplied by 0.7. This multiplication parameter was determined empirically by testing images from different experiments. Finally, the percent of pixels above the threshold was calculated for each cell. The script is available at https://github.com/nbul/Cytoskeleton.
Image analysis to quantify PatroninYFP distribution
Request a detailed protocolThe average projection of each zstack was done using Fiji software, and segmented using Packing Analyzer v2.0. The resulting binary images with coordinates of cellcell borders and vertices were used together with unprocessed average projection by custommade MATLAB scripts to extract the following values: apical cell area, cell eccentricity, average cell orientation within the image, the orientation of individual cellcell borders relative to average cell orientation, and mean intensity of each individual border. Only the cells that were completely within the image were taken for quantification. Bristle cells were excluded by their size. Only the borders that are between two cells that were completely within the image were quantified, and the borders adjacent to bristle cells were excluded. The background signal was determined by binarizing images with an adaptive threshold, which uses local firstorder image statistics around each pixel and is very efficient in detecting puncta. The mean background signal of cells that were completely within the image was subtracted from the mean intensities of cellcell borders. The values for each type of border within single embryos were averaged, and the average intensity of borders with 40–90° orientation was divided by the average intensity of borders with 0–10° orientation to produce a single value of asymmetry for each embryo. Finally, to produce distribution of signal intensity by angle, borders of all embryos at stage 15 of development were pulled and binned using 10° bins. The script is available at https://github.com/nbul/Intensity; Bulgakova and Brown, 2016.
Statistical analysis
Request a detailed protocolAll data was analyzed using Graphpad Prism 6.0c (http://www.graphpad.com). Samples from independent experiments corresponding to each genotype were pooled and tested for normality with the ShapiroWilk test. The αTubulin signal area was analyzed using ANOVA with a posthoc ttest.
Appendix 1
Data availability
At the time of the publication, all the biological data is available on https://doi.org/10.7488/ds/2642.

Edinburgh DataShareData file for the submitted publication "Robustness of the microtubule network selforganization in epithelia.https://doi.org/10.7488/ds/2642
References

Control of microtubule organization and dynamics: two ends in the limelightNature Reviews Molecular Cell Biology 16:711–726.https://doi.org/10.1038/nrm4084

A theory that predicts behaviors of disordered cytoskeletal networksMolecular Systems Biology 13:941.https://doi.org/10.15252/msb.20177796

Targeted gene expression as a means of altering cell fates and generating dominant phenotypesDevelopment 118:401–415.

Dynamic instability 30 years later: complexities in microtubule growth and catastropheMolecular Biology of the Cell 26:1207–1210.https://doi.org/10.1091/mbc.E13100594

Dynamic microtubules produce an asymmetric EcadherinBazooka complex to maintain segment boundariesJournal of Cell Biology 201:887–901.https://doi.org/10.1083/jcb.201211159

Drosophila p120catenin is crucial for endocytosis of the dynamic EcadherinBazooka complexJournal of Cell Science 129:477–482.https://doi.org/10.1242/jcs.177527

Active diffusion in oocytes nonspecifically centers large objects during prophase I and meiosis IJournal of Cell Biology 219:e201908195.https://doi.org/10.1083/jcb.201908195

Cytoskeletal selforganization in neuromorphogenesisBioArchitecture 4:75–80.https://doi.org/10.4161/bioa.29070

The integrin adhesion complex changes its composition and function during morphogenesis of an epitheliumJournal of Cell Science 122:4363–4374.https://doi.org/10.1242/jcs.055996

Optical flow analysis reveals that Kinesinmediated advection impacts the orientation of microtubules in the Drosophila oocyteMolecular Biology of the Cell 31:1246–1258.https://doi.org/10.1091/mbc.E19080440

Microtubulebased transport  basic mechanisms, traffic rules and role in neurological pathogenesisJournal of Cell Science 126:2319–2329.https://doi.org/10.1242/jcs.115030

Molecular motors: directing traffic during RNA localizationCritical Reviews in Biochemistry and Molecular Biology 46:229–239.https://doi.org/10.3109/10409238.2011.572861

Flexural rigidity of microtubules and actin filaments measured from thermal fluctuations in shapeJournal of Cell Biology 120:923–934.https://doi.org/10.1083/jcb.120.4.923

Tao controls epithelial morphogenesis by promoting fasciclin 2 endocytosisJournal of Cell Biology 199:1131–1143.https://doi.org/10.1083/jcb.201207150

Microtubule organization is determined by the shape of epithelial cellsNature Communications 7:13172.https://doi.org/10.1038/ncomms13172

Microtubules and MicrotubuleAssociated proteinsCold Spring Harbor Perspectives in Biology 10:a022608.https://doi.org/10.1101/cshperspect.a022608

Are microtubules tension sensors?Nature Communications 10:2360.https://doi.org/10.1038/s4146701910207y

Molecular motors and their role in membrane trafficAdvanced Drug Delivery Reviews 29:229–242.https://doi.org/10.1016/S0169409X(97)000811

Multiciliated cell basal bodies align in stereotypical patterns coordinated by the apical cytoskeletonJournal of Cell Biology 214:571–586.https://doi.org/10.1083/jcb.201601023

Dynamic instability of microtubules is regulated by forceJournal of Cell Biology 161:1029–1034.https://doi.org/10.1083/jcb.200301147

Which way to go? cytoskeletal organization and polarized transport in neuronsMolecular and Cellular Neuroscience 46:9–20.https://doi.org/10.1016/j.mcn.2010.08.015

Temperature dependence rigidity of nontaxol stabilized single microtubulesBiochemical and Biophysical Research Communications 402:66–69.https://doi.org/10.1016/j.bbrc.2010.09.112

Molecular motors: a theorist's perspectiveAnnual Review of Physical Chemistry 58:675–695.https://doi.org/10.1146/annurev.physchem.58.032806.104532

Cytoplasmic Linker proteins promote microtubule rescue in vivoJournal of Cell Biology 159:589–599.https://doi.org/10.1083/jcb.200208058

EB1 and EB3 control CLIP dissociation from the ends of growing microtubulesMolecular Biology of the Cell 16:5334–5345.https://doi.org/10.1091/mbc.e05070614

Mammalian end binding proteins control persistent microtubule growthJournal of Cell Biology 184:691–706.https://doi.org/10.1083/jcb.200807179

Dysfunction of heterotrimeric kinesin2 in rod photoreceptor cells and the role of opsin mislocalization in rapid cell deathMolecular Biology of the Cell 21:4076–4088.https://doi.org/10.1091/mbc.e10080715

Selforganisation and forces in the microtubule cytoskeletonCurrent Opinion in Cell Biology 15:118–124.https://doi.org/10.1016/S09550674(02)000145

Optimal dynamic instability of microtubulesDocumenta Mathematica, Extra Volume ICM 3:633–642.

A sensitive method for measuring polymerized and depolymerized forms of tubulin in tissuesJournal of Cell Biology 74:341–350.https://doi.org/10.1083/jcb.74.2.341

SoftwareCytoskeleton, version swh:1:rev:f94200246cc2389e3c165b5c55bb1ed87ebf8d25Software Heritage.

Quantitative analysis of tubulin and microtubule compartments in isolated rat hepatocytesJournal of Cell Biology 75:731–742.https://doi.org/10.1083/jcb.75.3.731

Drosophila EB1 is important for proper assembly, dynamics, and positioning of the mitotic spindleJournal of Cell Biology 158:873–884.https://doi.org/10.1083/jcb.200202032

The dynamics of microtubule/MotorProtein assemblies in biology and physicsAnnual Review of Fluid Mechanics 48:487–506.https://doi.org/10.1146/annurevfluid010814013639

Survival of the aligned: ordering of the plant cortical microtubule arrayPhysical Review Letters 104:058103.https://doi.org/10.1103/PhysRevLett.104.058103

Organization of Noncentrosomal microtubules in epithelial cellsCell Structure and Function 41:127–135.https://doi.org/10.1247/csf.16015

Expansion and polarity sorting in MicrotubuleDynein bundlesProgress of Theoretical Physics Supplement 173:17–25.https://doi.org/10.1143/PTPS.173.17

Motorinduced sliding of microtubule and actin bundlesPhysical Chemistry Chemical Physics 11:4821–4833.https://doi.org/10.1039/b818482h

Quantitative determination of the proportion of microtubule polymer present during the mitosisinterphase transitionJournal of Cell Science 107 ( Pt 4:881–890.

Regulation of microtubule dynamics by inhibition of the tubulin deacetylase HDAC6Journal of Cell Science 122:3531–3541.https://doi.org/10.1242/jcs.046813
Decision letter

Raymond E GoldsteinReviewing Editor; University of Cambridge, United Kingdom

Anna AkhmanovaSenior Editor; Utrecht University, Netherlands
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
Through a combination of experiments, theory, and computation, this paper uncovers several general rules that determine the relationship between cell shape and microtubule orientational order. It exploits the separation of time scales between fast microtubule dynamics (e.g. dynamic instability) and the much slower processes involved in cell shape changes to develop a minimal model that captures a wide range of experimental observations.
Decision letter after peer review:
Thank you for sending your article entitled "Robustness of bidirectional microtubule network selforganization" for peer review at eLife. Your article is being evaluated by Anna Akhmanova as the Senior Editor, a Reviewing Editor, and three reviewers.
Summary:
Plochocka et al., present a paper that deals with the affect of microtubule polymerization/depolymerization kinetics on the orientation of filaments in large arrays inside a closed cellular geometry. The authors find that the distribution of MT orientations is dependent on cell shape but independent of the kinetic parameters using both experiments and a geometric model of the MT array. Because of this independence, the authors write that this is an archetypal example of robustness in a complex biological system.
Essential revisions:
1) The authors use a purely geometric model of MT growth and the effects on MT kinetics. We are concerned that this is a step too far. First, crosslinking proteins and molecular motors play a huge part in the assemble of MT arrays and these are not considered at all in this work. E.g., does the orientation distribution change if motors are inhibited in the experiments? For the role of motors, see for example:
Zemel, Assaf and Mogilner, (2008).
Zemel, Assaf and Mogilner, (2009).
2) The authors' model treats the MTs as noninteracting which seems incorrect. In addition, it seems that the model is essentially guaranteed to give an alignment along the long axis of an elliptical enclosure because the catastrophe rate is governed by the contact angle between the MT and boundary. The particular form of that rate needs to be explained and justified better. The authors state that the simulation results are insensitive to the alphas and betas, but what about theta_{C}? What about proteins that lie on the cell surface and potentially affect MT kinetics in a spatial manner?
3) The authors analyze the MT organization only using the angular distribution across the whole cell. But in the cell, and in the more detailed simulation model, they could look at the local nematic order, and plot its spatial distribution. Surely this will reveal more details, which the detailed simulation should recover better than the "hairyball" model. Similarly, the spatial distribution of the zipped MTs, of the locations of catastrophes etc., as function of eccentricity.
4) The cells are approximated as ellipses, similar to the "spherical cow" simplification. However, cells are usually treated as polygons: do the multicell Yjunctions, play a role in the MT organization?
5) The assumption in the "hairyball" model that the MTs remain pointing to the ellipse center is not clear and not explained. is this an effective way to introduce the MTMT interactions that are missing from this model?
6) In Figure 6 the multiple lines for the cellshape/ellipse are difficult to discern. We suggest to plot them separately.
Except for the highest eccentricity, the simplest "hairyball" model fits best the experiments. Why? and does it include the nonuniform nucleation of MTs along the membrane?
7) We are concerned that the case for robustness is overstated. The authors should provide justification for the claims that their results (given the caveats above) would be generalizable to other kinds of MT array geometries (as claimed in the discussion). And it is important to emphasize what is being analyzed. The quantities of interest are the AVERAGED distributions, which by definition are taken on time scales long compared to the underlying microscopic dynamics.
As a general principle it is clear that once there is a separation of time scales
then the slow degrees of freedom dominate, and if reorientations occur at the boundary the law of reflection used will dictate the results. That there should be a strong dependence on the cell geometry (i.e. eccentricity) was already implicated in the work of Khuc Trong et al., (where it was MTs emanating from the boundary) and would follow from the dominance of the wall collisions.
In that sense, the present results are reminscent of much older work on wave pattern formation in the Faraday instability Gluckman et al., (1993) in which time averaging of the chaotic patterns revealed mean values that reflected the geometry of confinement. The work also reminds us of that of Dumais et al., on the relationship between cell geometry and cell division planes in plant cells, where microtubule orientations also play a role (Besson and Dumais, (2011)). But it seems to us that a claim of generality must be made after having demonstrated that the issues neglected above do not change the results.
https://doi.org/10.7554/eLife.59529.sa1Author response
Essential revisions:
1) The authors use a purely geometric model of MT growth and the effects on MT kinetics. We are concerned that this is a step too far. First, crosslinking proteins and molecular motors play a huge part in the assemble of MT arrays and these are not considered at all in this work. E.g., does the orientation distribution change if motors are inhibited in the experiments? For the role of motors, see for example:
Zemel, Assaf and Mogilner, (2008).
Zemel, Assaf and Mogilner, (2009).
We thank the reviewers for raising this important point. Indeed, we are aware of the role of molecular motors and crosslinkers in MT organization. While the model we use in stochastic simulations is purely geometrical, it includes the effect of crosslinking proteins via a largescale parameterization of MTMT interactions. It is done by introducing the probability of a polymerizing MT zipping along the MT it encounters, which depends on the angle between these MTs. Here, the concentration and the nature of the crosslinking proteins affect theta_{c} and p_{cat}. We show in Figure 2 that varying these two parameters does not affect the average MT angle distribution, although it affects the distribution of MT density per unit area inside individual cells (see the new Figure 3). The case of theta_{c}=0 corresponds to noninteracting MTs and could be interpreted as a model of an experiment in the absence of crosslinkers.
We expect that in the studies mentioned in the question, the effect of the crosslinkers is much stronger than in our system. In these papers, due to the fact that the MTs are free (as opposed to anchored) and stabilized (as opposed to dynamic) with either free or periodic boundary conditions, the crosslinkers play a dominant role in the MT organization. However, in our system, the MT minusends are anchored on the cell boundary, while the MT plusends are dynamic. Anchoring inhibits MTs sliding along each other in bundles, and dynamic plusends allow MTs to “sense” the scale of the cell by growing towards another cell boundary and, upon reaching it, either growing parallel to it or collapsing. Due to these two effects, the role of crosslinkers in our system is secondary, which we now support by in vivo experiments. For this reason, we did not directly include sliding as an additional degree of freedom in our mathematical model. We now mention the roles of crosslinkers and molecular motors in the Introduction; the clarifications of how actions of motors and crosslinkers are indirectly incorporated into the model; and the aforementioned interpretation of the results with theta_{c}=0 in subsection “Stochastic simulations demonstrate robustness of microtubule selforganization for a wide range of parameter values”.
Supporting robustness of MT selforganization and in agreement with our model, we show now that in vivo the loss of zygotic kinesin1 and its knockdown with RNAi  the motor protein known to “bundle” MTs (e.g., Yan et al., 2013)  did not affect MT angle distribution in our experimental system. Although testing all possible candidates is outside of the scope of this paper, we also demonstrated that the loss of zygotic MTactin crosslinker Shortstop  a protein known to regulate MTs in the Drosophila epidermis (Takacs et al., 2017)  similarly did not affect MT angle distribution, further substantiating our findings. These new results are included in subsection “in vivo manipulations of microtubule dynamics and stability alter microtubule density but not alignment”, in the new Figure 4, and in the Discussion section.
(2.1) The authors' model treats the MTs as noninteracting which seems incorrect.
We would like to highlight that the detailed stochastic simulations were performed with the interacting MTs, where the parameterization of the MTMT interaction was inspired by that used in plants (Tindemans et al., 2010; Allard et al., 2010). Furthermore, the simplified analytical model with noninteracting MTs was only proposed after we discovered that the MT angle distribution is independent of the MTMT interactions in the detailed simulations and confirmed this in vivo.
(2.2) In addition, it seems that the model is essentially guaranteed to give an alignment along the long axis of an elliptical enclosure because the catastrophe rate is governed by the contact angle between the MT and boundary.
We agree with the reviewer that the angledependence of the MT catastrophe rate at the boundary ensures alignment with the long axis of an elliptical enclosure. However, it does not explain how well individual MTs are aligned with each other  the full MT angle distribution  or how it depends on the cell geometry in nonelliptical enclosures. We have previously shown that cell geometry guides MT selforganization in epithelia in vivo (Gomez et al., 2016), and hence made understanding the nature of the full MT angle distribution the focus of this paper. We highlighted this point subsection “An analytical model shows that microtubule selforganization depends on the cell geometry and minusend distribution”. We found the analytical form of the MT angle distribution, how it depends to leading order on cell geometry and weakly depends on the MT dynamic instability rates. This analytical result further predicted that the density of the MT minusends and the average cellshape are required to obtain a quantitative agreement with the in vivo experimental results. We would add here that the alignment with the cell major axis is indeed guaranteed, and it contributes to the robustness phenomenon. This is one of the main points of our work, highlighting the independence of the average MT organization of most of the details, such as the dynamics of individual MTs.
Finally, we would like to mention that after receiving the reviews, we appreciate that the title could be misleading. We renamed the paper so that the title reflects the findings more precisely as "Robustness of the microtubule network selforganization in epithelia”.
(2.3) The particular form of that rate needs to be explained and justified better.
This form of the rate was originally inspired by the wellestablished induction of catastrophe when a MT grows against a barrier (Janson et al., 2003). The role of cell borders as barriers in epithelial cells is supported by an observation that MTs buckle at the cell cortex (Singh et al., 2018). Finally, in our previous work, we discovered experimentally that the MT catastrophe rate at the cell boundary in vivo is angledependent (see Figure 5 in Gomez et al., 2016). Altogether, we are confident that the use of this particular form of MTborder interaction is welljustified and closely reflect this interaction in vivo. We thank the reviewers for this comment and have included additional justification in subsection “Stochastic simulations demonstrate robustness of microtubule selforganization for a wide range of parameter values”.
(2.4) The authors state that the simulation results are insensitive to the alphas and betas, but what about theta_{C}?
Since MTs are stiff filaments, we only tested the cases of theta_{c} between 0 to 40 degrees (at 0, 10, 20, 30, and 40 degrees) to remain within the biologically relevant regime (Dixit and Cyr, 2004). We reported our results in Figure 2, where we chose to only plot them for the smallest nonzero and the largest values of theta_{c} (10 and 40 degrees) because there was no significant difference in the MT angle distribution curves. The case of theta_{c}=0 is the case of noninteracting MTs, for which the MT angle distribution is given by the analytical model presented later in the paper. We clarified this in subsection “Stochastic simulations demonstrate robustness of microtubule selforganization for a wide range of parameter values”. We also included a new Figure 3, where we show how local microtubule density depends on theta_{c}, in contrast to the statistics captured by the MT angle distribution, which is insensitive to the strength of MTMT interaction.
(2.5) What about proteins that lie on the cell surface and potentially affect MT kinetics in a spatial manner?
Here we interpret “the cell surface” in the question to refer to the apical membrane. Although we do not know which particular proteins the reviewers refer to, we tested effects of shotstop loss and did not see any changes of average MT alignment although the local organization appears altered – as expected, MTs are more likely to run along cell boundaries. We included a description of this finding in subsection “in vivo manipulations of microtubule dynamics and stability alter microtubule density but not alignment”. The discovery of other proteins “that lie on the cell surface and potentially affect MT kinetics in a spatial manner” is a subject of future studies and is outside of the scope of this paper.
At the same time, we have recently discovered that molecules localized in the plane of MTs (e.g.actin cables) alter the MT selforganization (Płochocka et al., 2019). The main effect of actin cables was to shift the mean of the MT angle distribution towards their direction. We included a more extensive discussion about the effects of spatially distributed regulators and crosslinkers on the MT selforganization and robustness in subsection “The importance of microtubule interactions”.
3) The authors analyze the MT organization only using the angular distribution across the whole cell. But in the cell, and in the more detailed simulation model, they could look at the local nematic order, and plot its spatial distribution. Surely this will reveal more details, which the detailed simulation should recover better than the "hairyball" model. Similarly, the spatial distribution of the zipped MTs, of the locations of catastrophes etc., as function of eccentricity.
We appreciate the reviewer’s question. Our main reason for focussing on the averaged qualifiers of the MT network was to enable comparison of our computational and analytical results with in vivo data. Since it is noisy, we could only perform comparisons of the averaged quantities. At first, we attempted to use the nematic order parameter S_{2} to quantify the degree of alignment in the detailed simulations, but the density of MTs in the simulations was too low to produce promising results. This was due to the fact that in our system the rescue rate at the MT minusends on the cell boundary is low, since it mimics the situation in vivo, where the rescue rate is the same on the boundary as in the cell interior. We, therefore, returned to working with the averaged MT angle distribution to quantify the MT alignment. We have included this clarification the Introduction.
As per reviewer’s request, we have included additional results of the stochastic simulations. In particular, we have included how local MT density is affected by MT zipping in a new Figure 3. We showed how MT density changes for a cell of a given eccentricity with increasing theta_{c} for a fixed value of p_{cat}, since we found that the degree of MT zipping increases with increasing theta_{c} and is not significantly altered by changing p_{cat} in the range used in our simulations. We have added the description of these results in subsection “Stochastic simulations demonstrate robustness of microtubule selforganization for a wide range of parameter values”. Our experimental results further show that in vivo zipping varies (Figure 4), and that loss of the crosslinker shot seems to lead to more zipping near cell boundaries consistently with a reduced theta_{c}. We have added a discussion of these findings in subsection “Robustness on the tissue scale”, and in subsection “The importance of microtubule interactions”.
Finally, we did not include the spatial map of averaged zipping in the detailed simulations since it was not very insightful  zipping increases with theta_{c}, and on average more zipping happens closer to the cell center. We also did not include the map of catastrophes since they mostly occur at the cell boundary.
4) The cells are approximated as ellipses, similar to the "spherical cow" simplification. However, cells are usually treated as polygons: do the multicell Yjunctions, play a role in the MT organization?
We agree with the reviewers that the cells are polygons. As the analytical MT angle distribution in Equation 11 shows, the timeaveraged MT angle distributions are robust but differ between cells, as they depend on the geometry of the cell boundary. When we computed the MT angle distributions in cells with corners, the MT angle distribution had peaks corresponding to the corners. This is due to the fact that to leading order, the MT angle distribution is proportional to the square of the cell crossection integrated over the cell boundary. However, in vivo, this effect averages out for the tissueaveraged MT angle distribution. We hypothesize that the Yjunction geometry matters for local processes that rely on the MT dynamics, particularly on the MT plus ends, and indeed we observed local irregularities in their organization. However, this was not the focus of this manuscript. We elaborated on this point in subsection “The analytical model accurately predicts microtubuleself organization given both the experimental cell shape and distribution of microtubule minusends”, and in subsection “The importance of microtubule interactions”.
5) The assumption in the "hairyball" model that the MTs remain pointing to the ellipse center is not clear and not explained. is this an effective way to introduce the MTMT interactions that are missing from this model?
We thank the reviewers for the comment. We would like to emphasize that the hairyball model is a 0th order conceptual model that cannot include MTMT interactions. The primary purpose of considering this model was to address the question we were often asked: “Would stretching the cell give you the experimental result?”. The answer was surprising  it would for eccentricity up to 0.95. However, this model's assumptions are not physical: the MTs are not dynamic, and their average directions pointing to the center of the cell after the cell stretching is a byproduct of stretching both the cell and the MT angle distribution in the same manner. We now clarify this in subsection “A conceptual geometric model accurately predicts in vivo alignment”.
It is an interesting suggestion that perhaps “MTs remain pointing to the ellipse center” is a “way to introduce the MTMT interactions”. However, in the stochastic simulations, the MT angle distribution does not depend on the parameterization of MTMT interactions; therefore, it is probably not why the distribution is a good fit for the experimental data.
The main benefit of including this conceptual model in the paper is that it provides the cell biology community with a simple formula to fit the MT angle distribution in epithelial cells with only oneparameter  the “effective aspect ratio”. This fit works much better than other widely used distributions such as von Mises distribution. Additionally, the fact the hairyball distribution does not depend on the MT dynamic instability parameters supports the robustness conclusion, further highlighting that geometry is one of the most significant determining factors of the MT selforganization.
Reading the reviews, we realize that perhaps it would improve the paper if we emphasized that the model is conceptual in the subsection title. Hence, we changed the corresponding subsection's title to “A conceptual geometric model accurately predicts in vivo alignment” and highlighted the model’s usefulness to the community of biologists both in subsection “A conceptual geometric model accurately predicts in vivo alignment”and in the Discussion.
(6.1) In Figure 6 the multiple lines for the cellshape/ellipse are difficult to discern. We suggest to plot them separately.
As we state in the manuscript, it is not known why the hairyball model fits the experiments up to ecc=0.95. We tried to find if the resulting distribution asymptotically matches the analytical formula or has underlying physical significance, and we did not find such connections.
As stated in the manuscript, the results of the hairyball model do not include nonuniform nucleation of MTs along the membrane. However, the nonuniform minusend density could be easily included as a multiplicative factor. We tested it, and it did not change the agreement with the experimental data. We noted this in the Discussion. That said, we would prefer not to include in detail the case of nonuniform nucleation since this geometric model is already purely conceptual and excludes the MT dynamics.
We would also like to add that we tried other conceptual models but did not succeed. In particular, we tested if the lengthaveraged angle distribution of random “sticks” placed inside an ellipse explains the in vivo MT angle distribution. This tests the hypothesis that it is possible to have the experimental MT angle distribution in the absence of internal MT alignment and only due to elongated cell geometry. This did not work, and hence there is something about the stretching deformation used in the hairyball distribution that captures the fact that the MTs “sense” the cell geometry due to the fact that they are dynamic and anchored at the cell boundary. We decided not to include other conceptual models in the manuscript as they did not produce meaningful results.
(6.2) Except for the highest eccentricity, the simplest "hairyball" model fits best the experiments. Why? and does it include the nonuniform nucleation of MTs along the membrane?
We altered this figure as suggested, which is now the new Figure 8.
7) We are concerned that the case for robustness is overstated. The authors should provide justification for the claims that their results (given the caveats above) would be generalizable to other kinds of MT array geometries (as claimed in the discussion). And it is important to emphasize what is being analyzed. The quantities of interest are the AVERAGED distributions, which by definition are taken on time scales long compared to the underlying microscopic dynamics.
As a general principle it is clear that once there is a separation of time scales
then the slow degrees of freedom dominate, and if reorientations occur at the boundary the law of reflection used will dictate the results. That there should be a strong dependence on the cell geometry (i.e. eccentricity) was already implicated in the work of Khuc Trong et al., (where it was MTs emanating from the boundary) and would follow from the dominance of the wall collisions.
In that sense, the present results are reminscent of much older work on wave pattern formation in the Faraday instability Gluckman et al., (1993) in which time averaging of the chaotic patterns revealed mean values that reflected the geometry of confinement. The work also reminds us of that of Dumais et al., on the relationship between cell geometry and cell division planes in plant cells, where microtubule orientations also play a role (Besson and Dumais, (2011) ). But it seems to us that a claim of generality must be made after having demonstrated that the issues neglected above do not change the results.
We agree with the reviewer that “the quantities of interest are the AVERAGED distributions” and reemphasized it in Introduction; in the text in subsection “Stochastic simulations demonstrate robustness of microtubule selforganization for a wide range of parameter value”, subsection “in vivo manipulations of microtubule dynamics and stability alter microtubule density but not alignment”, subsection “The analytical model accurately predicts microtubuleself organization given both the experimental cell shape and distribution of microtubule minusends”; and in Discussion.
We also appreciate the discussion of similarities between the results in different fields where the average of fast smallscale dynamics strongly depends on the slow degrees of freedom, which is the cell geometry in epithelia. We expect that since the systems discussed in subsection “Our system is a particular but generalizable scenario” are similar to the one we have considered, the robustness result holds, because in all of them the MTs are dynamic and anchored, while the separation of timescales is also valid. We clarified this further in the manuscript in subsection “Our system is a particular but generalizable scenario”.
Since it is impossible to test the results in all the possible MT array geometries, the authors meant the mentioned conclusion of the paper to be a suggestion (not a claim). However, we now support the findings of this paper with experiments in an additional tissue – pupal wings – subsection “in vivo manipulations of microtubule dynamics and stability alter microtubule density but not alignment”, and new Figure 5. Therefore, while we tuned down this statement, we are confident about the generality of our result.
https://doi.org/10.7554/eLife.59529.sa2Article and author information
Author details
Funding
Engineering and Physical Sciences Research Council (EP/L016508/01)
 Aleksandra Z Płochocka
Royal Society of Edinburgh (Royal Society of Edinburgh and the Scottish Government Personal Research Fellowship)
 Lyubov Chumakova
Scottish Government (Royal Society of Edinburgh and the Scottish Government Personal Research Fellowship)
 Lyubov Chumakova
Biotechnology and Biological Sciences Research Council (BB/P007503/1)
 Natalia A Bulgakova
Leverhulme Trust (RPG2017249)
 Natalia A Bulgakova
 Lyubov Chumakova
Scottish Funding Council
 Aleksandra Z Plochocka
HeriotWatt University
 Aleksandra Z Plochocka
University of Edinburgh
 Aleksandra Z Plochocka
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank the technical staff of the Wolfson Light Microscopy Facility (LMF) and the Drosophila Facility at the University of Sheffield for the support with in vivo experiments and Prof David Strutt, University of Sheffield, for critical reading of the manuscript. This research was supported by The Maxwell Institute Graduate School in Analysis and its Applications, a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council grant EP/L016508/01, the Scottish Funding Council, HeriotWatt University and the University of Edinburgh (AZP); BBSRC BB/P007503/1 (NAB); Royal Society of Edinburgh and the Scottish Government Personal Research Fellowship (LC). Additional support was provided by the Leverhulme Trust grant RPG2017–249 (to LC and NAB).
Senior Editor
 Anna Akhmanova, Utrecht University, Netherlands
Reviewing Editor
 Raymond E Goldstein, University of Cambridge, United Kingdom
Publication history
 Received: June 1, 2020
 Accepted: January 26, 2021
 Accepted Manuscript published: February 1, 2021 (version 1)
 Version of Record published: March 1, 2021 (version 2)
Copyright
© 2021, Płochocka et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,117
 Page views

 220
 Downloads

 0
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.