On the origin of universal cell shape variability in confluent epithelial monolayers
Abstract
Cell shape is fundamental in biology. The average cell shape can influence crucial biological functions, such as cell fate and division orientation. But cell-to-cell shape variability is often regarded as noise. In contrast, recent works reveal that shape variability in diverse epithelial monolayers follows a nearly universal distribution. However, the origin and implications of this universality remain unclear. Here, assuming contractility and adhesion are crucial for cell shape, characterized via aspect ratio (r), we develop a mean-field analytical theory for shape variability. We find that all the system-specific details combine into a single parameter α that governs the probability distribution function (PDF) of r; this leads to a universal relation between the standard deviation and the average of r. The PDF for the scaled r is not strictly but nearly universal. In addition, we obtain the scaled area distribution, described by the parameter μ. Information of α and μ together can distinguish the effects of changing physical conditions, such as maturation, on different system properties. We have verified the theory via simulations of two distinct models of epithelial monolayers and with existing experiments on diverse systems. We demonstrate that in a confluent monolayer, average shape determines both the shape variability and dynamics. Our results imply that cell shape distribution is inevitable, where a single parameter describes both statics and dynamics and provides a framework to analyze and compare diverse epithelial systems. In contrast to existing theories, our work shows that the universal properties are consequences of a mathematical property and should be valid in general, even in the fluid regime.
Editor's evaluation
In this important study, the authors unveil the reason for nearly universal shape fluctuations that have been reported earlier by theoretically analysing the energy of a confluent epithelial tissue. The comparison of their analytic results with simulations and experimental data is compelling, only the justification of the cell area distribution is somewhat incomplete. The manuscript is relevant for all people with an interest in tissue structure and dynamics.
https://doi.org/10.7554/eLife.76406.sa0Introduction
D’Arcy Thompson argued, in his book On Growth and Form, physical principles could explain tissue packing and cell shape (Thompson, 1917). Shape formation of tissues and organs during embryogenesis is a long-standing, fascinating problem of developmental biology. Since cells are the functional units of a tissue, shapes in the organs must originate at the cellular level (Paluch and Heisenberg, 2009; Wyatt et al., 2015; Pérez-González et al., 2021; Hannezo et al., 2013). Cell shapes are vital in both health and disease. As cancer progresses (Sailem and Bakal, 2017; Park et al., 2016), as asthma advances (Park et al., 2016; Park et al., 2015; Veerati et al., 2020; Atia et al., 2018), as wounds heal (Nnetu et al., 2012; Poujade et al., 2007), as an embryo develops (Farhadifar et al., 2007; Atia et al., 2018), cells progressively change their shape. Besides, cell shape may influence crucial biological functions, such as cell growth or selective programmed cell death (apoptosis) (Chen et al., 1997), the orientation of the mitotic plane (Wyatt et al., 2015; Bosveld et al., 2016; Xiong et al., 2014), stem cell lineage (McBeath et al., 2004; Wang et al., 2011), terminal differentiation (Watt et al., 1988; Roskelley et al., 1994), and division-coupled interspersion in many mammalian epithelia (McKinley et al., 2018). Moreover, the nuclear positioning mechanism in neuroepithelia depends on cell shape variation (Yanakieva et al., 2019). Thompson regarded cell-to-cell shape variability as a biologically unimportant noise Thompson, 1917; however, it is now known that shape variability is not an exception but a fundamental property of a confluent cellular monolayer (Graner and Riveline, 2017). In a seminal work, (Atia et al., 2018) showed that cell shape variability, quantified by the aspect ratio (), follows virtually the same distribution across different epithelial systems. But, the origin of this near-universal behavior, and whether it is precisely universal, remains unclear.
Previous works have shown the similarity in dynamics between cellular monolayers and glassy systems (Angelini et al., 2011; Park et al., 2015; Malinverno et al., 2017; Garcia et al., 2015). Additionally, crucial insights into the dynamics of cellular systems are obtained via simplified model systems, treating cells as polygons (Park et al., 2015; Bi et al., 2015; Bi et al., 2016; Chiang and Marenduzzo, 2016; Sadhukhan and Nandi, 2021; Sussman et al., 2018). One of these models, the vertex model (VM), shows a rigidity transition akin to the jamming transition (Bi et al., 2015; Merkel et al., 2019). However, this transition does not exist in other confluent models, such as the Voronoi model (Sussman and Merkel, 2018) or the cellular Potts model (CPM) (Sadhukhan and Nandi, 2021). But, all three models are similar from the perspective of glass transition (Bi et al., 2016; Sussman et al., 2018; Chiang and Marenduzzo, 2016; Sadhukhan and Nandi, 2021). A point of clarification on the terminology is quintessential here. The jamming (rigidity) transition is a zero-temperature, zero-activity phenomenon. It is a genuine phase transition characterized via an order parameter. The observed shape index, , the ratio of perimeter to the square root of cell area, is an order parameter of the jamming transition (Park et al., 2015; Bi et al., 2015; Merkel et al., 2019). By contrast, the glass transition refers to the extreme dynamical slowing down when the relaxation time, , becomes a certain value, usually taken as in experiments. It is not associated with any phase transition, and no static order parameter exists (Berthier and Biroli, 2011). Extensive research in the last couple of decades shows that although jamming and glass transitions can coexist, they are distinct phenomena controlled by entirely different physics (Mari et al., 2009; Biroli and Garrahan, 2013; Berthier and Witten, 2009; Ikeda et al., 2012). However, these terms are often used imprecisely while describing the dynamics of biological systems. The ‘jamming transition’ is sometimes used to describe the changing system behavior from fluid-like fast to solid-like slow dynamics (Berthier et al., 2019; Atia et al., 2021); this is actually the glass transition. Their synonymous use may lead to erroneous conclusions. Since both transitions exist in this field, distinguishing them is crucial to avoid confusion (Park et al., 2015; Chiang and Marenduzzo, 2016; Sussman et al., 2018; Sadhukhan and Nandi, 2021).
Atia et al., 2018 argued that the universal behavior of aspect ratio in an epithelial monolayer is related to the jamming transition. In a jammed system, one can carry out a voronoi tesselation centering each particle and obtain the tesselated volume, , of the individual particles. is known to follow a -Gamma distribution, , where is the Gamma function and the distribution is characterised by a single parameter (Aste and Di Matteo, 2008). From this association, (Atia et al., 2018) conjectured the viability of describing the distribution of via . But, the analytical derivation of the -Gamma distribution in granular packing (Aste and Di Matteo, 2008) relies on the fact that is additive, whereas, as the authors of Atia et al., 2018 rightly point out, , is not. Thus, there “exists no rigorous basis for the applicability of the k-Gamma distribution” (Atia et al., 2018). Yet, (Atia et al., 2018) and several subsequent works Lin et al., 2018; Li et al., 2021; Kim et al., 2020; Wenzel and Voigt, 2021 have shown that the probability distribution function (PDF) of scaled aspect ratio, rs (defined in ‘Aspects of universality’), in diverse biological and model systems can be described by , and the value of is nearly the same, around 2.5, for these diverse systems. Furthermore, the standard deviation, , vs the mean aspect ratio, , follows a universal relation (Atia et al., 2018; Kim et al., 2020; Ilina et al., 2020). What is the origin of this universality? What determines the value of ? How is the PDF of related to the microscopic properties of a system? Answers to questions like these are crucial for deeper insights into the cell shape variability and unveiling the implications of the universality. However, it requires an analytical theory that is rare in this field due to the inherent complexity of the problem and the presence of many-body interactions.
The existing theoretical works on cellular shapes assume solid like property of the system: specifically, the property that deformation of a solid is elastic in nature, that is the system resumes its shape when the external force is relaxed. For example, cell shape in the context of the jamming transition within the VM has been studied in Czajkowski et al., 2018. In a recent work, Li et al demonstrated the presence of a gamma distribution of in a wide variety of systems (Li et al., 2021). The analytical framework of Li et al., 2021 related this distribution to Boltzmann-like features and the elastic nature of the system. However, the origin of the universal properties and what dictates the value of remained unclear. Furthermore, in the context of cellular systems, the solid-like nature is only applicable deep in the glassy regime or in a jammed system (Czajkowski et al., 2018). But, most biological systems are fluid-like due to activity. The universal behavior in diverse systems seems to suggest that a more general mechanism must exist. In this work, we take a different approach. We develop a mean-field analytical theory for cell shape variability without consideration of solid-like nature or rigidity, and thus, our results should be applicable even in the fluid regime. Crucially, our work reveals that the origin of the universal features is a mathematical property.
The main results of this work are as follows: (1) We find that the aspect ratio distribution is described by a single parameter, α, containing all the system-specific parameters. Having a single parameter within the theory implies that determines the distribution. We demonstrate this in Figure 4f, illustrating the predictive power of the theory. This also implies a universal relation between and (Figure 1e, Figure 3f, Figure 4d). (2) The PDF of rs is not strictly, but nearly, universal; is a direct consequence of a mathematical property. The k-Gamma distribution for rs is obtained as a rough approximation of our analytical expression. Crucially, this nearly universal distribution of rs exists for systems even in the fluid regime; we demonstrate this by comparing our analytical results with our simulations (Figure 3c) and existing experimental data (Figure 4c). (3) We also obtain the PDF for the scaled area, , and show that it is not universal, in contrast to what has been proposed elsewhere (Wilk et al., 2014; Figure 2 and Figure 4e). (4) We demonstrate that simultaneous measurements of the PDFs for and can reveal the effects of changing physical conditions, such as maturation, on the individual model parameters. We have verified our theory via simulations of two distinct models of a confluent epithelial monolayer: the discrete lattice-based CPM on square and hexagonal lattices and the continuous VM (see Appendix 2 for details). Moreover, comparisons with existing experimental data on a wide variety of epithelial systems show excellent agreements.
One remarkable aspect of our work is as follows. It is often hard to control a specific property in a biological system as a perturbation can significantly affect other proteins. However, α includes all such effects. Therefore, even in the absence of detailed knowledge of the individual changes, our theory allows the characterization of different effects in diverse systems by treating α as a control parameter. Such a characterization can illuminate the mechanistic notion if functions and shapes are related irrespective of the molecular details. We further demonstrate in our simulations that α can also be a parameter for the dynamics. Thus, the same parameter describes both statics and dynamics, governs the origin and aspects of universality, and provides a framework to analyze and compare diverse epithelial systems.
Results
Analytical theory for the shape variability
Simplified model systems, representing cells as polygons, have been remarkably successful in describing both the static and dynamic aspects of an epithelial monolayer (Honda, 1978; Honda and Eguchi, 1980; Farhadifar et al., 2007; Bi et al., 2014; Fletcher et al., 2014; Bi et al., 2015; Park et al., 2015). The energy function, , governing these models is
where is the number of cells, the first term constrains area, , to a target area, A0, determined by cell height and cell volume, with strength . Cell heights in experiments remain almost constant in an epithelial monolayer (Farhadifar et al., 2007). The second term describes cortical contractility and adhesion (Figure 1a; Bi et al., 2015; Farhadifar et al., 2007; Prost et al., 2015). It constrains the perimeter, , to the target perimeter with strength . The energy function, Equation 1, can be numerically studied (Albert and Schwarz, 2016) via different confluent models, such as the VM (Farhadifar et al., 2007; Bi et al., 2015), the Voronoi model (Bi et al., 2016; Honda, 1978), or more microscopic models such as the CPM (Graner and Glazier, 1992; Hirashima et al., 2017; Hogeweg, 2000) and the phase-field model (Nonomura, 2012; Palmieri et al., 2015; Wenzel and Voigt, 2021). For concreteness, we mostly focus our discussions below within the CPM that has been demonstrated to be more appropriate for variable cellular shapes and sizes (Bosveld et al., 2016). However, our analytical results should generally apply for a confluent system independent of the microscopic details of the models, and we have verified them via numerical simulations in both the VM and the CPM on a square and hexagonal lattice. Furthermore, we have neglected cell division and apoptosis in our simulations for the results presented in the main text; their rates are usually low in epithelial monolayers (Poujade et al., 2007; Park et al., 2015). For example, they are of the order of per hour and per day, respectively, for an MDCK monolayer (Puliafito, 2017). Nevertheless, we show in Appendix 4, that the general conclusions of the theory remain unchanged when their rates are not very large.
Our starting point is the energy function, Equation 1, describing a confluent system of cells. We assume that the probability of a specific cellular configuration is given by a Boltzmann weight at an effective temperature, (see ‘Details of the analytical calculation’ for details). Note that in an active system includes contributions from all possible activities and the equilibrium temperature. An exact interpretation of depends on the system, and several definitions of exist, for example, the ratio of correlation to response function (Nandi and Gov, 2018; Petrelli et al., 2020; Nandi and Gov, 2017), from Einstein relation (Szamel, 2014), etc. The confluent models (such as the VM or the CPM) for epithelial systems, have two main variants: depending on the presence or absence of activity in the form of self-propulsion. The second variant represents equilibrium systems; is treated at the same footing as an equilibrium temperature and provides good agreements with experiments (Glazier and Graner, 1993; Hirashima et al., 2017; Sussman et al., 2018; Fletcher et al., 2014). Thus, the Boltzmann distribution is justified, at least within our simulations (Appendix 1). Excellent agreements of our results with experiments and analyses of the experimental systems in terms of an effective temperature (Atia et al., 2018; Kim et al., 2020) also validate this description. (Li et al., 2021) also finds a Boltzmann distribution description at an effective temperature applies to a wide variety of systems. An exact analytical calculation for the distribution of is impractical; therefore, we have made several simplifying assumptions. They are either motivated by or justified in our simulations. Here, we briefly discuss the main aspects of the calculation and relegate the technical details to the Materials and Methods (‘Details of analytical calculation’) . A detailed comparison of the analytical theory with our simulations and justifications of the assumptions are shown in ‘Comparison with simulations’.
One crucial aspect of these model systems (such as the VM, the CPM, or the Voronoi model) of epithelial monolayer is the constraint of confluency, which is area fraction is unity at all times; it enters the problem via the area term in Equation 1. This constraint is an intricate mathematical problem rendering a direct analytical calculation impractical. However, we can bypass this difficulty and gain valuable insights into the distribution of . First, a thin actomyosin layer, known as cortex, mainly governs the cellular mechanical properties (Prost et al., 2015). Therefore, the perimeter term should be dominant in determining shape. Second, shape fluctuation in these models can occur only via changes in the boundary. Third, being non-dimensional can vary independently of the cell area. Thus, in the regime of our interest, when the cells are compact objects (in contrast to being fractal-like in other regimes Sadhukhan and Nandi, 2021, see Appendix 2), we expect the area term in Equation 1 to be not crucial in determining . We have tested this assumption in two different ways in our simulations. First, if the area term is not paramount, the distribution of should not depend on ; as detailed later (Figure 2b and Figure 2—figure supplement 1), this is indeed true. Second, in the regime of our interest, the distribution of of a single cell (treating the rest of the system as medium Graner and Glazier, 1992; Glazier and Graner, 1993) is nearly the same as that of a confluent system (Figure 2a). Moreover, as discussed below, the energy from the area term varies only slightly as changes. Therefore, we assume that the area constraint is satisfied and not crucial for the aspect ratio distribution.
Then the energy function, Equation 1, becomes a sum of energies coming from individual cells. Since the perimeter of a cell is independent of that of others, we can concentrate on a particular cell, , with energy
where the first term represents contractility, and the second, effective adhesion. We have ignored the constant part as it does not affect any system properties. We first develop a coarse-grained description of cell-perimeter designating it via a set of representative points, as described in ‘Details of the analytical calculation’. To calculate the aspect ratio, , we first need to obtain the two radii of gyrations, s1 and s2, around the two principal axes (see ‘Radius of gyration’ for the definitions). Then , considering without loss of generality. However, a direct calculation of s1 and s2 is intricate due to their anisotropic natures. A slightly simpler calculation is possible for , the radius of gyration around the center of mass, and we have (Davis and Denton, 2018). Therefore, we first obtain the distribution of , , and then using this, obtain . As detailed in ‘Details of the analytical calculation’, using , with being the cell area, we obtain
where the normalization constant is determined via the constraint that total probability must be unity and with being a constant.
Additionally, as detailed in ‘Distribution for area’, Equation 18 together with the constraint of confluency (Weaire, 1986; Gezer et al., 2021), give the distribution for the scaled area , where is the average of area. It is a Gamma distribution, with a single parameter μ,
Since μ is related to the constraint of confluency, it should be independent of ; our simulations show that this is indeed true (Figure 3e). Therefore, α and μ together can distinguish how the model parameters and are affected by changing conditions such as maturation.
Note that the power 3/2 of the algebraic term in Equation 3 comes from the mathematical property of closed-looped objects. That is, for closed-loop objects, the lowest non-zero eigenvalue will have degeneracy 2, leading to the exponent 3 in the algebraic part of Equation 18, as detailed in ‘Details of the analytical calculation’. As shown below (next section), Equation 3 can be roughly approximated as a -Gamma distribution for rs that has been fitted with experimental and simulation data (Atia et al., 2018; Kim et al., 2020; Li et al., 2021; Wenzel and Voigt, 2021). The value of found in these fits comes from this mathematical property. Since the perimeter of a cell must be closed-looped, this mathematical property is inevitable. On the other hand, all the system specific details are contained in the parameter . We treat it as a free parameter in the theory and obtain its value via fits with data. Thus, Equation 3 provides a remarkable description, where all the system-specific details enter through a single parameter, α; it has profound implications leading to the universal behavior as we now illustrate.
Aspects of universality
We show the PDF of the aspect ratio, , at different values of α in Figure 1b, decays faster as increases, as expected from Equation 3. The plots look remarkably similar to the experimental results shown in Atia et al., 2018; we present detailed comparisons with experiments later. (Atia et al., 2018) has demonstrated that the PDFs of the scaled aspect ratio, , where is the ensemble-averaged , across different systems follow a near-universal behavior. We now plot the PDFs of rs in Figure 1c. The PDFs almost overlap, but they are not identical. A closer look at Equation 3 shows that if goes as , we can scale out of the equation and obtain a universal scaled distribution for rs. However, as shown in Figure 1d, there is a slight deviation in with the functional form of . This tiny deviation implies that the PDF of rs is not universal. If one ignores compared to , Equation 3 becomes a -Gamma function for rs with . However, since , this cannot be a good approximation, and the observed spread of around 2.5 is natural when fitted with this function (Atia et al., 2018; Li and Ciamarra, 2018; Kim et al., 2020; Wenzel and Voigt, 2021). On the other hand, since the deviation (Figure 1d) is minute, the PDFs of rs for different systems look nearly universal (Figure 1c): the α-dependence becomes so weak for the PDFs of rs that they seem to be independent of α. This result is a strong prediction of the theory and, as we show below, is corroborated by available experimental data on diverse epithelial systems.
Although the PDFs of rs are not strictly universal, there is another aspect, vs , which is universal. We show the as a function of in Figure 1e. Since there is only one parameter, α, in Equation 3, it determines both and . The monotonic dependence of on α (Figure 1d) implies a unique relationship between them. Therefore, we can express in terms of and, in turn, as a function of . Since there is no other system-dependent parameter in this relation, it must be universal. Note that at a constant , thus, α increases as increases or decreases. Both and become smaller as α increases, and the system on the vs plot moves towards the origin (Figure 1e). From the perspective of the dynamics, the relaxation time, , of the system grows as increases (Sadhukhan and Nandi, 2021). Thus, small and large , that is, less elongated cells and slow dynamics, follow each other, and the energy function, Equation 1, controls both behaviors. Finally, we show in Figure 1f some representative PDFs, Equation 4, at different values of μ for the scaled area . The PDF of has been argued to be universal Wilk et al., 2014; our theory shows that although the PDFs of for different systems follow the same functional form, they are not identical.
Comparison with simulations
We now compare our analytical theory with simulations of two distinct confluent models: the CPM and the VM. In our simulations, we use the original energy function, Equation 1, and other simulation details are presented in the Appendix 2. Unless otherwise stated, the CPM simulations, presented in the main text, are on the square lattice; the data for hexagonal lattice CPM simulations are shown in Figure 3—figure supplement 1. We first present tests of the crucial assumption that, in the regime of our interest here, we can ignore the constraint of confluency that enters via the area term in Equation 1. We simulate, within the CPM, single cells treating the rest of the system as medium and compare the distribution of with that in a confluent system. As shown in Figure 2a, the PDFs are nearly the same. Next, we have simulated the confluent systems with varying and find that the distribution of remains almost independent of . We have obtained that characterizes the aspect ratio distribution via Equation 3 at fixed , , and but varying . As shown in Figure 2b, within both the CPM and the VM, α remains almost constant with varying (the distributions are shown in Figure 2—figure supplement 2). Finally, the scaled area, , given by Equation 4, is a sharply peaked function around . We show the distribution of within the CPM in Figure 2c: the lines are fits with Equation 4 with μ as a fitting parameter. Thus, assuming that cells satisfy the target area, we can ignore the area term in Equation 1 to obtain the distribution of .
Figure 2c shows that becomes more sharply peaked as increases; this makes sense as greater values of ensure the area constraint is more effective and the distributions become sharply peaked around the average area. Thus, the standard deviation of decreases as increases. This implies that the variation in energy from the area term in Equation 1 is even less. We have checked that it is less than 10% for about a 300% change in . Thus, treating the area part of Equation 1 as a constant is justified. Figure 2d shows that μ almost linearly increases with . We present similar results for the VM in Figure 2—figure supplement 3. Since the area term is related to the cell height that remains nearly constant and the geometric constraint of confluency, we do not expect a substantial variation in in a particular system. However, μ also varies with (see Figure 2—figure supplement 4). Thus, in contrast to what has been proposed elsewhere (Wilk et al., 2014), as changes, though the PDF of is well-described by the same function, Equation 4, the values of μ can be different. Thus, the PDFs for different systems or the same system at different levels of activity and maturation need not be identical.
We now show that our analytical theory agrees well with simulation data. To highlight that the distribution of and its associated universal properties are also valid in the fluid regime, we have mostly simulated the systems in this regime. However, some of the simulations are also in the glassy regime (with relaxation time greater than 104). Since glass transition is not associated with any thermodynamic transition, we do not expect a drastic change in the static properties, such as the distribution of . Figure 3a and b show representative plots for the comparison of the PDFs of within the CPM, and the VM, respectively, where the lines represent fits with Equation 3. Figure 3a shows data with varying , and Figure 3b shows data with changing . As discussed above, our theory predicts nearly universal behavior for the PDFs of rs (Figure 1c). We plot the simulation data for the PDFs of rs for both the models at different parameters in Figure 3c; the PDFs almost overlap, consistent with the theory (see Figure 3—figure supplement 1 for more results). An important prediction of the theory is that the parameter , which governs the behavior of the cell shape variability, is linearly proportional to both and , hence with . This prediction also agrees with our simulations (Figure 3d and Figure 3—figure supplement 2). The slopes within the CPM and the VM are different; this possibly comes from the distinctive natures of the two models, but the qualitative behaviors are the same.
Figure 3e shows α vs μ within the CPM when we vary one of the parameters, , , and , keeping the other two fixed. First, when increases, the value of μ increases, but α remains almost constant (also see Figure 2b). Next, when increases, although α linearly increases, μ remains nearly the same. Finally, both parameters linearly increase with ; since higher implies more fluctuations, decreasing helps both and to become sharply peaked (see Figure 2—figure supplement 4, Figure 3—figure supplement 2 for their specific behaviors, and results within the VM). These results show when remains constant, varying and have distinctive effects on μ and α. These results are significant from at least two aspects: First, μ comes from the constraint of confluency (see ‘Distribution for area’), which should depend only on the area and be independent of the perimeter. Thus, the -independence of μ validates the phenomenological implementation Weaire, 1986; Gezer et al., 2021 of this constraint. Second, these results can provide crucial insights regarding the model parameters. The maturation of a monolayer can affect both and . Additional junctional proteins may be employed during maturation to increase . On the other hand, different forms of activity may reduce, decreasing . Since increases linearly with , alone is not enough to determine the dominant mechanism during the maturation process. However, assuming that remains constant in a particular system, simultaneous measurements of μ and α allow distinguishing effects of changing physical conditions, such as maturation, on the individual parameters.
We next verify the universal result of the theory: vs . Figure 3f shows vs within the CPM; we plot the theoretical prediction by the dotted line for comparison. The theory predicts that the state points move towards the origin as increases (Figure 1e); this is consistent with our simulations. Since , higher should correspond to slower dynamics. To test this hypothesis, we have simulated the CPM at different , , and to obtain the relaxation time, (see Appendix 2 for details). From these control parameters, we have calculated α and then , using our theory. We show as functions of and in Figure 3g; it is clear that indeed grows as increases or decreases. To show this behavior of on the same plot as vs , we plot in Figure 3f as a function of . A monolayer fluidizes under compressive or stretching experiments, where cell shape changes, but not cell area (Krishnan et al., 2012; Park et al., 2015; Atia et al., 2018). Such perturbations make the cells more elongated, increasing ; thus, our theory rationalizes the decrease in associated with fluidization under such perturbations. Finally, we show that our mean-field result that decreases linearly with agrees with simulations (Figure 3h). Further, to test our hypothesis that our main results remain unchanged in the presence of cell division and apoptosis, when their rates are small, we have simulated the CPM, including these processes. As shown in the Appendix 4, the simulation results justify our hypothesis.
Comparison with existing experiments
Having shown that our theory agrees well with both the CPM and the VM simulations, we next confront it with the existing experimental data. We first compare the theory with data taken from Atia et al., 2018 for three different confluent cell monolayers: the MDCK cells, the asthmatic HBEC, and the non-asthmatic HBEC. We chose the PDFs at three different times from Figure 3a and c in Atia et al., 2018. We fit Equation 3 with the data to obtain and present their values in Table 1; the corresponding fits for the MDCK cells are shown in Figure 4a (see Figure 4—figure supplement 1 for the other fits). Table 1 shows that α increases with maturation. Thus, progressive maturation can be interpreted as an increase in either or or both. The PDFs for rs corresponding to the MDCK cells are shown in the inset of Figure 4a together with one set of experimental data Atia et al., 2018; note that this is not a fit, yet the theory agrees remarkably well with the data.
We next calculate as a function of using the values of α, noted in Table 1 for the three systems. They are shown in Figure 4b along with the theory prediction. With maturation, the state points move towards lower , represented by the arrow in Figure 4b. As shown in Figure 3g, larger α corresponds to a system with higher . Thus, with maturation, as the PDFs become sharply peaked, as the cells look more roundish and becomes smaller, the system becomes more sluggish. This maturation effect is the same in all the systems (Figure 4b) and agrees with the interpretation presented in Atia et al., 2018. We have also examined that the theoretical prediction of vs agrees well with the experimentally measured values shown in Atia et al., 2018.
Our theory predicts that the PDF for rs, although not strictly universal, should be almost the same for different systems (Figure 1c). This prediction is a consequence of a crucial aspect of the theory: all the system-specific details enter via a single parameter, α in Equation 3. As shown in Figure 1d, deviates slightly from the behavior . This slight deviation implies that the PDF for rs can not be strictly universal and manifests as a variation in k when the PDF is fitted with the k-Gamma function in different experiments and simulations (Atia et al., 2018; Kim et al., 2020; Li and Ciamarra, 2018; Lin et al., 2020; Wenzel and Voigt, 2021). Nevertheless, since the deviation in Figure 1d is very weak, the values of are very close to each other. Therefore, the PDFs for rs in diverse epithelial systems–in experiments, simulations, and theory–should be nearly universal. To test this prediction, we have collected existing experimental and simulation data on different systems and show the PDFs of rs in Figure 4c. The variety in our chosen set is spectacular: it consists of various cancer cell lines (Kim et al., 2020), both asthmatic and non-asthmatic HBEC cells, MDCK cells (Atia et al., 2018), Drosophila wing disk (Lin et al., 2020), simulations data on both active (Li and Ciamarra, 2018) and equilibrium versions of the VM, the active Voronoi model (Atia et al., 2018), and the CPM. Yet, the PDFs shown in Figure 4c look nearly universal and in agreement with our analytical theory.
Additionally, our theory predicts a strictly universal behavior for vs . Since this relation does not have any system-specific details, data across diverse confluent monolayers must follow a universal relationship. We have collected existing experimental data for several systems: cancerous cell lines (Kim et al., 2020), human breast cancer cells (Ilina et al., 2020), and a jammed epithelial monolayer of MDCK cells (Fujii et al., 2019). Figure 4d shows the experimental data together with our theoretical prediction; the agreement with our theory, along with the aspect of universality, is truly remarkable. As α increases, dynamics slows down, and the points on this plot move towards lower . This result is consistent with the finding that cell shapes are more elongated and variable as the dynamics become faster in different epithelial systems (Atia et al., 2018; Park et al., 2015).
We have argued that simultaneous measurements of the PDFs of cell area and distinguish the effects of maturation on the two key parameters: and . The argument relies on the negligible effect of the perimeter constraint on μ (Equation 4). We now show a comparison of our theoretical result for the PDF of with existing experiments. Figure 4e shows experimental data for four different systems (Wilk et al., 2014; Puliafito, 2017) and the corresponding fits of Equation 4. Unlike what has been proposed elsewhere that epithelial monolayers have a universal area distribution (Wilk et al., 2014), we find, in agreement with experiments, that although the functional form remain the same, the distribution can vary.
What are the implications of these universal aspects of cell shape variability and our theory? Cell shape controls several crucial biological functions such as the mitotic-orientation (Wyatt et al., 2015; Bosveld et al., 2016; Xiong et al., 2014) and cell fate (McBeath et al., 2004; Wang et al., 2011; Roskelley et al., 1994). Our theory shows that the microscopic system properties are encoded via a single parameter, . Consequently, knowledge of one of the observables, such as , contains the information of the entire statistical properties in a monolayer. We now illustrate this predictive aspect of the theory. Experimental measurement of an average property is usually less complex and more reliable. We have collected the data for from the supplementary material of Kim et al., 2020 for three different systems: 10Ca1A, 10AT, and 10 A.ErbB2, shown by the hexagrams in Figure 4d. From these average values, we obtain α, which we use to theoretically calculate the PDFs for , as shown in Figure 4f. The inset of Figure 4f shows our theoretical PDFs for rs, together with the corresponding experimental data for comparison. The excellent agreement demonstrates that cell shape variability results from the geometric constraint imposed by the energy function, Equation 1, and is not a choice but inevitable for such systems. This result, we believe, will foster analysis of diverse epithelial systems to understand the interrelation between geometric properties and biological functions within a unified framework.
Discussion
We have obtained a mean-field theory for cell shape variability through the energy function , (Equation 1; Farhadifar et al., 2007; Honda, 1978; Bi et al., 2015; Bi et al., 2016; Sadhukhan and Nandi, 2021). We have used simplifying assumptions for analytical tractability and justified them in detailed simulations of the VM and the CPM on a square and a hexagonal lattice. The geometric restriction of confluency is a strong constraint on cell area. Considering that the area constraint is satisfied and that the cell cortex, described by the perimeter term, is crucial in determining the cell shape, allowed us to ignore the area term and obtain the distribution of . We have justified this assumption in our simulations in the regime of our interest where cells are compact objects. A detailed comparison of our simplified analytical theory with simulations and experiments shows excellent agreements. Recent experiments and simulations have revealed that cell shape variability is nearly universal in confluent epithelial monolayers Atia et al., 2018; Wenzel and Voigt, 2021; Kim et al., 2020; Ilina et al., 2020; our work provides the theoretical basis for such behavior. We have shown that the universal properties are associated with a mathematical property and valid in general, even in the fluid regime; this is significant since most biological systems are in the fluid regime due to activity. Our analytical theory reveals that the microscopic system properties enter the distribution via a single parameter, α: this leads to the universal behavior for vs and a nearly universal distribution for the scaled aspect ratio rs. Thus far, the PDF of rs has been fitted with a k-Gamma function with being around 2.5. We show that a rough approximation of our expression leads to the k-Gamma function; the slight variation in comes from the fact that the PDF is not strictly but nearly universal. On the other hand, is a direct consequence of a mathematical property: the lowest degeneracy of the eigenvalues of the connectivity matrix being two for a closed-looped object, here the perimeter.
A better understanding of the connection between the theoretical parameters and different system properties is crucial to exploit the universal aspects for deeper insights. Since all of the parameters combine into α, the effects of changing physical conditions on the individual parameters are difficult to determine from the measurements of alone. describes the cortical properties, and parameterizes different biological activities, including temperature. These are effective parameters, and their direct estimation in biological systems is impractical. Our theory provides an indirect way to estimate these parameters. The cell area in a confluent system is geometrically constrained. We have used the phenomenological implementation of the constraint of confluency and obtained the PDF for (Weaire, 1986). It is a Gamma function, described by a single parameter μ. (Wilk et al., 2014) has proposed that the PDF of is universal in various epithelial monolayers. However, we show that though the area distribution follows the same function, there is a variation in μ. Our work connects μ to the microscopic model parameters of Equation 1. In particular, μ should be independent of , whereas α varies linearly with both and . This distinction, assuming remains constant, allows inferring the effects of maturation on the individual model parameters.
We have neglected cell division, growth, and apoptosis in our theory. These are nonequilibrium processes, and the assumption of effective equilibrium becomes questionable. It will be interesting to see if the nonequilibrium statistical mechanics approach of Grossman and Joanny, 2022 can be applied to obtain the distribution in the presence of these processes. However, more generally, since the rate of these processes are low, we expect the basic form of the distribution of , Equation 3 to remain the same. The algebraic part of Equation 3 comes from the geometric property that remains the same. Therefore, we hypothesized that the effect of these additional processes should enter the distribution via . We have included these processes in our simulations in Appendix 4. The results suggest that the general form of the distribution, and the prediction that all system-specific details enter via the single parameter α, remain valid when these processes are not dominant (Appendix 4).
Our work demonstrates that a single parameter, α, describes both the cell shape statistics and the dynamics. We have shown in our simulations that the relaxation time, , grows as α increases or decreases. Experiments on confluent cellular monolayer have also reported similar results (Atia et al., 2018; Park et al., 2015; Kim et al., 2020). Thus, as cells become more compact and their shape variability reduces, the monolayer becomes more sluggish. This result may have far-reaching consequences. Most experiments usually measure and analyze biological functions via (Wyatt et al., 2015; Bosveld et al., 2016). Our theory implies that such knowledge contains a wealth of information. One can obtain α from , and all other properties, such as the distribution, the standard deviation, and the dynamics, can be analytically calculated. Pattern formation in a biological system is omnipresent during development. Pattern formation is generally controlled by gradients of different morphogens, such as the wingless or Dpp in the Drosophila wing (Jaiswal et al., 2006). Such gradients can be incorporated within our model. Since biological functions are related to cell shape, our results provide a statistical way to describe the system from the perspective of cellular functions, such as division or apoptosis. Such a description, in turn, will allow us to study the pattern formation in the presence of morphogen gradients. We are currently developing, in collaboration with colleagues, such a formalism to study the mechano-chemical pattern formation in the Drosophila wing disc.
It is well-recognized that different levels of organizations in biology are mechanistically related. One fundamental open question is how molecular-level events are related to cellular machines that control the cell shape (Gilmour et al., 2017). It is a difficult question as varying a specific property in a biological system is nearly impossible. Any perturbation will have significant impacts on several other proteins. Our work shows that all these perturbations enter the cell shape variability via a single parameter. The striking predictability, demonstrated by our theory, where determines the PDF, shows that the statistical distribution of cell shape is unavoidable. How do different cells respond to this inevitable distribution? Is cellular response similar across diverse systems? How is it related to organ-level morphogenesis? Having a single parameter that describes the static and dynamic aspects at the cellular level should help compare and analyze different systems and answer these questions.
We have emphasized that the jamming and the glass transitions are distinct phenomena controlled by different physics (Mari et al., 2009; Biroli and Garrahan, 2013; Berthier and Witten, 2009; Berthier et al., 2019; Atia et al., 2021). The glass transition is not associated with any thermodynamic transition; therefore, we expect our results to remain valid even in the glassy regime. Although most of our simulations are in the fluid regime, some are also in the glassy regime. Glassiness is relevant only in the study of dynamics, and no static order parameter exists to date (Berthier and Biroli, 2011). Since our focus in this work is the static properties, we have not discussed the glass transition. By contrast, the shape index, , has been shown to be an order parameter of the jamming transition (Bi et al., 2015; Merkel et al., 2019). Our formalism can also provide the distribution of . These results will be presented elsewhere. In addition, our expression for the distribution of area is the same as that obtained for the Voronoi tessellated volume of a jammed granular system (Aste and Di Matteo, 2008). Although the Boltzmann distribution aspect at an effective temperature is similar, the bases of the two theories are quite different. This connection provides an alternative way to think about the constraint of confluency.
In conclusion, we have developed a simple mean-field theory for the aspect ratio distribution in a confluent epithelial monolayer. We show that the universal properties of a biological system, whose physics is controlled by the energy function Equation 1, come from a mathematical property. We have analytically derived the cell shape variability, characterized via . The PDF of is described by a single parameter, α. As a result, vs becomes universal, and the PDF for the scaled aspect ratio, rs, is nearly universal. A rough approximation of our analytical form for the PDF of rs leads to the k-Gamma distribution (Aste and Di Matteo, 2008) that has been fitted to date with the existing experimental data (Atia et al., 2018; Kim et al., 2020). The distribution is valid in general, even in the fluid regime. The near-universal value of is the consequence of a mathematical property; the variation results from the fact that the PDFs are not strictly universal. α can also provide information on dynamics. Having a single parameter for the statistical and dynamical aspects of an epithelial monolayer should foster a detailed comparison of diverse epithelial systems, provide insights on the relation of biological functions to shapes, and elucidate the detailed cellular responses to the inevitable shape variability.
Materials and methods
Details of the analytical calculation
Request a detailed protocolTo set up the calculation, we represent the cell perimeter in a coarse-grained description via points, where lj is the infinitesimal line-element between jth and (j+1)th points (see Appendix 5). Note that our final results are independent of this discretization of the perimeter. To calculate the aspect ratio, , we first obtain the gyration tensor related to the moment of inertia tensor, a tensor in spatial dimension two, in a coordinate system whose origin coincides with the center of mass (CoM) of the cell. Diagonalization of this tensor gives the two principal eigenvalues, and , the squared-radii of gyrations around the respective principal axes, and . However, as discussed in the main text, a direct calculation of s1 and s2 is nontrivial due to their anisotropic natures. Therefore, we calculate the distribution of the radius of gyration, , around the center of mass, and we have (Davis and Denton, 2018). The distribution of is
where is the partition function, , the volume element, , the Boltzmann constant, and , the temperature. The first -function in Equation 5 ensures that the origin coincides with the CoM of the cell; the second -function selects specific values of , giving the distribution function.
A precise mathematical description of adhesion remains unclear (Graner and Riveline, 2017; Hilgenfeldt et al., 2008; Käfer et al., 2007). The VM represents it via a line tension: , where is the length between two consecutive vertices, and , and gives the line tension. Since the degrees of freedom in the VM are the vertices, considering constant, we obtain Equation 1. The constant implies a regular cell perimeter between vertices. However, the entire cell boundary is in contact with other cells (Figure 1a), and the perimeter is often irregular in experiments (Käfer et al., 2007; Atia et al., 2018). Therefore, we take a more general description, where the tension in a line-element lj is proportional to lj with strength . Thus, the adhesion part becomes , where is a constant. At the cellular level description, this microscopic difference can be accounted for by a renormalized . Since it is unclear how to measure in experiments, we mainly restrict our discussions on and . Unless otherwise stated, we assume remains constant. However, the theory does capture the variation in , as we show within both the VM and the CPM simulations (Figure 3h).
Within our coarse-grained description, the contractile term is . Now we can evaluate the integral, Equation 5, via a Gaussian approximation. However, an exact calculation, even at this level, is complicated. For analytical tractability and to gain insights, we use the Cauchy-Schwartz inequality (Arfken et al., 2018) and write , where is a constant. Note that this inequality is exact for the CPM, where we can consider as the line element and , the number of sides comprising the perimeter (see Appendix 5). In the experiments, the perimeter is usually obtained via a similar discretization (Atia et al., 2018). As shown in Appendix 5, this inequality, with the largest , is a reasonable approximation for the perimeter term since the variation in perimeter is not too strong due to the perimeter constraint in Equation 1. Use of this inequality makes the evaluation of the integral slightly easier. In any case, it only affects the constant in the exponential that we treat as a fitting parameter within our theory. Since we are not investigating any transition here, the use of this relation is justified. Crucially, as we describe below, the parameter has a different origin that is not affected by this inequality. Then, the perimeter part of becomes , where . The contractility and adhesion act as two competing effects (Figure 1a).
Thus, in the regime of our interest, the cell perimeter with the energy given by Equation 2 governs the distribution of that we calculate via that of . In the field of polymer physics, the distribution of has been calculated (Fixman, 1962; Eichinger, 1977; Eichinger, 1980). The mathematical structure of the two problems at this stage becomes equivalent. However, note that their physics are quite different. Specifically, our assumptions will not hold in the regime of interest of the polymer physics problem. To take advantage of an established notation and provide a connection between two disparate fields, we present our calculation in the notation of Eichinger, 1977; Eichinger, 1980. As discussed earlier, we describe the cell perimeter by the vector , representing a set of points on the perimeter. Then, we have , where is the column vector, the transpose of , , the Kirchhoff’s matrix (Eichinger, 1977; Eichinger, 1980) with and , and , the two-dimensional identity tensor. Thus, we have
where . The distribution of is
where the squared radius of gyration is . Since the radius of gyration does not depend on the coordinate system, we are allowed to chose one that diagonalizes . Say the diagonal matrix is , and represents the normal coordinates in this system.
The radius of gyration can be defined as the root-mean-square distance of different parts of a system either from its center of mass or around a given axis. We have designated the former as , defined as
where is the total volume element in the system with coordinates , and is the center of mass (CoM) of the system. The other two radii of gyration can be defined around the two principal axes (since we are in spatial dimension two) passing through the CoM. We calculate these two radii of gyration by writing the inertia tensor in a coordinate system whose origin coincides with the CoM and diagonalizing the tensor. The eigenvalues and are the squared-radii of gyrations around the respective principal axes. Thus, the aspect ratio, , is obtained as , assuming . As discussed in the main text, due to the anisotropic nature of s1 and s2, a direct calculation for their distributions is more complex than that of . So, we first calculate the distribution for and then, using this result, we obtain the distribution of .
Equation 5 can be written in the normal coordinate system as.
where we have used that corresponds to the zero-eigenvalue mode of the matrix. Integrating over , we get rid of this zero-eigenvalue that gives translation. Thus,
where we have defined as the dimensional vector excluding the coordinates corresponding to the zero-eigenvalue. The normalization factor, , can be calculated exactly through the integration as
Note that the integration in the calculation of is around the boundary of the cell; to separate out the radial part, we now write the volume element in polar coordinate . Then and . Thus, we obtain from Equation 9
with is the identity matrix of rank . Carrying out the integration over , we obtain
Using the value of , we obtain
where ’s are the eigenvalues of . The integral in Equation 13 can be performed via the contour integral and the resultant solution can be written as.
where are the distinct eigenvalues of and gives the residue at the pole . As we show below, the residues will have a term and in the limit , only the smallest will contribute.
Since the cell perimeter must be closed-looped, is a tridiagonal matrix with periodicity. Therefore, the number of zero-eigenvalues must be one, and the lowest degeneracy of the non-zero eigenvalues must be two (Kulkarni et al., 1999; Witt et al., 2009; Eichinger, 1977; Eichinger, 1980). We have already integrated out the coordinate corresponding to the zero-eigenvalue. Let us designate the lowest non-zero eigenvalue as λ. The pole corresponding to λ is located at , and of order 2. Thus, we obtain the residue as
Let’s first take the derivative, with respect to β, of the numerator and write part of the residue as
Next, differentiating the denominator, we obtain the other part of the residue as
A comparison of term1 and term2, given by Equation 16 and Equation 17, respectively, shows that there is an extra factor of in the denominator of term2. Thus, in the limit term2 can be ignored compared to term1. Therefore, we obtain the distribution function for as.
where is the normalization constant that we will fix later, and . Note that the lowest eigenvalue for the -dimensional Kirchoff’s matrix is proportional to , thus .
Now, and the aspect ratio . Moreover, we have , where is the area of the cell. Since the distribution of cell area is sharply peaked (Figure 2c in the main text), as the cell division and apoptosis are slow processes, can be taken as a constant. Therefore, using the last two relations in the first, we obtain . Thus, we obtain from Equation 18, the distribution of as.
where is the normalization constant: , where , is the Gamma function, and is the Kummer’s confluent Hypergeometric function (Erdelyi et al., 1953), .
Note that the exponent 3 in the algebraic part of Equation 18 comes from the mathematical property, therefore, must be true for any system. This is the source of the value when Equation 19 is approximated as a -Gamma distribution, as discussed in the main text. For the analysis of the theory, a symbolic software, such as ‘Mathematica’ (Wolfram Research Inc, 2019) is helpful. We note here that the relation between and is quite complex, though for most practical purposes it can be taken as a straight line: .
Radius of gyration
Request a detailed protocolThere are different ways to define the radius of gyration depending on the application. Mathematically, it is the root mean square distance from the centre of mass or a given axis. We have designated the first as (Equation 7). For a two-dimensional object, there are two independent axes, therefore, there will be two radii of gyration, s1 and s2, in the second method. As detailed in the Discussion, s1 and s2 are obtained as the two eigenvalues of the gyration tensor. The significance of these two eigenvalues are as follows: the geometric properties of the system are equivalent to those of an ellipse with s1 and s2 as the major and minor axes. Considering that the mass is uniformly distributed throughout the system, the area is proportional to or . We have used this property to obtain the distribution of in the analytical derivation of Equation 19.
Distribution for area
Request a detailed protocolEquation 18 gives the distribution of the squared radius of gyration of different cells in a monolayer. Using this equation, we now derive the distribution of area, . Note that Equation 18 is valid irrespective of whether the system is confluent or not. We have argued in the main text that the constraint of confluency is not crucial to obtain the distribution of the aspect ratio. But, this argument is not valid when we are interested in the distribution of cellular area, . Since confluency is a strong geometric constraint on , we must include this constraint in the derivation of . Now, as detailed in ‘Radius of gyration’, . Therefore, Equation 18 gives.
where β is a constant related to . We now consider the constraint of confluency. It is actually a long-standing hard mathematical problem. Even for random patterns, no exact result exists yet. (Weaire, 1986) proposed a phenomenological implementation of this constraint as a polynomial function of area; this remains one of the simplest possible ways to date to deal with this constraint (Gezer et al., 2021). Keeping only one term of this polynomial for simplicity, we can write this constraint as (Weaire, 1986; Gezer et al., 2021). Using this in Equation 20, we obtain
where we have defined . From Equation 21, we obtain the average area, . Thus, we obtain the normalized distribution for the scaled area, , as
This is the well-known k-Gamma function, defined in Aste and Di Matteo, 2008, and usually denoted with the variable . This same function has been used in fitting the scaled aspect ratio, rs, data in different existing experiments and simulations. Therefore, to avoid confusion with , which is obtained fitting the rs data, we have used μ to define the distribution function for . Since μ comes from the constraint of confluency, it should be independent of . Our simulation results within both the CPM and the VM (Figure 3e, in the main text) support this hypothesis.
Appendix 1
Boltzmann Distribution
As we discussed in the main text, there are two types of models of confluent systems, depending on the presence or absence of self-propulsion. In the absence of self-propulsion, the models are in equilibrium. It has been shown in the literature, that the equilibrium variants of the models do capture the static and dynamical properties of a confluent epithelial system (Farhadifar et al., 2007; Atia et al., 2018; Park et al., 2015). For concreteness, we have considered the equilibrium variants of the models in this work. Therefore, we expect that the Boltzmann distribution should hold, that is, a configuration with energy , given by Equation 1, should have the weightage of . However, as discussed below in Appendix 2, we restrict ourselves in the low- regime, where the cells, on average, cannot satisfy the perimeter constraint. This implies a nontrivial dependence of the distribution around coming from the constraint of confluency that in turn restricts the minimum value of the perimeter; the energy cannot go below a certain value. However, we show in Appendix 1—figure 1 that we can still approximate the distribution via a Boltzmann distribution.
If the system follows the Boltzmann distribution, that is, , plots of as a function of for different should follow a master curve that is a linearly decreasing function. As shown in Appendix 1—figure 1, except at very close to , as a function of at different indeed follows a master curve, implying the existence of the Boltzmann distribution.
Appendix 2
Simulation details
We have verified our theory via simulations of two distinct models: the CPM and the VM. In both the simulations, we have used the energy function , given in Equation 1 in the main text. The CPM is a lattice-based model. The underlying lattice has some effects on the quantitative aspects of the model; for example, on the square lattice, the object with the minimum possible perimeter for a given area is a square. However, the qualitative behaviors are independent of the lattice. To ascertain that the simulation results are independent of the lattice, we have simulated the CPM on two different lattices: the square lattice and the hexagonal lattice.
For the glassy dynamics, the geometric restriction leads to two different regimes, the low-, and the large- regime. However, the large- regime characterizes quite large adhesion compared to the cortical contractility; this regime, we feel, is not relevant for the experiments. Therefore, we present most results when is not too large. Here, we briefly discuss the different models.
Cellular Potts Model (CPM)
In CPM, dynamics proceeds by stochastically updating one boundary lattice site at a time. In a Monte-Carlo simulation (MC), we accept a move with a probability, , where we have set Boltzmann constant to unity, is the change in energy going from one configuration to the other. and are the cell indices of the current cell and the target cell indices, respectively. Unit of time (1 MC time) refers to such elementary moves, where is the total number of lattice sites in the simulation.
With square lattice
We have implemented CPM with square lattice in Fortran 90 and followed a Connectivity Algorithm developed in Durand, 2016 to prevent fragmentation of cells (Sadhukhan and Nandi, 2021). We have chosen a simulation box of size with 360 total cells in the system. The average area of cells in the system is 40, and the minimum possible perimeter on a square lattice with this area is 26. Unless otherwise specified, we always start with an initial configuration where each cell is rectangular with a size . We first equilibrate the system for MC time before collecting the data of aspect ratio. The variables that characterize the system are , , , and . We have ignored cell division and apoptosis, as discussed in the main text. However, to test the effects of these processes, we have included them within this model and present a specific set of results in Appendix 4.
To obtain the relaxation time, , we have calculated the overlap function, , defined as
where is the center of mass at time of a cell with index , is a Heaviside step function
denotes averaging over initial times and the overline implies ensemble average. Relaxation time is defined as . The parameter is related to the vibrational motion of the particles (here cells) inside the cage formed by its neighbors. We have set (Sadhukhan and Nandi, 2021). We have taken 50 t0 averaging and 20 configurations for ensemble averaging.
With Hexagonal lattice
To simulate CPM on the hexagonal lattice, we have used an open-source application CompuCell3D (Swat, 2012; Zajac et al., 2003). CPM is the core of CompuCell3D. We have simulated the 2D confluent cell monolayer with periodic boundary conditions using the same energy function, Equation 1, in the main text. In this simulation, we have taken simulation box size with 361 cells. Initially, Cells are of three different types of areas and perimeters. We have started with a rectangular slab of the cell monolayer. We first equilibrate the system for 105 MC time steps before collecting the data. Fragmentation is allowed in these simulations; however, we have restricted ourselves in the low regime, where fewer cells are fragmented. While calculating the PDFs, we have eliminated the fragmented cells.
Vertex Model (VM)
Vertex models (VM) can be viewed as the continuum version of the CPM. Within VM, vertices of the polygons are the degrees of freedom, and cell edges are defined as straight lines (or lines with a constant curvature) connecting between vertices (Farhadifar et al., 2007; Fletcher et al., 2014). Within the Monte-Carlo simulatoin (Wolff et al., 2019), dynamics proceeds by stochastically updating each vertex position of all the cells by a small amount . We accept the move with a probability where we have set Boltzmann constant to unity, is the change in energy going from . Initially, we start will 1024 cells having equal area and perimeter . We have equilibrated the system for MC time before the start of collecting data.
Two different regimes
As we discussed in the main text, theoretical studies have shown the existence of two distinct regimes as a function of within these models of confluent epithelial cells (Sadhukhan and Nandi, 2021; Bi et al., 2015). The existence of these two regimes is controlled by a geometric constraint. Given a certain average cell area, the perimeter can have a minimum value. If in Equation 1 is less than this value, cells cannot achieve this target perimeter: the cells in this regime are compact. It seems that this is the regime of experimental interest. In contrast, when is large compared to this minimum value, cells can achieve this value to minimize energy . Within the CPM, the cells look like fractal objects in this regime (Sadhukhan and Nandi, 2021). This behavior is expected to be similar within both the VM and the Voronoi model. But there are technical difficulties to simulate these latter models in this regime due to numerical instabilities.
Appendix 3
Procedure to calculate aspect ratio () in our simulation
As stated in the main text, we characterize the cell shape via aspect ratio, . To calculate , we follow (Atia et al., 2018) and briefly discuss the method here. Consider one particular cell, represented by the set of points , as schematically shown in Appendix 3—figure 1. We first calculate the gyration tensor, , (proportional to the moment of inertia tensor) in a frame of reference whose origin coincides with the center of mass of the cell. is a tensor in spatial dimension two. We next diagonalize and obtain the two eigenvalues, and . Without loss of generality, we assume that , and obtain . This process can be viewed as approximating the cell with an ellipse with the major and minor axes given by s1 and s2, respectively. For the CPM on the square lattice and the VM, we have written our own codes; for the CompuCell3D (Swat, 2012; Zajac et al., 2003), we have used <Plugin Name="MomentOfInertia"/>in the.XML file and using the example of Demos/MomentOfInertia, we have calculated the lengths of semiaxes in the Steppables.py file. We have obtained the aspect ratio of a cell by taking the ratio of lengths of semi-major and semi-minor axes.
Appendix 4
Including cell division and apoptosis
We have ignored cell division and apoptosis in our theory since the rates of these processes are extremely low. However, even when these processes are present, as long as their rates are not very high, we expect the form of the distribution, Equation 3 in the main text, should remain valid. There are two distinct parts in . The algebraic part is determined by the geometric property, which is a cell perimeter is a closed-loop object. This property cannot change. Therefore, the effects of these processes can enter the distribution via alone. Thus, the main conclusions of the theory should remain valid even when division and apoptosis are present.
We have included cell division and apoptosis in our simulations to test this hypothesis. We keep the rate of division, , and apoptosis, , the same such that the average number of cells remains the same. Every kd time step, we randomly select a cell and divide it into two, with a randomly chosen division plane. To decide the division plane, we first calculate the center of mass (CoM) and then chose a point on the perimeter in a random direction. The line connecting the CoM and this point gives the division plane (a line for the two-dimensional system). The area of these two daughter cells then grows till it becomes of the order of A0 (Equation 1 in the main text). To avoid dividing a cell that has just undergone division, we impose a cut-off on the area () for selecting a cell for division. For apoptosis, every ka time step, we randomly choose a cell and assign it a target area . We also relax the constraint that does not allow fragmentation in our simulation for the cells undergoing apoptosis.
We show the PDF of at three representative values of ka above. As shown in Appendix 4—figure 1a PDFs of r at different values of can be fitted with Equation 3, the values of α for the best fits are quoted in the figure. Crucially, an important prediction of the theory that the PDF for rs becomes nearly universal remains also valid, as shown in Appendix 4—figure 1b. We find that α as a function of kd or ka quickly reaches a constant with increasing ka (that is, with decreasing rate) (Appendix 4—figure 1c).
Appendix 5
Validity of the approximation of perimeter
We have used the Cauchy-Schwartz inequality, , where is a constant. Within the CPM, the perimeter consists of elementary sides, as schematically shown in Appendix 5—figure 1a above, then . Since , . The inequality is exact in that case. The data from the simulation confirms this (Appendix 5—figure 1). In the case of the continuum models, such as the vertex model, we can use a similar discretization as schematically shown by the red line in the inset of Appendix 5—figure 1. In fact, the estimate of cell perimeter in the experiments is obtained via a similar discretization method. With such a discretization representing the continuous cell perimeter via a discrete line, we again find that the inequality is very close to equality (Appendix 5—figure 1b).
The approximation of replacing the inequality by equality now depends on the fluctuation of the total perimeter or . This fluctuation depends on the values of the model parameters, such as , , , etc. We find that for the parameter values used in this work to simulate the systems, the highest fluctuation of the perimeter, or in the value of , is less than 5%. Thus, we are justified to approximate the inequality by equality.
Data availability
We have uploaded the source files of the simulation data and the Mathematica analysis files in Dryad.
-
Dryad Digital RepositoryData from: On the origin of universal cell shape variability in a confluent epithelial monolayer.https://doi.org/10.5061/dryad.xsj3tx9h0
References
-
Modeling cell shape and dynamics on micropatternsCell Adhesion & Migration 10:516–528.https://doi.org/10.1080/19336918.2016.1148864
-
Emergence of gamma distributions in granular materials and packing modelsPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 77:021309.https://doi.org/10.1103/PhysRevE.77.021309
-
Geometric constraints during epithelial jammingNature Physics 14:613–620.https://doi.org/10.1038/s41567-018-0089-9
-
Are cell jamming and unjamming essential in tissue development?Cells & Development 168:203727.https://doi.org/10.1016/j.cdev.2021.203727
-
Glass transition of dense fluids of hard and compressible spheresPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 80:021502.https://doi.org/10.1103/PhysRevE.80.021502
-
Theoretical perspective on the glass transition and amorphous materialsReviews of Modern Physics 83:587–645.https://doi.org/10.1103/RevModPhys.83.587
-
Glassy dynamics in dense systems of active particlesJ Chem Phys 150:200901.https://doi.org/10.1063/1.5093240
-
Energy barriers and cell migration in densely packed tissuesSoft Matter 10:1885–1890.https://doi.org/10.1039/c3sm52893f
-
A density-independent rigidity transition in biological tissuesNature Physics 11:1074–1079.https://doi.org/10.1038/nphys3471
-
Motility-driven glass and jamming transitions in biological tissuesPhysical Review. X 6:021011.https://doi.org/10.1103/PhysRevX.6.021011
-
An efficient cellular potts model algorithm that forbids cell fragmentationComputer Physics Communications 208:54–63.https://doi.org/10.1016/j.cpc.2016.07.030
-
An approach to distribution functions for Gaussian moleculesMacromolecules 10:671–675.https://doi.org/10.1021/ma60057a035
-
Configuration statistics of Gaussian moleculesMacromolecules 13:1–11.https://doi.org/10.1021/ma60073a001
-
Vertex models of epithelial morphogenesisBiophysical Journal 106:2291–2304.https://doi.org/10.1016/j.bpj.2013.11.4498
-
Statistical properties of poisson-voronoi tessellation cells in bounded regionsJournal of Statistical Computation and Simulation 91:915–933.https://doi.org/10.1080/00949655.2020.1836184
-
Simulation of the differential adhesion driven rearrangement of biological cellsPhysical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 47:2128–2154.https://doi.org/10.1103/physreve.47.2128
-
Simulation of biological cell sorting using a two-dimensional extended Potts modelPhysical Review Letters 69:2013–2016.https://doi.org/10.1103/PhysRevLett.69.2013
-
Instabilities and geometry of growing tissuesPhysical Review Letters 129:048102.https://doi.org/10.1103/PhysRevLett.129.048102
-
Cellular potts modeling of complex multicellular behaviors in tissue morphogenesisDevelopment, Growth and Differentiation 59:12358.https://doi.org/10.1111/dgd.12358
-
Evolving mechanisms of morphogenesis: on the interplay between differential adhesion and cell differentiationJournal of Theoretical Biology 203:317–333.https://doi.org/10.1006/jtbi.2000.1087
-
Description of cellular patterns by Dirichlet domains: the two-dimensional caseJournal of Theoretical Biology 72:523–543.https://doi.org/10.1016/0022-5193(78)90315-6
-
How much does the cell boundary contract in a monolayered cell sheet?Journal of Theoretical Biology 84:575–588.https://doi.org/10.1016/s0022-5193(80)80021-x
-
Unified study of glass and jamming rheology in soft particle systemsPhysical Review Letters 109:018301.https://doi.org/10.1103/PhysRevLett.109.018301
-
Unjamming and collective migration in MCF10A breast cancer cell linesBiochemical and Biophysical Research Communications 521:706–715.https://doi.org/10.1016/j.bbrc.2019.10.188
-
Fluidization, resolidification, and reorientation of the endothelial cell in response to slow tidal stretchesAmerican Journal of Physiology. Cell Physiology 303:C368–C375.https://doi.org/10.1152/ajpcell.00074.2012
-
Eigenvalues of tridiagonal pseudo-toeplitz matricesLinear Algebra and Its Applications 297:63–80.https://doi.org/10.1016/S0024-3795(99)00114-7
-
Role of cell deformability in the two-dimensional melting of biological tissuesPhysical Review Materials 2:045602.https://doi.org/10.1103/PhysRevMaterials.2.045602
-
E2 and gamma distributions in polygonal networksPhysical Review Research 3:L042001.https://doi.org/10.1103/physrevresearch.3.l042001
-
Dynamic migration modes of collective cellsBiophysical Journal 115:1826–1835.https://doi.org/10.1016/j.bpj.2018.09.010
-
Universal statistical laws for the velocities of collective migrating cellsAdvanced Biosystems 4:e2000065.https://doi.org/10.1002/adbi.202000065
-
Endocytic reawakening of motility in jammed epitheliaNature Materials 16:587–596.https://doi.org/10.1038/nmat4848
-
Jamming versus glass transitionsPhysical Review Letters 103:025701.https://doi.org/10.1103/PhysRevLett.103.025701
-
Effective temperature of active fluids and sheared soft glassy materialsThe European Physical Journal. E, Soft Matter 41:117.https://doi.org/10.1140/epje/i2018-11731-7
-
Biology and physics of cell shape changes in developmentCurrent Biology 19:R790–R799.https://doi.org/10.1016/j.cub.2009.07.029
-
Unjamming and cell shape in the asthmatic airway epitheliumNature Materials 14:1040–1048.https://doi.org/10.1038/nmat4357
-
Collective migration and cell jamming in asthma, cancer and developmentJournal of Cell Science 129:3375–3383.https://doi.org/10.1242/jcs.187922
-
Cell-size distribution in epithelial tissue formation and homeostasisJournal of the Royal Society Interface 14:20170032.https://doi.org/10.1098/rsif.2017.0032
-
BookGlazierIn: Swat MH, editors. Computational Methods in Cell Biology. Academic Press. pp. 325–366.
-
Airway mechanical compression: its role in asthma pathogenesis and progressionEuropean Respiratory Review 29:190123.https://doi.org/10.1183/16000617.0123-2019
-
On the distribution of cell areas in a voronoi networkLetters Section 53:L101–L105.https://doi.org/10.1080/13642818608240647
-
Universal area distributions in the monolayers of confluent mammalian cellsPhysical Review Letters 112:138104.https://doi.org/10.1103/PhysRevLett.112.138104
-
SoftwareMathematica, version version 12.0Wolfram research inc.
-
Cell and tissue morphology determine actin-dependent nuclear migration mechanisms in neuroepitheliaThe Journal of Cell Biology 218:3272–3289.https://doi.org/10.1083/jcb.201901077
-
Simulating convergent extension by way of anisotropic differential adhesionJournal of Theoretical Biology 222:247–259.https://doi.org/10.1016/s0022-5193(03)00033-x
Article and author information
Author details
Funding
Department of Atomic Energy, Government of India (RTI 4007)
- Souvik Sadhukhan
- Saroj Kumar Nandi
Science and Engineering Research Board (SRG/2021/002014)
- Saroj Kumar Nandi
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Satyabrata Bandyopadhyay, Mustansir Barma, Tamal Das, Manish Jaiswal, Sritama Datta, Basil Thurakkal. and Sumedha for many enlightening discussions and Sam A Safran for comments. We acknowledge the support of the Department of Atomic Energy, Government of India, under Project Identification No. RTI 4007 and the computational facility of TIFR Hyderabad. We also thank SERB for grant via SRG/2021/002014.
Copyright
© 2022, Sadhukhan and Nandi
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
-
- 2,387
- views
-
- 322
- downloads
-
- 11
- citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
-
- Developmental Biology
Neuronal stem cells generate a limited and consistent number of neuronal progenies, each possessing distinct morphologies and functions, which are crucial for optimal brain function. Our study focused on a neuroblast (NB) lineage in Drosophila known as Lin A/15, which generates motoneurons (MNs) and glia. Intriguingly, Lin A/15 NB dedicates 40% of its time to producing immature MNs (iMNs) that are subsequently eliminated through apoptosis. Two RNA-binding proteins, Imp and Syp, play crucial roles in this process. Imp+ MNs survive, while Imp−, Syp+ MNs undergo apoptosis. Genetic experiments show that Imp promotes survival, whereas Syp promotes cell death in iMNs. Late-born MNs, which fail to express a functional code of transcription factors (mTFs) that control their morphological fate, are subject to elimination. Manipulating the expression of Imp and Syp in Lin A/15 NB and progeny leads to a shift of TF code in late-born MNs toward that of early-born MNs, and their survival. Additionally, introducing the TF code of early-born MNs into late-born MNs also promoted their survival. These findings demonstrate that the differential expression of Imp and Syp in iMNs links precise neuronal generation and distinct identities through the regulation of mTFs. Both Imp and Syp are conserved in vertebrates, suggesting that they play a fundamental role in precise neurogenesis across species.
-
- Computational and Systems Biology
- Developmental Biology
The initially homogeneous epithelium of the early Drosophila embryo differentiates into regional subpopulations with different behaviours and physical properties that are needed for morphogenesis. The factors at top of the genetic hierarchy that control these behaviours are known, but many of their targets are not. To understand how proteins work together to mediate differential cellular activities, we studied in an unbiased manner the proteomes and phosphoproteomes of the three main cell populations along the dorso-ventral axis during gastrulation using mutant embryos that represent the different populations. We detected 6111 protein groups and 6259 phosphosites of which 3398 and 3433 were differentially regulated, respectively. The changes in phosphosite abundance did not correlate with changes in host protein abundance, showing phosphorylation to be a regulatory step during gastrulation. Hierarchical clustering of protein groups and phosphosites identified clusters that contain known fate determinants such as Doc1, Sog, Snail, and Twist. The recovery of the appropriate known marker proteins in each of the different mutants we used validated the approach, but also revealed that two mutations that both interfere with the dorsal fate pathway, Toll10B and serpin27aex do this in very different manners. Diffused network analyses within each cluster point to microtubule components as one of the main groups of regulated proteins. Functional studies on the role of microtubules provide the proof of principle that microtubules have different functions in different domains along the DV axis of the embryo.