E2 distribution and statistical regularity in polygonal planar tessellations

From solar supergranulation to salt flat in Bolivia, from veins on leaves to cells on Drosophila wing discs, polygon-based networks exhibit great complexities, yet similarities persist and statistical distributions can be remarkably consistent. Based on analysis of 99 polygonal tessellations of a wide variety of physical origins, this work demonstrates the ubiquity of an exponential distribution in the squared norm of the deformation tensor, $E^{2}$, which directly leads to the ubiquitous presence of Gamma distributions in polygon aspect ratio. The $E^{2}$ distribution in turn arises as a $\chi^{2}$-distribution, and an analytical framework is developed to compute its statistics. $E^{2}$ is closely related to many energy forms, and its Boltzmann-like feature allows the definition of a pseudo-temperature. Together with normality in other key variables such as vertex displacement, this work reveals regularities universally present in all systems alike

Polygonal networks are one of nature's favorite ways of organizing the multitude -from supergranulation on the solar surface [29] to cracked dry earth [37]; from ice wedges in northern Canada [61] to the scenic Salar de Uyuni in Bolivia [23]; and from veins on leaves [11] to cells on Drosophila wing discs [53] (Fig. 1).These systems are driven by distinctive physical mechanisms, yet they share common features.Individual constituents, namely, "cells" appear to "randomize" into statistical distributions, and only interact with their immediate neighbors.On the collective level, especially in the dynamic and active systems, rich phenomena are observed, including unjamming and jamming, fluid-to-solid phase transition, and flow and migration [7,19,22,34,42,50,54,63,68].
Despite the complexity and variabilities involved in these phenomena, similarity patterns emerge.One particularly interesting instance is provided recently by Atia et al. [3].Within the context of confluent biological tissue and based on extensive experiments both in vitro and in vivo, the authors found that data  [37]; (b) Ice wedges from northern Canada [61]; (c) Veins on a Ficus lyrata leaf [11]; (d) Desiccation pattern of an ancient lake on Mars [43]; (e) A snapshot of Salar de Uyuni in Bolivia, world's largest salt flat [23]; (f) Supergranulation on solar surface [29]; (g) Plated MDCK cells (this work).(h) A processed image of a developing Drosophila wing disc (this work) with the aspect ratio of cells color-mapped.(i) A regular hexagon P R (blue, with vertices x j ), a deformed hexagon P (red, with vertices y j ), and a uniform deformation as mean-field approximation P (black dashed outline, with vertices y j ). on cell aspect ratio collapse and follow a normalized Gamma distribution, implying a universal principle underlying the geometric configuration and pertinent processes.
What is the basis of this universality?Does it carry beyond the biological context?This work composes of a two-part discovery to answer these questions.In the first, we analyze a total of 99 data sets in 8 groups that include convection patterns (solar supergranulation), landforms (salt flats, on Mars, and in or near the Arctic), cracked dry earth, and biological patterns (veins on leaves, cells on Drosophila wing discs, and plated MDCK cells).(Fig. 1 and Table 1.)We demonstrate that the Gamma distribution in polygon aspect ratio derives from an exponential distribution in the squared strain tensor norm, E 2 , based on which we develop a unifying solution for the former.Both distributions persist in all data examined.In the second, we tackle the origin of the E 2 distribution, and demonstrate that it arises as a χ 2 -distribution owing to asymptotic normality of vertex displacement and other key variables.We present a theoretical framework to accurately compute E 2 from vertex statistics.Importantly, E 2 is closely related to various definitions of system energy including all of the bulk-, perimeter-, and moment-based (known as the Quantizer) forms [7,15,27,34,36,44].The Boltzmann-like feature of E 2 enables the definition of a pseudo-temperature that can be considered as a consistent quantifier [4,17,18,47].The exponential and normal distributions we reveal in this work are hidden regularities that transcend the specific physical systems, and are embedded universal features of random polygonal tessellations.

E 2 and its relationship with aspect ratio
In this first part of this work, we define E 2 and analytically establish its relationship with the polygon aspect ratio.We demonstrate that a k-Gamma distribution in the latter is derivable from an exponential distribution in the former.We begin by defining the mean-field deformation tensor, E.An exemplary processed image of a Drosophila wing disc 120h after egg laying (AEL) is shown in Figure 1h, where the color scale indicates magnitude of the cell aspect ratio (defined in Methods).We choose as our reference frame a regular n-polygon centered at the origin, with vertices x j = r 0 e j , e j = [cos(j2π/n), sin(j2π/n)], j = 1, ..., n. ( This polygon is denoted by P R , and r 0 is a scaling factor.Figure 1k uses n = 6, a hexagon as an illustrative example.We consequently regard any n-sided polygon P with vertices y j as a "deformation" from P R , Note that this deformation is to be understood as a morphological deviation from the reference P R , rather than an actual physical deformation, although the latter is a possibility.We also require that the centroid of P is aligned with that of P R .The deformation is in general non-uniform, in the sense that for n ≥ 3, a single deformation tensor of F ∈ R 2×2 cannot be identified by y j = Fx j for all of j = 1, ..., n.Nevertheless we can introduce an approximation, namely, This approximation is analogous to a Taylor expansion in which only the leading order term is retained. The use of a uniform deformation to approximate the local and non-uniform deformation field is effectively coarse-graining, reducing the degree of freedom from 2n to 4 and suppressing the fluctuations.This idea is illustrated in Figure 1i (right), where P is the approximate and uniformly deformed polygon.From F we define the usual strain tensor and its squared Frobenius norm, where I is the identity tensor, and the approximation in the first equation is valid in the small-to-moderate deformation regime.In this work we use |E| 2 (mathematical representation) and E 2 (terminology) interchangeably.The restriction to isochoric (area-conserved) deformation requires that det F = 1 and equivalently (to leading order), TrE = 0, which is satisfied by choosing r 0 in (1).Consequently, E has a small degree of freedom of 2, which eventually leads to its strong regularity that we demonstrate later.
We pursue an analytical expression for E via minimizing the difference between y j 's and y j 's in (3), from which we obtain (SI) Here In addition to satisfying the trace-free condition, the scaling also ensures that v j 's are dimensionless.
Further details on the computation of |E| 2 is deferred until later.For now we show that if E 2 follows an exponential distribution (that shall be validated below both via data and analysis), namely, then the aspect ratio follows a k-Gamma distribution such as shown in [3].Here ρ(•) denotes a (normalized) probability density function (PDF), and β is similar to an inverse temperature as the PDF is Boltzmannlike.The aspect ratio a r of the polygon is calculated via the second area moments (Methods), and is related to E 2 by (SI) where for convenience we define a shape factor x in relation to a r using x = a r − 1.Based on (8, 9), a where Here we have considered the normalization condition.Equation ( 10) is our prediction of the distribution in the aspect ratio.
Data validate E 2 and a r distributions Both distributions and their relationship are extensively validated with a total of 99 data sets spanning 8 groups summarized in Table 1; detailed descriptions, data sources, and method of analysis are presented in the Methods section.The center column shows the PDFs of the shape factor, a r − 1 (symbols).The theoretical predictions per Eq. ( 10) are shown in dashed, and exhibit excellent agreement with data.They are generated per Eq.(10), with the single input parameter, β, extracted from the analysis of |E| 2 distribution.
Lastly in the right column, the PDFs for a r − 1 are normalized with a r − 1, where • denotes a mean, e.g., Both data and theoretical predictions are normalized following this practice.The dot-dashed are best fits Type (abbreviation)  2); R 2 for a r − 1 indicates quality of agreement between theory and data (e.g., in middle column, Fig. 2).Data sources are presented in Methods.
using a k-Gamma distribution defined as where Γ(k) is the Gamma function, and k is the single fitting parameter.The agreement is evident, and the k-values are found to vary between 2 and 3.
Overall corroboration between theory and data is quantified by the coefficient of determination, R 2 , and are listed in Table 1 for all cases (see Methods for definition).The values are uniformly close to 1 ("perfect agreement") with minimal variations from case to case.Validity of the theoretical prediction (10) is also attested by Fig. 3a.Here we define a pseudo temperature T as the inverse of β, namely, T = β −1 .
To compare with data, a theoretical prediction is generated by using (10) in the integration (11).
The above results validate that the E 2 does follow an exponential, Boltzmann-like distribution.The universality of this distribution in all data sets, according to our theory, necessarily leads to a universality of k-Gamma distributions for the aspect ratio.That is, the validity extends beyond the confluent tissues studied in [3], and to all systems we analyzed.As a corollary, (10) provides a fundamental solution for the aspect ratio, of which the k-Gamma distributions ( 12) are convenient approximations.This solution does not normalize to a single curve with a single k value if fitted with Gamma distributions.In fact, it predicts a positive correlation between β and k (Fig. S8).Qualitatively, this means that PDFs of aspect ratio (a r , or equivalently, x) with wider spread and greater mean (corresponding to lower values of β) present themselves relatively to the left after normalization (corresponding to lower values of k and noting that maxima occur at 1 − 1/k per (12)).This trend is fully corroborated with data from our own work (Fig. 2, and Fig. S2 in SI, center and right columns) and Atia et al. [3] (Fig. 3 therein), as well as predictions from a self-propelled Voronoi model in the supplemental information of the latter.In summary, the variability in k arises from the variability in β.
We remark that the agreement between our theoretical prediction and data can be even better if we use a variation of the deformation tensor computed as the square-root of the area moment tensor, M, via Eq.( 21) in the Methods section.This is not surprising, as now both E and a r share the same origin, the agreement shown in Fig. 3a inset is near perfect.The slight differences between the two definitions are  ) for all 99 data sets (symbols); the theoretical prediction (dashed) is generated per (10).The inset alternatively presents E from a moment-based calculation which further improves agreement.(b) A comparison between data and predicted temperature, T n and T ("theory").The inset shows that sub-ensemble temperature T n 's are quantitatively similar to tessellation temperature, T .
due to higher order effects that we theoretically and numerically demonstrate in the SI.Here and below we focus on using (5) for its apparent analytical simplicity.
In the second part of this work, we demonstrate the origin of the highly regular statistical distribution in E 2 .Figure 3b shows results comparing the pseudo-temperature computed from Eq. ( 6) (denoted "data") with the theoretical prediction we develop below (denoted "theory").Here the subscript n denotes a restriction to the sub-ensemble of n-gons, T n := |E| 2 n .Polygons other than n = 4-8 are of statistically insignificant occurrences and not included in the evaluation.
The key relationship we utilize is a quadratic form to compute |E| 2 given vertex displacement, v j , Here v is a concatenated vector in R 2n , Other vectors such as ŷ, û, and ê are similarly defined from their two-dimensional counterparts, and all vectors and tensors in the 2n-dimensional space are denoted by a hat to differentiate from the planar quantities.Ĉ ∈ R 2n×2n sym is a second-order tensor with block components C ij given in (6).Not surprisingly, Ĉ has 2 non-trivial eigenvalues (SI), matching the degree of freedom of E (note that E = (E 11, E 12 ; E 12 , −E 11 ) in a general component form): PT Ĉ P = diag(2/n, 2/n, 0, ..., 0).
The diagonalization above with the orthonormal tensor P helps us express |E| 2 in a particularly simple form, We realize that ŵk = pk • v, pk being an eigenvector of Ĉ.If v is characterized by a covariance matrix, Σ, then the variance of ŵk is [33] Var and T n is readily calculated as Note we have used Ĉ = 2 n (p 1 ⊗ p1 + p2 ⊗ p2 ).Here and after and as a good approximation, we assume v has a zero mean.Equation ( 16) is a precise expression to compute T n given Σ, and is used to generate the theoretical predictions in Fig. 3b, main panel.The tessellation average T can be computed by taking the weighted sum of T n , namely, T = n (N n T n )/ n N n , where N n is the number of n-gons.On the other hand, sub-ensemble temperatures are typically quantitatively similar to the tessellation temperature, as shown in Fig. 3b inset.
If we further assume that ŵ1,2 follow identical normal distributions, immediately we have In other words, the exponential distribution arises actually as a χ 2 -distribution with 2 degrees of freedom.
On the other hand, if the variances Var( ŵ1,2 ) are not identical but quantitatively similar, which is true for all tessellations we study (see Fig. S5), Eq. ( 17) still holds to the leading order.(This point is straightforward to prove via Taylor expansion and not shown here for brevity.)Note that even in this situation, per ( 16) the formula for T n is still accurate without approximation.This provides an essential illustration of the origin of the E 2 distribution, and Eq. ( 17) is a main result of the current work.It remains to be shown below that ŵ1,2 distributions are indeed approximately normal and independent.The apparent candidate to explain the resulting normality is the central limit theorem in the generalized version for dependent and identical random variables [39], noting that ŵk derives from vk via a linear combination (15).It is peculiar to note that ûk and vk themselves also demonstrate approximate normality, shown in Fig. 4c and d.The normality in ûk is again be explained by the central limit theorem.We can write û, the concatenated vector for u j 's as (SI)

Asymptotic normality contributes to statistical regularity
In the absence of apparent anisotropy, components of ŷ − ŷ are approximately dependent yet identical, satisfying precondition of the theorem.Hence ûk asymptotes to normality.On the other hand, from Eq. ( 7) we have v = û/r 0 and r 0 = (ŷ • ê)/n.The normality in v is difficult to theoretically prove.However, it is reasonable to speculate the loss of the apparent scale would create similarity to preserve or even enhance normality -see also Figs 4f, g below.In addition, it is extensively confirmed by the data as shown in Fig. 4d.
The asymptotic normality can be better illustrated via a simple Monte Carlo simulation following the schematic in Fig. 1k, where we temporarily restrict to an isolated hexagon, and displacements ûk,0 (k = 1, 2, ..., 2n) are prescribed according to independent, identical distributions as shown in Fig. 4e (Methods, Eqs.(22,23)).Two representative cases are examined, the first with a steeper than Gaussian initial descent, and the second non-monotonic.In Figs.4f and g, both ûk and vk already demonstrate In a summary, the above exercises demonstrate that asymptotic normality is prevalent in planar tessellations, as key variables derive from linear combinations of statistically similar components.As a result E 2 distributions become highly regular due to combined normality and its low-dimensionality.

Physical Implications
The physical meaning of E is self-evident: it represents deformation, and hence is typically associated with energy in one form or another.In the wide range of phenomena we studied, the constitutive relations come in different forms (some are yet unknown).However, some usual possibilities can be contemplated.
If the energy is bulk-elastic in nature, then any physically reasonable elastic model of a polygon, valid in the small-to-moderate deformation regime, must follow the form [25] where µ is the first Lamé constant, whereas the second constant is not needed as TrE = 0 (SI).On the other hand, if energy is associated with edge lengths or perimeters, such as in the case of models for 2D confluent tissues [7,34,44], Eq. ( 18) is still a formally valid approximation, as the change in perimeter is also proportional to |E| 2 to leading-order approximation (SI).Last but not least, in the Quantizer problem [15,27,36] the cell-wise energy functional is the moment of inertia, which is in both two and three dimensions (SI), and 2m 0 is the moment of inertia of the regular reference polygon, P R .Thus its distribution can be computed via knowing both the area and E 2 distributions.These examples of constitutive relations cover a reasonably wide range of physical systems.
Above we have taken the liberty in naming a pseudo-temperature, T (or T n for the sub-ensembles).
Indeed, such definition is both tempting and appropriate in the presence of a Boltzmann-like distribution.
The tests by Dean and Lefèvre [13] and McNamara et al. [47] become trivial: the ratio of two overlapping exponential distributions will necessarily give another exponential distribution.We therefore name this pseudo-temperature the "E 2 temperature", and propose it as a candidate for a thermodynamically meaningful temperature owing to its consistent regularity and connection with physical quantities.This temperature quantifies the overall deformation, and possibly also system energy.To further explore a thermodynamic framework would require system-specific physical principles, e.g., energy minimization, which we shall explore in future work.
We have thus demonstrated three main points in this work: i) An exponential distribution in E 2 leads to a k−Gamma distribution in the aspect ratio.In fact, k-Gamma distributions are mere approximations to the more basic solution we develop.ii) E 2 distribution is a χ 2 -distribution with two degrees of freedom, arising from combined effects of asymptotic normality and the small dimensionality of E, which is analogous to the small dimensionality of the volume function in granular assembly [8].We have developed a formula to compute E 2 from vertex statistics.iii) E 2 and aspect ratio distributions as well as normality in displacements are true universal features as we have shown via both a large collection of data and theoretical derivations illustrating their mechanistic origins.The strong regularity in E 2 and vertex displacements are "hidden patterns" revealed by this work.The mean-field strain tensor, with its clear physical and geometric meaning, is an ideal quantity connecting the conservation principles, the energy (or pseudo-energy) landscape, and the geometric distributions.It is a powerful quantifier to describe polygonal networks randomized by active agitations, structural defects, and noises, among others.Analysis may also be extended to polytopes in three and higher dimensions.
Drosophila wing disc, fixed (Droso) Drosophila were cultured at 25Â°C.To obtain fixed wing discs at different stages, eggs were laid for 2 to 4h, and larvae were dissected at 72, 84, 96, 108 and 120h after egg laying (AEL).Dissected wing discs were fixed in 4% paraformaldehyde for 15 min at room temperature.
Staining of fixed wing discs was performed essentially as described in [58] using rat anti-E-cad (1:400 DCAD2; DSHB) and anti-rat Alexa Fluor 647 (Jackson ImmunoResearch, 712-605-153).Images were captured on a Leica SP8 confocal microscope.To compensate for aberrations due to the curvature of wing disc and signals from the peripodial epithelium, we used the Matlab toolbox ImSAnE [28] to detect and isolate a slice of the wing disc epithelium surrounding the adherens junctions, which was then projected into a flat plane, as described previously [53].
Drosophila wing disc, live (Droso) For live imaging of cultured wing discs, larvae expressing GFPlabelled E-cadherin from a Ubi-Ecad:GFP transgene were dissected at 96h  [30].Images were acquired using LAS X software on a Leica TCS SP8 confocal microscope system using a HC PL APO 63×/1.40objective.

Image and data analysis
Fluorescent images are analyzed using Tissue Analyzer, a plug-in of ImageJ (version 1.52j), from which the cells are sectioned and cell centroids, edges, and vertices are identified.Post-processing is then performed with MATLAB.For each cell, the second area moment tensor, defined with respect to the cell centroid c, where the integration is over the polygon (cell) P. Note that here we ignore the curvature of cell edges and assume (by approximation) that they are straight lines connecting vertices.Standard and exact formulae are available for polygons which we use to compute the components of M with only the coordinates of the vertices, y j 's.The aspect ratio a r is a r = max(λ 1 , λ 2 ) min(λ 1 , λ 2 ) , where λ 1 and λ 2 are eigenvalues of M.
As an alternative approach to calculate E, we could bypass F and make use of the moment.We denote this definition E M , In the SI we demonstrate that to the leading order the two definitions are approximately equal.We note that while ( 21) is an area-based calculation, (5) in the proper text is vertex-based.
The coefficient of determination, R 2 , follows the standard definition, Here 'Var' denotes variance, f is the data presented in array form, and f is the corresponding array generated via fitting (such as for |E| 2 ) or a theoretical prediction (such as for a r ).

Monte Carlo simulation
Results shown in Fig. 4e-h are generated via a simple Monte Carlo simulation.We generate deformation from regular hexagons using Eq. ( 2).The initial displacements u 0 i follow independent and identical distributions where θ j is uniformly distributed between [0, 2π], and d i follows We set β to ensure the variances of the components are 1.The exponent p = 4/3 gives test case 1 (red) shown in Fig. 4e, whereas p = 10 gives test case 2 (blue), which has a non-monotonic distribution of d i , and hence ûk,0 .These random initial displacements lead to a centroid-translation (SI), Center-correction of u i,0 provides u i and y i in (2).Once y i 's are generated, other quantities such as vi , ŵi and |E| 2 are computed according to formula provided in the proper text.

Figure 1 :
Figure 1: (a)-(h) Examples of randomized polygonal networks in nature.(a) Land cracks due to desiccation [37]; (b) Ice wedges from northern Canada [61]; (c) Veins on a Ficus lyrata leaf [11]; (d) Desiccationpattern of an ancient lake on Mars[43]; (e) A snapshot of Salar de Uyuni in Bolivia, world's largest salt flat[23]; (f) Supergranulation on solar surface[29]; (g) Plated MDCK cells (this work).(h) A processed image of a developing Drosophila wing disc (this work) with the aspect ratio of cells color-mapped.(i) A regular hexagon P R (blue, with vertices x j ), a deformed hexagon P (red, with vertices y j ), and a uniform deformation as mean-field approximation P (black dashed outline, with vertices y j ).

Figure 2 uses 4
representative cases to demonstrate the agreement.(More are shown in Fig. S2 in the SI.)The left column shows the PDFs of |E| 2 , which are very well fitted by the exponential form, exp(−β|E| 2 ), where β is extracted as a fitting parameter.(See also Fig. 4b where all |E| 2 profiles are presented and the value of β is theoretically predicted.)

Figure 2 :
Figure 2: Universality in strain and aspect ratio distributions.Left column shows the PDFs of |E| 2 , fitted with an exponential form exp(−β|E| 2 ) to extract β.This β value is used in Eq. (10) to generate the theoretical curves in the middle column (dash), in comparison to the aspect ratio data (symbols).The coefficients of determination, R 2 are shown within the panels.Right column: both data and theoretical curves from the center column are normalized using the mean values, and fitted with a k-Gamma function (12) (thick dashed).A single parameter k is extracted and shown in the figure legends.Data are from[21],[2], and[60], respectively, for the top 3 rows; and from this work for Droso.

Figure 3 :
Figure 3: (a)  Overall correlation between a r and T ( β −1 ) for all 99 data sets (symbols); the theoretical prediction (dashed) is generated per(10).The inset alternatively presents E from a moment-based calculation which further improves agreement.(b) A comparison between data and predicted temperature, T n and T ("theory").The inset shows that sub-ensemble temperature T n 's are quantitatively similar to tessellation temperature, T .

Fig. 4a presents ŵ1, 2
Fig.4apresents ŵ1,2 in the hexagon sub-ensemble (n = 6) for all 99 tessellations, whereas the cases for n = 5 and 7 are included in the SI.PDFs are all normalized for comparison with standard Gaussian, N (0, 1) (dark solid lines).Although the PDFs exhibit appreciable fluctuations due to the relatively small sample size in the n-sub-ensemble, the approximate normalities are evident.Quantitative similarities of ŵ1,2 are demonstrated in Fig.S5.In addition, ŵ1,2 are indeed only weakly dependent, as Cov( ŵ1 , ŵ2 )/Var( ŵ1 ) = 0.078 ± 0.046 for all cases, consistent with the anticipated 2 degrees of freedom.All E 2 distributions (normalized by the predicted temperature T n ) are shown in Fig.4b.
[16] and then cultured based on the procedure of Dye et al.[16].Live wing discs were imaged using a Perkin Elmer Ultraview spinning disc confocal microscope every 8 mins for 12 hours.
Plated MDCKIIG cells (MDCK) MDCKIIG (a gift from W. James Nelson, Stanford University) cells were cultured in low-glucose Dubecco's modified Eagle's medium (DMEM) (Life Technologies) supplemented with 10% fetal bovine serum (FBS) and antibiotic-antimycotic.Cells were used at low passage number, checked regularly for contamination by cell morphology and mycoplasma testing.Cells were plated at different densities (1.5, 3, 4.5, 6, and 7.5×10 4 cells/cm 2 ) on coverslips coated with 0.6 mg/ml of collagen for 15 min at room temperature and washed with PBS.After 48 hours, cells were fixed with 4% paraformaldehyde in PBS++ (phosphate-buffered saline supplemented with 100 mM MgCl 2 and 50 mM CaCl 2 ) for 10 min at room temperature.Immunostaining was performed as in Ibar et al. using mouse anti-ZO1 (1:1000, Life Technologies #33-9100) and anti-mouse Alexa Fluor 647 (Jackson ImmunoResearch)