On the origin of universal cell shape variability in confluent epithelial monolayers

  1. Souvik Sadhukhan
  2. Saroj Kumar Nandi  Is a corresponding author
  1. Tata Institute of Fundamental Research, India

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.sa0

Introduction

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 (r), 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, q, 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 10-100s 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, x, of the individual particles. x is known to follow a k-Gamma distribution, P(x,k)=[kk/Γ(k)]xk1exp[kx], where Γ is the Gamma function and the distribution is characterised by a single parameter k (Aste and Di Matteo, 2008). From this association, (Atia et al., 2018) conjectured the viability of describing the distribution of r via P(x,k). But, the analytical derivation of the k-Gamma distribution in granular packing (Aste and Di Matteo, 2008) relies on the fact that x is additive, whereas, as the authors of Atia et al., 2018 rightly point out, r, 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 P(rs,k), and the value of k is nearly the same, around 2.5, for these diverse systems. Furthermore, the standard deviation, sd, vs the mean aspect ratio, r¯, 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 k2.5? How is the PDF of r 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 r 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 k2.5 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 r¯ determines the distribution. We demonstrate this in Figure 4f, illustrating the predictive power of the theory. This also implies a universal relation between sd and r¯ (Figure 1e, Figure 3f, Figure 4d). (2) The PDF of rs is not strictly, but nearly, universal; k2.5 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, a, 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 a and r 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.

Theoretical results for cell shape variability.

(a) Schematic illustration of a confluent model of an epithelial monolayer. Cortical contractility and adhesion impose competing forces. (b) PDF of aspect ratio, r, is governed by the parameter α. (c) PDF of the scaled aspect ratio, rs=(r-1)/(r¯-1) is nearly universal. (d) The PDF of rs can be universal if r¯+1/r¯ is proportional to 1/α (Equation 3), but there is a slight deviation showing that it is not strictly universal. Inset: Same as in the main figure, but in log-log scale. Lines are fits with the 1/α form. (e) Standard deviation (sd) vs r¯ follows a universal relation; the state points move towards the origin as α increases. (f) The PDF of a follows Gamma distribution with a single parameter, μ, Equation 4.

Figure 2 with 4 supplements see all
Tests of theoretical assumptions.

(a) PDFs of aspect ratio, r, within the CPM for a single cell in a medium are nearly identical with those for a confluent system when the parameters are the same: λA=1.0, λP=0.5, P0=38 and A0=90. Symbols and lines are data for the confluent system and single cell in a medium, respectively. (b) α, giving the distribution of r remains independent of λA. For the CPM, λP=0.5, and T=25.0; for the VM T=2.5×10-3, λP=0.02, A0=1, and P0=3.7. (c) PDF of the scaled area, a, within the CPM at different λA, but fixed T=12.5 and λP=0.5. The lines are fits with Equation 4 with the values of μ as shown. (d) μ almost linearly increases with λA. For the CPM simulations in (b) and (c), we have A0=40 and P0=26. The PDFs are calculated from at least 4×104 independent configurations.

Table 1
Values of α from fits of Equation 3 with the PDFs of r.

The data for the three systems are taken as function of maturation time, in units of h (hours) and d (days). The experimental data are taken from Atia et al., 2018.

Cell typeTimeValue of α
32 h7.535
MDCK47 h10.998
63 h12.496
6 d3.060
Asthmatic HBEC14 d4.472
20 d5.389
6 d4.798
Non-asthmatic HBEC14 d7.358
20 d9.602

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

(1) H=i=1N[λA(AiA0)2+λP(PiP0)2],

where N is the number of cells, the first term constrains area, Ai, to a target area, A0, determined by cell height and cell volume, with strength λA. 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, Pi, to the target perimeter P0 with strength λP. 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 10-2 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, T (see ‘Details of the analytical calculation’ for details). Note that T in an active system includes contributions from all possible activities and the equilibrium temperature. An exact interpretation of T depends on the system, and several definitions of T 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; T 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 r 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 r. 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, r 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 r. We have tested this assumption in two different ways in our simulations. First, if the area term is not paramount, the distribution of r should not depend on λA; 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 r 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 λA 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, i, with energy

(2) HP=λPPi22λPP0Pi,

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, r, 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 r=s1/s2, considering s1>s2 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 s, the radius of gyration around the center of mass, and we have s2=s12+s22 (Davis and Denton, 2018). Therefore, we first obtain the distribution of s2, P(s2), and then using this, obtain P(r). As detailed in ‘Details of the analytical calculation’, using s2=A(r+1/r), with A being the cell area, we obtain

(3) P(r)=1N(r+1r)3/2(11r2)eα(r+1r),

where the normalization constant N is determined via the constraint that total probability must be unity and αλP(1-KP0)/T with K 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 a=A/A¯, where A¯ is the average of area. It is a Gamma distribution, with a single parameter μ,

(4) P(a)=μμΓ(μ)aμ1exp[μa].

Since μ is related to the constraint of confluency, it should be independent of λP; our simulations show that this is indeed true (Figure 3e). Therefore, α and μ together can distinguish how the model parameters λP and T are affected by changing conditions such as maturation.

Figure 3 with 2 supplements see all
Comparison with simulations.

(a) PDF of r at different λP within the CPM, with λA=1.0, and T=25.0. Symbols are simulation data, and lines represent the fits with Equation 3. (b) PDF of r within the VM with varying T, λP=0.02, λA=0.5, P0=3.7. Lines are fits with Equation 3 with α as shown. (c) PDFs for the rs, for the same data as in (a) and (b) for the CPM and the VM, respectively, show a virtually universal behavior. (d) Theory predicts α linearly varies with λP/T, simulation data within both the models agree with this prediction. λA=1.0 for the CPM; P0=3.7 and λA=0.5 for the VM simulations. Dotted lines are fits with a linear form. (e) α vs μ within the CPM when we have varied only one of the three variables: T, λA, and λP does not depend on λA, and μ does not depend on λP. (f) Our theory predicts a universal behavior for sd vs r¯, symbols are the CPM data at different parameter values, and the line is our theory (not a fit). We have also plotted the scaled relaxation time, 2.5/lnτ0.25, to show them on the same figure. τ increases as r¯ decreases. (g) τ as functions of α (lower axis) and r¯ (upper axis). (h) Theory predicts α linearly decreases with P0, this agrees with simulations (symbols); parameters in the CPM: λA=1.0, λP=0.5, T=16.67; the VM: λA=0.5, λP=0.02, T=0.0025. [CPM simulations here are on square lattice, A0=40, P0=26 in (a–f); PDFs are calculated over at least 4×104 independent configurations].

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 k-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 k2.5 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, P(r), at different values of α in Figure 1b, P(r) 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, rs=(r-1)/(r¯-1), where r¯ is the ensemble-averaged r, 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 r¯+1/r¯ goes as 1/α, 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 r¯+1/r¯ with the functional form of 1/α. This tiny deviation implies that the PDF of rs is not universal. If one ignores 1/r compared to r, Equation 3 becomes a k-Gamma function for rs with k=2.5. However, since rO(1), this cannot be a good approximation, and the observed spread of k 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, sd vs r¯, which is universal. We show the sd=r2¯r¯2 as a function of r¯ in Figure 1e. Since there is only one parameter, α, in Equation 3, it determines both sd and r¯. The monotonic dependence of r¯ on α (Figure 1d) implies a unique relationship between them. Therefore, we can express α in terms of r¯ and, in turn, sd as a function of r¯. Since there is no other system-dependent parameter in this relation, it must be universal. Note that αλP/T at a constant P0, thus, α increases as λP increases or T decreases. Both sd and r¯ become smaller as α increases, and the system on the sd vs r¯ 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 r¯ 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 a=A/A¯. The PDF of a has been argued to be universal Wilk et al., 2014; our theory shows that although the PDFs of a 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 r 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 λA and find that the distribution of r remains almost independent of λA. We have obtained α that characterizes the aspect ratio distribution via Equation 3 at fixed λP, P0, and T but varying λA. As shown in Figure 2b, within both the CPM and the VM, α remains almost constant with varying λA (the distributions are shown in Figure 2—figure supplement 2). Finally, the scaled area, a, given by Equation 4, is a sharply peaked function around a=1. We show the distribution of a 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 r.

Figure 2c shows that P(a) becomes more sharply peaked as λA increases; this makes sense as greater values of λA ensure the area constraint is more effective and the distributions become sharply peaked around the average area. Thus, the standard deviation of A decreases as λA 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 λA. Thus, treating the area part of Equation 1 as a constant is justified. Figure 2d shows that μ almost linearly increases with λA. 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 λA in a particular system. However, μ also varies with T (see Figure 2—figure supplement 4). Thus, in contrast to what has been proposed elsewhere (Wilk et al., 2014), as T changes, though the PDF of a 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 r 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 r. Figure 3a and b show representative plots for the comparison of the PDFs of r within the CPM, and the VM, respectively, where the lines represent fits with Equation 3. Figure 3a shows data with varying λP, and Figure 3b shows data with changing T. 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 λP and 1/T, hence with λP/T. 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, λA, λP, and T, keeping the other two fixed. First, when λA increases, the value of μ increases, but α remains almost constant (also see Figure 2b). Next, when λP increases, although α linearly increases, μ remains nearly the same. Finally, both parameters linearly increase with 1/T; since higher T implies more fluctuations, decreasing T helps both r and a 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 λA remains constant, varying λP and T 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 λP-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 λP and T. Additional junctional proteins may be employed during maturation to increase λP. On the other hand, different forms of activity may reduce, decreasing T. Since α increases linearly with λP/T, r alone is not enough to determine the dominant mechanism during the maturation process. However, assuming that λA 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: sd vs r¯. Figure 3f shows sd vs r¯ 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 αλP/T, higher α should correspond to slower dynamics. To test this hypothesis, we have simulated the CPM at different λP, P0, and T to obtain the relaxation time, τ (see Appendix 2 for details). From these control parameters, we have calculated α and then r¯, using our theory. We show τ as functions of α and r¯ in Figure 3g; it is clear that indeed τ grows as α increases or r¯ decreases. To show this behavior of τ on the same plot as sd vs r¯, we plot 2.5/lnτ-0.25 in Figure 3f as a function of r¯. 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 r¯; 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 P0 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 λP or 1/T 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.

Figure 4 with 1 supplement see all
Comparison with existing experiments.

(a) PDF for aspect ratio, r. Symbols are experimental data at three different times for MDCK cells, and lines represent the fits with our theory, Equation 3. The values of α are quoted in Table 1. Inset: PDF of the rs, the lines are theory, and the symbols are data. (b) Using the values of α, obtained for the three sets of data as quoted in Table 1, we obtain sd vs r¯ using our theory. Different symbols represent the types of the systems, and the colors blue, red, and black represent early, intermediate, and later time data, respectively. With maturation, the system moves towards lower r¯ and smaller sd. (c) PDF for rs for a wide variety of systems, shown in the figure, seems to be nearly universal, consistent with our theory. (d) The theory predicts a universal relation for sd vs r¯. Symbols are data for different systems, and the line is our theoretical prediction. (e) PDF of the scaled area for different epithelial systems and the lines represent the fits with Equation 4. (f) Predictive power of the theory: we use the r¯ for the three sets of cells from Kim et al., 2020, as marked by the hexagrams in (d), and obtain the corresponding values of α and obtain the PDFs for r using Equation 3. The colors correspond to the type of cells in (d), and the continuous line corresponds to the lower r¯ data. Inset: Lines are the theoretical PDFs for rs using the values of α for different cells, and the symbols show the experimental data. We have collected the experimental and some of the simulation data from different papers. Data taken from other papers are marked with a blue superscript in the legends. The sources are as follows: (1) Atia et al., 2018, (2) Li et al., 2021, (3) Lin et al., 2018, (4) Kim et al., 2020, (5) Fujii et al., 2019, (6) Ilina et al., 2020, (7) Wilk et al., 2014 and (8) Puliafito, 2017.

We next calculate sd as a function of r¯ 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 r¯, 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 r¯ 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 sd vs r¯ 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, r¯+1/r¯ deviates slightly from the behavior 1/α. 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 k 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 sd vs r¯. 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 r¯. 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 r distinguish the effects of maturation on the two key parameters: λP and T. 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 a 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 r¯, 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 r¯ 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 r, 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 H, (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 r. 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 sd vs r¯ 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 k being around 2.5. We show that a rough approximation of our expression leads to the k-Gamma function; the slight variation in k comes from the fact that the PDF is not strictly but nearly universal. On the other hand, k2.5 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 r alone. λP describes the cortical properties, and T 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 a (Weaire, 1986). It is a Gamma function, described by a single parameter μ. (Wilk et al., 2014) has proposed that the PDF of a 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 λP, whereas α varies linearly with both λP and 1/T. This distinction, assuming λA 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 r, 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 r¯ 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 r¯ and analyze biological functions via r¯ (Wyatt et al., 2015; Bosveld et al., 2016). Our theory implies that such knowledge contains a wealth of information. One can obtain α from r¯, 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 r¯ 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, q, 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 q. 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 r. The PDF of r is described by a single parameter, α. As a result, sd vs r¯ 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 k2.5 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 protocol

To set up the calculation, we represent the cell perimeter in a coarse-grained description via n 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, r, we first obtain the gyration tensor related to the moment of inertia tensor, a 2×2 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, s12 and s22, the squared-radii of gyrations around the respective principal axes, and r=s1/s2. 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, s, around the center of mass, and we have s2=s12+s22 (Davis and Denton, 2018). The distribution of s2 is

(5) P(s2)=1Zρ=12δ(j=1nxjρ)δ(11ns2xx)eHP/kBTx˙ds,

where Z is the partition function, x˙=ρ=12j=1ndxjρ, the volume element, kB, the Boltzmann constant, and T, 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 s2, 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: ijΛijij, where ij is the length between two consecutive vertices, i and j , and Λij gives the line tension. Since the degrees of freedom in the VM are the vertices, considering Λij constant, we obtain Equation 1. The constant Λij 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 P0. Thus, the adhesion part becomes -λPK~P0jlj2, where K~ is a constant. At the cellular level description, this microscopic difference can be accounted for by a renormalized P0. Since it is unclear how to measure P0 in experiments, we mainly restrict our discussions on λP and T. Unless otherwise stated, we assume P0 remains constant. However, the theory does capture the variation in P0, as we show within both the VM and the CPM simulations (Figure 3h).

Within our coarse-grained description, the contractile term is Pi2=(jlj)2. 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 (jlj)2njlj2=νjlj2, where νO(n) is a constant. Note that this inequality is exact for the CPM, where we can consider lj=1 as the line element and ν=n, 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 n, 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 k2.5 has a different origin that is not affected by this inequality. Then, the perimeter part of becomes νλ~Pjlj2=νλP(1-KP0)jlj2, where K=K~/ν. 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 r that we calculate via that of s. In the field of polymer physics, the distribution of s 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 x={x11,x12,x21,x22,xn1,xn2}, representing a set of n points on the perimeter. Then, we have P=νλ~Px(KI2)x, where x is the column vector, the transpose of x, K, the Kirchhoff’s matrix (Eichinger, 1977; Eichinger, 1980) with Kii=2 and K(i-1)i=Ki(i-1)=-1, and I2, the two-dimensional identity tensor. Thus, we have

(6) HP/kBT=γx(KI2)x,

where γ=νλP(1-KP0)/kBT. The distribution of s2 is

P(s2)=1Zρ=12δ(j=1nxjρ)δ(11ns2xx)eγx(KI2)xx˙ds,

where the squared radius of gyration is s2=n-1xx. Since the radius of gyration does not depend on the coordinate system, we are allowed to chose one that diagonalizes K. Say the diagonal matrix is Λ, and q 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 s, defined as

(7) s=1Ni=1N(xixCM)2,

where N is the total volume element in the system with coordinates xi, and xCM 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 s12 and s22 are the squared-radii of gyrations around the respective principal axes. Thus, the aspect ratio, r, is obtained as r=s1/s2, assuming s1s2. 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 s. So, we first calculate the distribution for s2 and then, using this result, we obtain the distribution of r.

Equation 5 can be written in the normal coordinate system as.

(8) P(s2)=1Zρ=12δ(qnρ)δ(11ns2qq)exp(γq(ΛI2)q)q˙ds,

where we have used qnρxjρ that corresponds to the zero-eigenvalue mode of the matrix. Integrating over qnρ, we get rid of this zero-eigenvalue that gives translation. Thus,

(9) P(s2)=1Zδ(11ns2q0q0)exp(γq0(Λ0I)q0)q˙0ds,

where we have defined q0 as the 2(n-1) dimensional vector excluding the coordinates corresponding to the zero-eigenvalue. The normalization factor, Z, can be calculated exactly through the integration as

(10) Zexp(γq0(Λ0I2)q0)q˙0=(πγ)(n1)|Λ0|1.

Note that the integration in the calculation of P(s2) is around the boundary of the cell; to separate out the radial part, we now write the volume element in polar coordinate u. Then q˙0=n(n-1)s2(n-1)-1dsu˙ and q0q0/ns2=uu=1. Thus, we obtain from Equation 9

(11) P(s2)=Adβeiβeγns2u[(Λ0iβnγs2In1)I2]uu˙,

with A=(γπ)(n1)|Λ0|12πn(n1)s2(n1)1In1 is the identity matrix of rank n-1. Carrying out the integration over u, we obtain

(12) P(s2)=Adβeiβ(πγns2)(n1)|Λ0iβγns2I|.

Using the value of A, we obtain

(13) P(s2)=12πsdβeiβj=1n1(1iβγns2λj)=(γns2)n12πs|Λ0|dβeiβj=1n1(γns2λjiβ),

where λj’s are the eigenvalues of K. The integral in Equation 13 can be performed via the contour integral and the resultant solution can be written as.

(14) P(s2)=(γns2)n12πs|Λ0|2πikRes(λk),

where λk are the distinct eigenvalues of K and Res(λk) gives the residue at the pole λk. As we show below, the residues will have a term exp[-nγs2λk] and in the limit s2, only the smallest λk will contribute.

Since the cell perimeter must be closed-looped, K 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 β=-iγns2λ, and of order 2. Thus, we obtain the residue as

(15) Res=ddβ[eiβj=1n3(γns2λjiβ)]β=iγns2λ.

Let’s first take the derivative, with respect to β, of the numerator and write part of the residue as

(16) term1=ieγns2λ(γns2)n3j=1n3(λjλ).

Next, differentiating the denominator, we obtain the other part of the residue as

(17) term2=eiβ[i(γns2λ1iβ)2j=2n3(γns2λjiβ)+i(γns2λ2iβ)2j=1,j2n3(γns2λjiβ)+]|β=iγns2λ=ieγns2λ(γns2)n2j=1n3(λjλ)[1λ1λ+1λ2λ+1λ3λ+].

A comparison of term1 and term2, given by Equation 16 and Equation 17, respectively, shows that there is an extra factor of s2 in the denominator of term2. Thus, in the limit s2 term2 can be ignored compared to term1. Therefore, we obtain the distribution function for s2 as.

(18) P(s2)=|Λ0|n2γ2j=1n3(λjλ)s3eγnλs2Cs3eα~s2

where C is the normalization constant that we will fix later, and α~=γnλ. Note that the lowest eigenvalue for the n-dimensional Kirchoff’s matrix is proportional to 1/n, thus nλO(1).

Now, s2=s12+s22 and the aspect ratio r=s1/s2. Moreover, we have s1s2=A, where A 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, A can be taken as a constant. Therefore, using the last two relations in the first, we obtain s2=A(r+1/r). Thus, we obtain from Equation 18, the distribution of r as.

(19) P(r)=1N(r+1r)3/2(11r2)eα(r+1r),

where N is the normalization constant: N=Γ(5/2)/α5/2W(5/2)1F1(5/2,7/2,2α), where W(x)=2x/x, Γ(x) is the Gamma function, and F11(a,b,c) is the Kummer’s confluent Hypergeometric function (Erdelyi et al., 1953), αλP(1-KP0)/T.

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 k2.5 when Equation 19 is approximated as a k-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 sd and r¯ is quite complex, though for most practical purposes it can be taken as a straight line: sd0.71r¯-0.75.

Radius of gyration

Request a detailed protocol

There 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 s (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 s1s2 or s2. We have used this property to obtain the distribution of r in the analytical derivation of Equation 19.

Distribution for area

Request a detailed protocol

Equation 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, A. 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, A. Since confluency is a strong geometric constraint on A, we must include this constraint in the derivation of P(A). Now, as detailed in ‘Radius of gyration’, s2A. Therefore, Equation 18 gives.

(20) P(A)A3/2exp[-βA]

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 f(A)Aν (Weaire, 1986; Gezer et al., 2021). Using this in Equation 20, we obtain

(21) P(A)Aμ1exp[βA],

where we have defined μ=ν+5/2. From Equation 21, we obtain the average area, A¯=μ/β. Thus, we obtain the normalized distribution for the scaled area, a=A/A¯, as

(22) P(a)=μμΓ(μ)aμ1exp[μa].

This is the well-known k-Gamma function, defined in Aste and Di Matteo, 2008, and usually denoted with the variable k. 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 k, which is obtained fitting the rs data, we have used μ to define the distribution function for a. Since μ comes from the constraint of confluency, it should be independent of λP. 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 E, given by Equation 1, should have the weightage of exp[-E/kBT]. However, as discussed below in Appendix 2, we restrict ourselves in the low-P0 regime, where the cells, on average, cannot satisfy the perimeter constraint. This implies a nontrivial dependence of the distribution around E0 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, P(E)exp(-E/T), plots of log[P(E)] as a function of E/T for different T should follow a master curve that is a linearly decreasing function. As shown in Appendix 1—figure 1, except at very close to E0, P(E) as a function of E/T at different T indeed follows a master curve, implying the existence of the Boltzmann distribution.

Appendix 1—figure 1
The systems in our simulation follow a Boltzmann distribution.

If the system follows a Boltzmann distribution at a temperature T, plot of the logarithm of distribution log[P(E)] as a function of E/T at different T should follow a master curve that linearly decreases with E/T. In the main figure, we show the distribution, P/E at different T and the inset shows the data collapse for P/E as a function of E/T. The data presented in this figure are obtained from the CPM simulation of a total system size 120X120 with 360 cells. We have set P0=26 and A0=40. We have equilibrated the system for 105 time steps before collecting data.

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-P0, and the large-P0 regime. However, the large-P0 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 P0 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, P(σσ)=e-ΔT, where we have set Boltzmann constant kB 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 M such elementary moves, where M 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 120×120 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 5×8. We first equilibrate the system for 8×105 MC time before collecting the data of aspect ratio. The variables that characterize the system are λA, λP, P0, and T. 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, Q(t), defined as

(23) Q(t)=1Nσ=1NW(a|Xcmσ(t+t0)Xcmσ(t0)|)t0¯,

where Xcmσ(t) is the center of mass at time t of a cell with index σ, W(x) is a Heaviside step function

(24) W(x)={1if x00if x<0,

t0 denotes averaging over initial times and the overline implies ensemble average. Relaxation time is defined as Q(t=τ)=0.3. The parameter a is related to the vibrational motion of the particles (here cells) inside the cage formed by its neighbors. We have set a=1.12 (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 128×128 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 T 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 δr. We accept the move with a probability P(CC)=eΔHT where we have set Boltzmann constant kB to unity, Δ is the change in energy going from CC. Initially, we start will 1024 cells having equal area (=1.0) and perimeter (=3.72). We have equilibrated the system for 2×106 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 P0 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 P0 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 P0 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 (r) in our simulation

As stated in the main text, we characterize the cell shape via aspect ratio, r. To calculate r, we follow (Atia et al., 2018) and briefly discuss the method here. Consider one particular cell, represented by the set of points {xi1,xi2}, as schematically shown in Appendix 3—figure 1. We first calculate the gyration tensor, I, (proportional to the moment of inertia tensor) in a frame of reference whose origin coincides with the center of mass (xc1,xc2) of the cell. I is a 2×2 tensor in spatial dimension two. We next diagonalize I and obtain the two eigenvalues, s12 and s22. Without loss of generality, we assume that s1s2, and obtain r=s1/s2. 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 3—figure 1
The process to calculate the aspect ratio of a Cell.

We first obtain the moment of inertia, I, in a frame of reference whose origin coincides with the center of mass of the cell, then diagonalize I, and calculate r, as the square root of the ratio of the two eigenvalues, s21 and s22, assuming s1 ≥ s2.

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 P(r). 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, kd-1, and apoptosis, ka-1, 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 (32) for selecting a cell for division. For apoptosis, every ka time step, we randomly choose a cell and assign it a target area A0=0. We also relax the constraint that does not allow fragmentation in our simulation for the cells undergoing apoptosis.

We show the PDF of r at three representative values of ka above. As shown in Appendix 4—figure 1a PDFs of r at different values of kd=ka 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 4—figure 1
PDF of r at three representative values of ka.

(a) PDFs of r at different values of kd = ka (quoted in b) within CPM for λP = 0.5 and T = 12.5. (b) PDF of rs corresponding to the same data as in (a). The lines are fits with Equation 3 of the main text. (c) α vs kd or ka quickly reaches a constant with increasing ka.

Appendix 5

Validity of the approximation of perimeter

We have used the Cauchy-Schwartz inequality, (jlj)2njlj2=νjlj2, where νO(n) is a constant. Within the CPM, the perimeter consists of n elementary sides, as schematically shown in Appendix 5—figure 1a above, then Pi2=[j=1nj]2. Since j=1, Pi2=[j=1nj]2=nj=1nj2. 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 n. This fluctuation depends on the values of the model parameters, such as T, λP, P0, 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 n, is less than 5%. Thus, we are justified to approximate the inequality by equality.

Appendix 5—figure 1
Test of the perimeter inequality.

(a) In the CPM, the approximation of P2i as njlj2, where n is the number of elementary lengths, ℓj is exact as can be seen from plot where the two estimates are equal. (b) For the continuum models, such as the vertex model, we can use a similar discretization as schematically shown in the figure. Experimental estimate of the perimeter comes from such a discretization, where the continuous cell perimeter is represented by a discrete line, as schematically shown with the red line. Even in this case, we find that the inequality is very close to an equality.

Data availability

We have uploaded the source files of the simulation data and the Mathematica analysis files in Dryad.

The following data sets were generated
    1. Nandi S
    (2022) Dryad Digital Repository
    Data from: On the origin of universal cell shape variability in a confluent epithelial monolayer.
    https://doi.org/10.5061/dryad.xsj3tx9h0

References

  1. Book
    1. Arfken GB
    2. Weber HJ
    3. Harris FE
    (2018)
    Mathematical Methods for Physicists (7th edition)
    Elsevier.
  2. Book
    1. Erdelyi A
    2. Magnus W
    3. Oberhettinger F
    4. Tricomi FG
    (1953)
    Higher Trascendental Functions
    New York: McGraw-Hill.
    1. Glazier JA
    2. Graner F
    (1993) Simulation of the differential adhesion driven rearrangement of biological cells
    Physical Review. E, Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 47:2128–2154.
    https://doi.org/10.1103/physreve.47.2128
  3. Book
    1. Swat MH
    (2012)
    Glazier
    In: Swat MH, editors. Computational Methods in Cell Biology. Academic Press. pp. 325–366.
  4. Software
    1. Wolfram Research Inc
    (2019)
    Mathematica, version version 12.0
    Wolfram research inc.

Decision letter

  1. Karsten Kruse
    Reviewing Editor; University of Geneva, Switzerland
  2. Aleksandra M Walczak
    Senior Editor; CNRS LPENS, France
  3. Jean-François Joanny
    Reviewer; Institut Curie, CNRS UMR168, France

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "On the origin of universal cell shape variability in a confluent epithelial monolayer" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Jean-François Joanny (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

1) The authors consider a regime where the area is equal to the target or preferred area of the cells and the area part of the free energy can be neglected. Quantitatively check the validity of this approximation.

2) Justify, why effectively equilibrium statistics can be used in the analysis.

3) In Appendix A you replace the square of the perimeter by the sum of the squares of the lengths of the edges of the cell. Check the accuracy of this approximation quantitatively.

4) Make the text more accessible to biologists.

Reviewer #1 (Recommendations for the authors):

I do not have any recommendations for further analysis or computation. However, the way the manuscript is written is rather poor. The authors often use jargon and, at places, the text is repetitive and confusing. Below I give some examples that might help the authors to make their work more accessible to the broad readership of eLife.

The authors very often switch between AR and r for denoting the aspect ratio. The figures sometimes, but not always contain both notions. I would strongly suggest introducing one symbol only and using it throughout the manuscript.

2nd paragraph is very long and partly repetitive.

p.2

"In a jammed granular system, the statistics of tessellated volume, x, follows a k-Γ distribution," This sentence is not very biologist friendly.

"Almost universal" Other expressions used are "nearly universal", "virtually universal". Neither of them is well-defined and it is not clear if they all refer to the same feature. If yes, use only one notion. Also, try to make it scientific and not just intuitive.

"The existing theoretical works on cellular shapes are based on the elasticity framework." Again, hard to understand for a biologist.

p.3

On this page there a several jumps in the presentation. Try to streamline.

"One crucial aspect of these systems" Which ones?

"the cortex, a thin layer of cytoplasm" Better: the cortex is a thin actomyosin layer

"we mostly focus our discussions below within the CPM " but then "cells are compact objects (in contrast to being fractal-like in the CPM simulation [30])" Can you clarify the compact vs fractal nature of the cells in the CPM you use?

"we first need to obtain the two radii of gyrations, s1 and s2, around the two principal axes." You should probably explain 'radius of gyration' for biologists as well as their relation to cell area.

"as detailed in the Appendix B, Equation (??), presented in the Appendix," Number is missing.

"Since μ is related to the constraint of confluency" Here Ii was confused, because you consider an isolated cell, right?

p.4

"Note that the power 3/2 of the algebraic term in Equation (3) comes from the mathematical property of closed-looped objects." Again, this requires more explanation for biologists

Figure 1d you could maybe show a log-log plot.

p.5

"Figure 3g" It seems to be premature to refer to this figure here.

"Finally, we show" Where?

"although they follow the same distribution, they need not be identical" What are 'they'? Also, I find it hard to understand this sentence.

"we have simulated the confluent systems with varying λ_A and find that the AR distribution remains almost independent of λ_A" Can you show the distribution?

"as T changes, the PDF of a, though well-described by a single parameter μ via Equation (4), can still be different" Yo mean the values of \mu can be different?

Figure 4c-e Why do you not show all curves for all cell types?

Figure 4d: Hexagons? Maybe stars?

p.10

last paragraph of the main text is just a repetition of the previous paragraphs

Reviewer #2 (Recommendations for the authors):

– The authors consider a regime where the area is equal to the target or preferred area of the cells and the area part of the free energy can be neglected. Can the validity of this approximation be made more quantitative? A naive guess would be that the modulus \λ_a is very large but in the end the authors calculate fluctuations of the area due to shape fluctuations and they would cost a large energy which is not taken into account. In the vertex model, as explained by the authors there is a transition between a soft solid and a hard solid behavior and the transition only depends on the ratio between the perimeter and the square root of the area. Which regime do the authors consider?

– The authors define the tissue by an energy which is used in all the models that they quote. In my mind this energy is an efficient way to write that the forces in the epithelial layer are balanced. However, they then use a Boltzmann Boltzmann statistics with a well defined temperature for the configuration of the tissue. This temperature seems to play an important role as it enters the expression for the parameter \α. Is this equilibrium statistics also used in the simulations? In my mind, a tissue is an intrinsically non-equilibrium system and I do not see how the equilibrium distribution can be justified. The authors also claim that the results are not changed if cell division is considered, which is difficult for me to understand. The manuscript would be much stronger if the authors were discussing this effective equilibrium assumption in details.

– Another approximation is made in appendix A which replaces the square of the perimeter by the sum of the squares of the lengths of the edges of the cell. I understand that this strongly simplifies the calculations. Could the authors make more quantitative this approximation which is only justified by an inequality. How accurate is it?

– On more general grounds, the results obtained by the authors are interesting but do they have any biological relevance? Could the authors point out a biological phenomenon for which cell shape fluctuations would play an important role?

https://doi.org/10.7554/eLife.76406.sa1

Author response

Essential revisions:

Reviewer #1 (Recommendations for the authors):

I do not have any recommendations for further analysis or computation. However, the way the manuscript is written is rather poor. The authors often use jargon and, at places, the text is repetitive and confusing. Below I give some examples that might help the authors to make their work more accessible to the broad readership of eLife.

We sincerely thank the reviewer for the insightful comments, appreciation of the importance of the work, and suggestions to make the manuscript accessible to a broader audience. We have revised the manuscript accordingly and removed the repetitions. We hope the reviewer and the editor finds the revised manuscript accessible to a broad audience. Here we briefly respond to the comments of the reviewer.

The authors very often switch between AR and r for denoting the aspect ratio. The figures sometimes, but not always contain both notions. I would strongly suggest introducing one symbol only and using it throughout the manuscript.

Thanks, we have now used the notation r to denote the aspect ratio.

2nd paragraph is very long and partly repetitive.

We have now shortened this paragraph and removed the repetitions.

p.2

"In a jammed granular system, the statistics of tessellated volume, x, follows a k-Γ distribution," This sentence is not very biologist friendly.

We have expanded this sentence, explaining what are these volumes: “In a jammed system, one can carry out …. by a single parameter k [41].”

"Almost universal" Other expressions used are "nearly universal", "virtually universal". Neither of them is well-defined and it is not clear if they all refer to the same feature. If yes, use only one notion. Also, try to make it scientific and not just intuitive.

We thank the reviewer for this suggestion. They do refer to the same feature and we now use the term “nearly universal” alone. We have also explained what it means (Page 5, left column): “the α-dependence becomes so weak … independent of α.”

"The existing theoretical works on cellular shapes are based on the elasticity framework." Again, hard to understand for a biologist.

We have now rephrased this sentence explaining the elasticity framework: “The existing theoretical works … when external force is relaxed.” (Page 2, first column).

p.3

On this page there a several jumps in the presentation. Try to streamline.

Thanks for pointing this out. We have now made it more contiguous by providing the details, either in the Appendix or within the main text itself. We have also streamlined the presentation by providing more details in the Materials and methods Section and including additional details in the Appendix.

"One crucial aspect of these systems" Which ones?

We have explicitly mentioned the models now: “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”.

"the cortex, a thin layer of cytoplasm" Better: the cortex is a thin actomyosin layer

We have rephrased it as “First, a thin actomyosin layer, known as cortex, mainly governs the cellular mechanical properties”.

"we mostly focus our discussions below within the CPM " but then "cells are compact objects (in contrast to being fractal-like in the CPM simulation [30])" Can you clarify the compact vs fractal nature of the cells in the CPM you use?

There are two regimes of these models, the system behaves like soft and hard solids in these two regimes. We mainly focussed on the CPM, because from the theoretical perspective, this model has a lot of advantage and easier to understand many of the complex behaviours. Within the CPM, the cells look like fractal objects and compact objects, respectively in these two regimes. We have added a subsection in Appendix C, titled “Two different regimes”, explaining this behaviour (page 15, last subsection of Appendix 2).

"we first need to obtain the two radii of gyrations, s1 and s2, around the two principal axes." You should probably explain 'radius of gyration' for biologists as well as their relation to cell area.

We have added a section, “Radius of gyration”, in the Materials and methods Section with the various definitions of radii of gyration and their relation with the area (Page 13, Section IV B).

"as detailed in the Appendix B, Equation (??), presented in the Appendix," Number is missing.

Thanks for pointing this out, we have fixed it in the revised manuscript.

"Since μ is related to the constraint of confluency" Here Ii was confused, because you consider an isolated cell, right?

We apologise for this confusion. Actually, there are two different calculations presented in this work: the distribution of aspect ratio and the distribution of the cell area. The first one is the main result of the work, but the second result is also important as it connects the derivation to several other problems.

Our interest is the cellular aspect ratio in a confluent monolayer. But the constraint of confluency is an extremely hard mathematical problem and no exact results exist till date to treat it analytically. Motivated by our simulation results, we have argued that the constraint of confluency is not crucial for the aspect ratio distribution. We have also verified this in our simulations. This slightly simplifies the derivation of P(r).

But, for the distribution of cellular area, the constraint of confluency is crucial. Equation

(18) is valid for a general system. We then included this constraint in our analytical derivation following Weaire et al. (Ref. 67) and obtained the distribution of cell area.

We have clarified this in the text in section Materials and methods, IV C.

p.4

"Note that the power 3/2 of the algebraic term in Equation (3) comes from the mathematical property of closed-looped objects." Again, this requires more explanation for biologists.

We have detailed this point in the Appendix A. In the main text (page 4, second paragraph), we have now stated the property that leads to this value and has referred the reader to the appendix: “That is, for closed-loop objects, …. As detailed in Sec. IV A”.

Figure 1d you could maybe show a log-log plot.

We have now shown the plot in log-log scale in the inset of Figure 1(d). p.5

p.5

"Figure 3g" It seems to be premature to refer to this figure here.

We have removed the reference to Figure 3(g) here.

"Finally, we show" Where?

Thanks, we have now added the reference [Figure 1(f)].

"although they follow the same distribution, they need not be identical" What are 'they'? Also, I find it hard to understand this sentence.

We have now clarified this sentence: “The PDFs of a has been argued… they are not identical”.

"we have simulated the confluent systems with varying λ_A and find that the AR distribution remains almost independent of λ_A" Can you show the distribution?

We have now included these distributions in the supplementary figure, Figure 2 —figure supplement 2 and added in the main text: “(the distributions are shown in Figure 2 —figure supplement 2).”

"as T changes, the PDF of a, though well-described by a single parameter μ via Equation (4), can still be different" Yo mean the values of \mu can be different?

Yes. We have now rephrased this sentence for better clarity: “… by the same function, Equation (4), …. need not be identical”.

Figure 4c-e Why do you not show all curves for all cell types?

We are a bit confused with this question. We do have the data corresponding to all the cell types mentioned in the figure.

On the other hand, if the reviewer means plotting the lines for all the cell types: we represented most data with symbols as plotting the lines for all the cell types makes the figure extremely hard to understand.

Figure 4d: Hexagons? Maybe stars?

Thanks for pointing this out. Actually these symbols are called Hexagrams that we incorrectly wrote Hexagons, we have corrected this in the revised manuscript.

p.10

last paragraph of the main text is just a repetition of the previous paragraphs

We have now removed the repetitions. The last paragraph now summarises the achievements of the work.

Reviewer #2 (Recommendations for the authors):

– The authors consider a regime where the area is equal to the target or preferred area of the cells and the area part of the free energy can be neglected. Can the validity of this approximation be made more quantitative? A naive guess would be that the modulus \λ_a is very large but in the end the authors calculate fluctuations of the area due to shape fluctuations and they would cost a large energy which is not taken into account. In the vertex model, as explained by the authors there is a transition between a soft solid and a hard solid behavior and the transition only depends on the ratio between the perimeter and the square root of the area. Which regime do the authors consider?

We thank the reviewer, Prof. Joanny, for the insightful comments, questions, and suggestions. We have included all the suggestions in the revised manuscript. Here, we briefly respond to the comments and questions.

There are two distinct analytical results in the paper: the distributions of aspect ratio and cell area. We have argued that the area constraint (or the constraint of confluency) is not crucial for the aspect ratio distribution, but it is extremely essential for the distribution of the area. As detailed below, we have made this argument more quantitative with additional results in the revise manuscript.

As we discussed in the main text, this approximation is motivated by our simulation results:

– First, as shown in Figure 2(c), the area distribution is sharply peaked around the average value. As the reviewer correctly points out, the sharpness of this distribution grows as λA increases.

– Second, we find that even at moderate values of λA (=0.1 — 0.8), where the area distribution is not a δ-function, the value of α, which determines the aspect ratio distribution, does not change (Figures 2b). We now show these distributions as well in the new supplementary Figure 2 —figure supplement 2.

– These results are consistent with the expectation that aspect ratio, being a dimensionless parameter, should be independent of the actual value of cellular area. Our simulation results, presented in Figure 2(b), and the additional results presented in the supplementary figures, Figure 2 —figure supplement 1 and Figure 2 —figure supplement 2, show that the area fluctuation does not significantly affect the shape fluctuation. Therefore, we are justified to neglect the area term in the energy, Equation (1), to obtain the aspect ratio distribution.

– On the other hand, as you also pointed out, for the other result, that is the area fluctuation due to shape fluctuation, there is indeed a large cost. We have taken this into account in an indirect manner (the discussion between Eqs. (20) and (21) in Sec. IV C), considering the very strong geometric constraint, that is the constraint of confluency. As you know, this is a hard mathematical problem and no exact results exist. We have used the results of Weaire et al. [Ref. 67] in the form of a constraint. This is a phenomenological approach, our analysis for the area distribution can be viewed as a support of this approach.

– We fixed the target perimeter, P0, in our simulations, except for the data in Figure 2(h), such that the system is in low-P0 regime. We believe that the low-P0 (hardsolid-like) regime is more relevant experimentally. But, as shown in Figure 2(h), the theory seems to capture the variation with P0 as well. We have now highlighted these points in the revised manuscript, and have specifically added a subsection (Two different regimes) in Appendix 2 discussing this.

– The authors define the tissue by an energy which is used in all the models that they quote. In my mind this energy is an efficient way to write that the forces in the epithelial layer are balanced. However, they then use a Boltzmann Boltzmann statistics with a well defined temperature for the configuration of the tissue. This temperature seems to play an important role as it enters the expression for the parameter \α. Is this equilibrium statistics also used in the simulations? In my mind, a tissue is an intrinsically non-equilibrium system and I do not see how the equilibrium distribution can be justified. The authors also claim that the results are not changed if cell division is considered, which is difficult for me to understand. The manuscript would be much stronger if the authors were discussing this effective equilibrium assumption in details.

The confluent models used in this work are well-established models for epithelial monolayers. These models have two distinct variants: depending on the presence or absence of activity. When there is no activity, the models are in equilibrium. The nonequilibrium models can have several interesting complex behaviours, as shown by you in the recent PRL paper (Ref. 73). But, the equilibrium models are also significant and capture the key static and dynamic properties of such systems, as shown by several groups including Jülicher, Fredberg, etc (Refs. 31, 52, 56, 65, etc.). How far can these equilibrium models represent a cellular system is an important question that we did not address in this work. This is also an important question in general. The role of activity is also significant, but for a cleaner understanding, we have not considered the effects of activity in the current paper. One crucial question we are currently exploring is whether an effective equilibrium behaviour is possible even under activity, similar to what we found in some other contexts of glassy dynamics. In fact, a single particle effective temperature description seems to capture many aspects in a dense glassy system [Nandi and Gov, Soft Matter (2017)]. We intend to explore this question in the future.

In the current work, in the absence of cell division and apoptosis, we are justified to use the Boltzmann distribution, as the temperature in these models has the same status as an equilibrium temperature. We did not consider any Boltzmann distribution in our simulations, but we now show in the additional data ( Appendix 1) that this is a reasonable approximation.

For the second point, we are sorry for the confusion. It is not yet clear to us how to develop an analytical description in the presence of cell-division and apoptosis, one possibility is to include these rates following this work: https://www.pnas.org/doi/ 10.1073/pnas.1011086107. We are also trying to apply the non equilibrium approach that you developed in the recent PRL (Ref. 73). Our point in the current manuscript is simply the following: when the rates of these processes are small, the basic characteristic of the system should remain similar. The analytical derivation suggests that the aspect ratio distribution will have two different parts: the algebraic part is determined by the geometric property that cell perimeter is closed-looped.

This cannot change. So, all system specific details should enter via α alone. What we showed in the manuscript is that this expectation seems to be reasonable when we include these processes in our simulations, with their rates being not very large. We have now clarified this in the manuscript (Page 9, Column 2, second paragraph): “We have neglected …. processes are not dominant (Appendix 4).”

– Another approximation is made in appendix A which replaces the square of the perimeter by the sum of the squares of the lengths of the edges of the cell. I understand that this strongly simplifies the calculations. Could the authors make more quantitative this approximation which is only justified by an inequality. How accurate is it?

Thanks for raising this question. This is indeed an important trick in our analytical derivation. We have now discussed this approximation in detail in Appendix 5 justifying this approximation. We have also discussed this point in Sec. IV A, Page 11: “Note that this inequality … integral slightly easier.”

– On more general grounds, the results obtained by the authors are interesting but do they have any biological relevance? Could the authors point out a biological phenomenon for which cell shape fluctuations would play an important role?

We believe these results can pave the way to understand pattern formation in a new light. Apart from providing explanations for several existing experimental results, the analytical result is important as we can now formulate a theoretical framework to understand the mechano-chemical pattern formation in confluent epithelial tissues in terms of this analytical result. Many past works have shown that cellular aspect ratio is directly correlated with biological functions, such as differentiation and division. But all these studies are either at the single cell level or in terms of the averages in a monolayer. Having the analytical distribution of the aspect ratio, we can now develop a statistical description of the monolayer. For example, we are now studying how to couple the model with different gradients of morphogens (such as the wingless and Dpp in the Drosophila wing disc). Having such a description will then allow us to study the pattern formation in such systems.

We have included a brief discussion of this aspect in the conclusion (Pages 9-10): “Pattern formation in a biological system … in the Drosophila wing disc.”

https://doi.org/10.7554/eLife.76406.sa2

Article and author information

Author details

  1. Souvik Sadhukhan

    Tata Institute of Fundamental Research, Hyderabad, India
    Contribution
    Software, Formal analysis, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0926-6601
  2. Saroj Kumar Nandi

    Tata Institute of Fundamental Research, Hyderabad, India
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Investigation, Visualization, Writing – original draft, Writing – review and editing
    For correspondence
    saroj@tifrh.res.in
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0977-1947

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.

Senior Editor

  1. Aleksandra M Walczak, CNRS LPENS, France

Reviewing Editor

  1. Karsten Kruse, University of Geneva, Switzerland

Reviewer

  1. Jean-François Joanny, Institut Curie, CNRS UMR168, France

Publication history

  1. Preprint posted: August 21, 2021 (view preprint)
  2. Received: December 15, 2021
  3. Accepted: December 22, 2022
  4. Accepted Manuscript published: December 23, 2022 (version 1)
  5. Version of Record published: January 11, 2023 (version 2)

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

  • 1,177
    Page views
  • 165
    Downloads
  • 0
    Citations

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

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

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)

  1. Souvik Sadhukhan
  2. Saroj Kumar Nandi
(2022)
On the origin of universal cell shape variability in confluent epithelial monolayers
eLife 11:e76406.
https://doi.org/10.7554/eLife.76406

Further reading

    1. Developmental Biology
    Greta Rossi, Gabriele Ordazzo ... Vania Broccoli
    Research Article Updated

    Wolfram syndrome 1 (WS1) is a rare genetic disorder caused by mutations in the WFS1 gene leading to a wide spectrum of clinical dysfunctions, among which blindness, diabetes, and neurological deficits are the most prominent. WFS1 encodes for the endoplasmic reticulum (ER) resident transmembrane protein wolframin with multiple functions in ER processes. However, the WFS1-dependent etiopathology in retinal cells is unknown. Herein, we showed that Wfs1 mutant mice developed early retinal electrophysiological impairments followed by marked visual loss. Interestingly, axons and myelin disruption in the optic nerve preceded the degeneration of the retinal ganglion cell bodies in the retina. Transcriptomics at pre-degenerative stage revealed the STAT3-dependent activation of proinflammatory glial markers with reduction of the homeostatic and pro-survival factors glutamine synthetase and BDNF. Furthermore, label-free comparative proteomics identified a significant reduction of the monocarboxylate transport isoform 1 (MCT1) and its partner basigin that are highly enriched on retinal glia and myelin-forming oligodendrocytes in optic nerve together with wolframin. Loss of MCT1 caused a failure in lactate transfer from glial to neuronal cell bodies and axons leading to a chronic hypometabolic state. Thus, this bioenergetic impairment is occurring concurrently both within the axonal regions and cell bodies of the retinal ganglion cells, selectively endangering their survival while impacting less on other retinal cells. This metabolic dysfunction occurs months before the frank RGC degeneration suggesting an extended time-window for intervening with new therapeutic strategies focused on boosting retinal and optic nerve bioenergetics in WS1.

    1. Developmental Biology
    Kavitha Chinnaiya, Sarah Burbridge ... Marysia Placzek
    Research Article

    The tuberal hypothalamus controls life-supporting homeostatic processes, but despite its fundamental role, the cells and signalling pathways that specify this unique region of the CNS in embryogenesis are poorly characterised. Here we combine experimental and bioinformatic approaches in the embryonic chick to show that the tuberal hypothalamus is progressively generated from hypothalamic floor plate-like cells. Fate-mapping studies show that a stream of tuberal progenitors develops in the anterior-ventral neural tube as a wave of neuroepithelial-derived BMP signalling sweeps from anterior to posterior through the hypothalamic floor plate. As later-specified posterior tuberal progenitors are generated, early-specified anterior tuberal progenitors become progressively more distant from these BMP signals and differentiate into tuberal neurogenic cells. Gain- and loss-of-function experiments in vivo and ex vivo show that BMP signalling initiates tuberal progenitor specification, but must be eliminated for these to progress to anterior neurogenic progenitors. ScRNA-Seq profiling shows that tuberal progenitors that are specified after the major period of anterior tuberal specification begin to upregulate genes that characterise radial glial cells. This study provides an integrated account of the development of the tuberal hypothalamus.