{'fraction_non_alphanumeric': 0.03200589970501475, 'fraction_numerical': 0.014995083579154375, 'mean_word_length': 4.659710628825821, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hierarchical structure of the energy landscape in the Voronoi model of dense tissue

The Voronoi model is a popular tool for studying confluent living tissues. It exhibits an anomalous glassy behavior even at very low temperatures or weak active self-propulsion, and at zero temperature the model exhibits a disordered solid structure with no evidence of a rigidity transition. Here we investigate the properties of the energy landscape in this limit. We find two disordered solid phases that have similar structural features but that differ in the ultrametricity of their energy landscapes; the crossover between these two states shares phenomenological properties with a Gardner transition. We further highlight how the metric used to calculate distances between configurations influences the ability to detect hierarchical arrangements of basins in the energy landscape.

The Voronoi model is a popular tool for studying confluent living tissues. It exhibits an anomalous glassy behavior even at very low temperatures or weak active self-propulsion, and at zero temperature the model exhibits a disordered solid structure with no evidence of a rigidity transition. Here we investigate the properties of the energy landscape in this limit. We find two disordered solid phases that have similar structural features but that differ in the ultrametricity of their energy landscapes; the crossover between these two states shares phenomenological properties with a Gardner transition. We further highlight how the metric used to calculate distances between configurations influences the ability to detect hierarchical arrangements of basins in the energy landscape.
Understanding the collective behavior of cells in biological tissues has become one of the major interdisciplinary challenges of recent years, with applications ranging from wound healing to cancer treatment [1][2][3][4][5][6]. Both experimental and theoretical efforts have been crucial in understanding the properties of these tissues and the mechanisms by which tissue properties are regulated, for example in the way that tissues can transition from rigid to flexible as the properties of individual cells are regulated [3,[7][8][9][10][11][12]. Rigidity transitions are also seen in particulate systems, such as granular materials and colloidal suspensions, in which changes in particle density and temperature can lead to disordered materials in a kinetically arrested jammed or glassy state [13][14][15]. In living tissues, the nature and properties of the rigid states are still under debate [3,[7][8][9], owing both to the explicitly non-equilibrium nature of cellular motion and the many-body interactions found in confluent tissue [16]. These differences raise several challenges to the generalization of ideas and methods developed in the context well-studied particulate matter [4,5,7,17,18].
Several models have been proposed to understand the collective behavior of cellular systems, from single particle descriptions to density field models [19][20][21][22][23][24]. The Voronoi model represents a confluent tissue as a space filling polygonal tiling, where each positional degree of freedom corresponds to a cell whose shape is obtained by an instantaneous Voronoi tessellation [25][26][27]. The dynamics is controlled by an energy functional that is quadratic in the area and perimeter of each cell (described in more detail below), and the mechanical properties of the tissue can be either solid-like or fluid-like depending on the temperature and a shape parameter, p 0 , which quantifies the target shape of the individual cells [25].
At zero temperature (or in the absence of cellular activity) it has been argued that the Voronoi model possesses a finite shear modulus over its entire range of model parameters [28]. This is in sharp contrast with particulate systems, in which a zero-temperature rigidity transition can be observed by changing the system density [29][30][31][32]. The particulate jamming transition is typically interpreted in In the middle is the original minimized configuration (red) and the perturbed one (green). On the right is the original minimized configuration (red) and the one minimized after the perturbation (blue). the context of constraint counting, in which the transition occurs when the number of independent particle-particle contacts equals the number of degrees of freedom [33,34]. As described below, the 2D Voronoi model is always at this point of marginal stability [28], and thus an analysis based only on the balance between constraints and degrees of freedom is insufficient. It has been proposed instead that energetic rigidity, in which not only simple constraints but also residual stresses play a controlling role [17,18,35,36], is a better framework for understanding Voronoi model rigidity in the athermal limit [36].
Here we explore this unusual athermal regime of the Voronoi model and show that, even in the absence of a zero-temperature rigidity transition, there is nevertheless a profound change in the statistics of the energy landscape in different regions of the model parameter space. We find evidence for a transition to an ultrametric, hierarchical arrangement of basins in the energy landscape, suggesting two different phases of a disordered solid [39][40][41][42]. The ultrametric state is characterized by energy minima forming a tree-like structure in phase space where minima within a given sub-basin are much closer to one another than they are to minima in any other sub-basin [29,31]. This is consistent with a Gardner phenomenology [29,31]; the phenomenological , and the subdominant ultrametric constructed from it using a minimum spanning tree [37,38]. Matrices corresponding to N = 4096 and p0 = 3.75, 3.83 are shown, where the different distances are grouped using a single-linkage clustering algorithm which clusters the minima sequentially by distance. All results are averages of 10 initial configurations subject to 100 perturbations and minimizations each.
We model the confluent tissue as a monolayer of N cells [19,20,25] in a square domain of side-length L with periodic boundary conditions. Each cell i is represented by its center r i with a shape given by an instantaneous Voronoi tessellation of the space. We choose the unit of length to be given by the square root of the average area of all the cells. We can then write a dimensionless version of the contribution of each cell to the energy functional as [17,21,22,45] Here a i and p i are the dimensionless area and perimeter of cell i, p 0 is the target perimeter, and k A represents the ratio between the relative stiffness of the area and perimeter elasticity of the cell. Biologically, the first term models cellular incompressibility and the resistance of the cellular monolayer to height fluctuations; the second term models the competition between active contractility of the actomyosin subcellular cortex and the effective cell membrane tension due to cell-cell adhesion and cortical tension.
For each set of parameters, we start with N cells distributed at random positions. The configuration is then minimized using the FIRE algorithm [46,47], which we halt when the maximum net force on each cell is less than 10 −12 . Due to numerical constraints, we consider p 0 3.85, since it has been shown for athermal systems that, above this value, configurations with multi-fold vertices are obtained which lead to numerical instabilities in the minimization protocol [28]. Further details of the simulations can be found in the Supplemental Material.
To probe the structure of the energy landscape, we start from an initial configuration that corresponds to a local minimum and perturb it to find new stable configurations. In our primary perturbation protocol, we displace the position of each cell according to the vector is the perturbation to cell i along the x-axis, X 2i+1 = ε sin(θ i ) is that along the y-axis, and θ i is a random angle uniformly distributed between 0 and 2π. We considered a random length ε drawn from a uniform distribution between 0 and ε max , to guarantee the possibility of visiting minima in the same and different top-level basins. The norm of the perturbation vector is After the perturbation, we subtract the global translation of the tissue, and then let the tissue relax to a new minimum. An example of this process is shown in Fig. 1; details of alternate "perturband-minimize" schemes can be found in the Supplemental Material, where we show that our results are not qualitatively sensitive to these details.
We will be exploring different metrics to characterize distances between minima in the energy landscape. To begin quantifying these distances we consider the contact metric (denoted by the superscripted X) discussed in Refs. [40][41][42], where d X (a, b) is the distance between configuration a and b, and C a ij is the 2D contact vector between two cells, where each component C a ij,x = x a i − x a j is the distance along the respective axis between cells i and j if those cells share and edge, and C a ij = 0 otherwise. In the Supplemental Material, we show how the normalized distance to the original minimum, d X N (a, b) = d X (a, b)/ |a||b|, where the norm of a configuration corresponds to |a| = d X (a, 0), scales with the number of cells and the norm of the perturbation vector. We observe that d X (a, b) ∼ √ N and choose ε max = 0.5 √ N for all p 0 , which is large enough so that the perturbed configuration does not always relax to the initial minimum but small enough that nearby minima are accessible. Figure 2 depicts matrices where the color corresponds to the distance between minima, for all pairs of minima found for p 0 = 3.75 and 3.83. These matrices were constructed for 10 3 minima obtained for a tissue of 4096 cells. Each element of the matrix corresponds to the distance between two minima, a and b, given by Eq. (2). The distances are all sorted using the single-linkage clustering algorithm on the metric of the fully minimized systems, which groups them sequentially based on their relative distance [48]. The colors represent different distances. In white are the minima that are closest to each other. We find groups of minima that are all at this minimum distance, forming a white region that corresponds to subbasins. The matrices do not show any substantial visual change with p 0 .
A Gardner phase is characterized by an ultrametric phase space consisting of a tree-like structure, where minima within a given sub-basin are all much closer to one another than to minima in any other sub-basin. This property is codified by an ultrametric inequality, where a, b and c are three different configurations in phase space and d X (a, b) is the distance between configurations.
To verify if the properties of the tissue are consistent with a Gardner phase, we compute how close the metric is to being ultrametric. To do so, we first find the subdominant ultrametric, d < (a, b), i.e., the ultrametric that is closest to d X (a, b) itself. The subdominant ultrametric can be found by first computing the distances in the minimum spanning tree of the space of minima. Then, for each pair of minima a and b, we compute the path between them in the minimum spanning tree and define d < (a, b) as the largest distance between two neighboring minima along the path [38]. The corresponding matrices are shown in Fig. 2. Having found d < (a, b), we finally calculate the generalized distance between the metric and the subdominant ultrametric using where · denotes the average over all configuration pairs a and b. If D X = 0 then the energy landscape is ultrametric, while D X > 0 quantifies how far it is from ultrametricity. In the inset of Fig. 2, we show that D X increases slightly with p 0 , but more importantly it depends strongly on N . In the main plot we re-scale D X / √ N and obtain a reasonable collapse of the data. Since the typical distance between minima and the distance to ultrametricity both scale with √ N , the space with this contact metric is not ultrametric in the thermodynamic limit [42].
Using distances based on the contact vectors suggests that the landscape of the Voronoi model is not ultrametric, but does the choice of metric itself influence this result? We note that in the Voronoi model, the contact network of the tissue does not change significantly even as the tissue rigidity changes substantially [18,36]. We further note that the energy functional in Eq. (1) is a simple collection of harmonic springs, in a coordinate basis of shape space rather than in the positional basis of the degrees of freedom generating the shapes. This suggests a different metric might be more appropriate, and in this context we propose one based on the contribution of each cell i to the total energy of the tissue, E i . We take the same form for the metric as Eq. (2), but where C a ij → C a ij = E a i − E a j , if i and j are neighbors and zero otherwise. We call this metric the "energy metric", d E (a, b). We adopt the same perturb-and-minimize protocol as before, using ε max = 0.1 √ N for all p 0 . As such, we keep biasing the perturbations to nearby minima in the new metric, otherwise more simulations would be needed to probe the same volume of configuration space. In the Supplemental Material we show that the results using the contact metric remain qualitatively the same using the new ε max . Just like the contact metric, the energy metric scales with system size, d E (a, b) ∼ √ N , since it depends on the total number of cell-cell contacts. We observe in the inset of Fig. 3 that, for the energy metric, the distance to ultrametricity (D E ) does not scale with N and decreases with p 0 (inset of Fig. 3). Since the distance between minima scales as d E (a, b) ∼ √ N , while the generalized distance does not depend on the system size, this suggests that the system does, in fact, become ultrametric in the thermodynamic limit: In the case of the contact metric the calculated values are already normalized since we increase the box size with N , while the typical cell size is fixed. In the case of the energy metric this is no longer the case. Thus, to properly compare the system properties at different p 0 (since E i varies with p 0 ), we also consider a normalized version of the energy metric: d E N (a, b) = d E (a, b)/ |a||b|, for which the typical distance between configurations does not depend on either p 0 or N . Figure 3 shows the schematic representation of the normalized energy metric, d E N (a, b), and its subdominant ultrametric. From this representation, we can already observe changes in the structure of the energy landscape. As p 0 increases, the color gradient is less smooth and the boxes corresponding to the different sub-basins become sharper. Furthermore, it is also observed that more minima fall into the same sub-basin. These properties suggest that, when the value of p 0 decreases, the structure of the energy landscape becomes more hierarchical [42].
We also compute the generalized distance to ultrametricity when using d E N . Again we see that ultrametricity is approached with increasing system size. Furthermore, with the scaling shown in the main panel of Fig. 3 we can collapse the different curves. The values of p * 0 and ν were chosen such that we obtain the best collapse. Recent work in particulate systems has interpreted a similar scaling as a distance to jamming [30,42]. This suggests not only that the Voronoi model has an ultrametric structure in the thermodynamic limit, but that there may be a transition between two different solid phases. Previously, it was shown that the athermal Voronoi model did not have a rigidity transition [28]. Nevertheless, we show that there is a clear difference between the energy landscape for low and high p 0 : at low p 0 the glass state is characterized by large residual stresses and an ultrametric energy landscape, and at high p 0 the energy landscape is not ultrametric. We find that the change in the structure of the energy landscape occurs for preferred shape parameters in the range p 0 = 3.75 − 3.83, which is close to where the zero-temperature shear modulus changes markedly [28] and where the dynamics at finite temperature changes in character [25,27,49].
As p 0 increases, fewer sub-basins are found inside each basin (as represented by the different unconnected clusters in Fig. 3), suggesting that the energy landscape flattens out. Another way of exploring this flattening of the energy landscape is by studying the behavior of the model as the relative area modulus k A is varied. In the limit k A = 0 the Voronoi model is no longer marginally constrained, and it acquires a residual-stress-based rigidity transition as a function of p 0 at T = 0 [28]. As shown in Fig. 4, for p 0 < 3.79 we find that the tissue is both rigid and the energy landscape is ultrametric. For slightly larger p 0 the energy landscape deviates from ultrametricity and the variance of D N E increases significantly. In this regime, the energy landscape consists of a mixture of hierarchical basins and several nearly flat basins and so, in many cases, small perturbations will not drive the system to a different minimum. Due to finite size effects it is difficult to assess if a new solid phase exists. In the Supplemental Material we use a simple technique to estimate the transition from the solid to the fluid phase. Using the fraction of configurations with zero energy we estimate a transition point around p * 0 = 3.8022 ± 0.0001, while in Fig. 4 p * 0 = 3.798±0.001 seems to be more appropriate value to collapse the data. We used ν = 0.4 ± 0.01, which was the best value for the collapse. More simula-tions would be needed to conclude whether a new solid phase at k A = 0 exists before the fluid phase, or if the observations are finite size effects and both transitions actually coincide. For p 0 > 3.8 the energy landscape is flat, characteristic of a fluid-like tissue. For any k A > 0, different energy minima are found for all p 0 , consistent with previous work showing a finite shear modulus for the whole range of model parameters investigated [28]. Although we have not explored the thermal case, recent studies with the thermal 2D Voronoi model also suggest a change in the energy landscape close to p 0 ≈ 3.81, where there is an emergence of a fractal-like energy landscape and cells become virtually free to diffuse in specific phase space directions up to a small distance [49].
In summary, we have found indications of a hierarchical structure of the energy landscape in a model of dense biological tissue whose zero-temperature rigidity is quite different from that of constraint-based particulate systems. Strikingly, we find that the choice of metric to characterize distances between minima is crucial: defining distance based on changes in neighboring contact vectors vs contributions to the energy give qualitatively different interpretations of the structure of the energy landscape. In particulate systems the contact vectors enter explicitly in the relevant energy functional -i.e., the energy of a soft harmonic repulsion or a Lennard-Jones interaction is a simple function of the contact vector between interacting particles. In Voronoi models the energy cannot be decomposed into independent pairwise contributions, which we speculate is the reason that choosing a distance metric based on the total energy associated with each degree of freedom is required to uncover the hierarchical structure of the landscape. We further speculate that this may point more generally to the importance of the choice of metric for systems in which many-body interactions dominate over pairwise ones. An avenue for future research could be relating these tissue-like systems to particulate ones, such as soft or hard spheres [40,42]. This could be done by establishing the relation between the effects of p 0 in the Voronoi model and pressure in particulate systems. Since both exhibit an ultrametric landscape, this could allow a generalization of glassy physics outside of particulate systems and glass-forming materials [14,15,29].

I. ACKNOWLEDGMENTS
The authors acknowledge financial support from the Portuguese Foundation for Science and Technology (FCT) under Contracts no.