Inferring the flow properties of epithelial tissues from their geometry

Amorphous materials exhibit complex material properties with strongly nonlinear behaviors. Below a yield stress they behave as plastic solids, while they start to yield above a critical stress Σc. A key quantity controlling plasticity which is, however, hard to measure is the density P(x) of weak spots, where x is the additional stress required for local plastic failure. In the thermodynamic limit P(x) ∼ x θ is singular at x = 0 in the solid phase below the yield stress Σc. This singularity is related to the presence of system spanning avalanches of plastic events. Here we address the question if the density of weak spots and the flow properties of a material can be determined from the geometry of an amorphous structure alone. We show that a vertex model for cell packings in tissues exhibits the phenomenology of plastic amorphous systems. As the yield stress is approached from above, the strain rate vanishes and the avalanches size S and their duration τ diverge. We then show that in general, in materials where the energy functional depends on topology, the value x is proportional to the length L of a bond that vanishes in a plastic event. For this class of models P(x) is therefore readily measurable from geometry alone. Applying this approach to a quantification of the cell packing geometry in the developing wing epithelium of the fruit fly, we find that in this tissue P(L) exhibits a power law with exponents similar to those found numerically for a vertex model in its solid phase. This suggests that this tissue exhibits plasticity and non-linear material properties that emerge from collective cell behaviors and that these material properties govern developmental processes. Our approach based on the relation between topology and energetics suggests a new route to outstanding questions associated with the yielding transition.


Introduction
A fascinating aspect of biological systems is their ability to grow into well-defined shapes [1]. Thinking about tissues as materials, what should their properties be to allow for robust morphogenesis? One view is that tissues are viscoelastic fluids, molded into desired shapes by surface tension and active forces [2][3][4][5][6][7][8][9]. An alternative picture is that they are yield stress materials [10] similar to clay. Such materials allow for great control, since shape is changed only if the magnitude of shear stress Σ is above the threshold yield stress Σ c . These approaches can be thought of as two extremes of a continuous spectrum of models, since at finite temperature, or at finite level of active stress and cell divisions in biological systems [3,11], materials always eventually flow. Experimental evidence of glassy behavior [12][13][14] indeed suggests relevance of intermediate cases. Quantitatively, an interesting observable to distinguish these regimes is the ratio between the strain rate increment δγ and the stress increment δΣ causing it. This ratio is simplyγ/Σ for a Newtonian liquid,

Flow and loading curves of vertex model
We use the standard vertex model of epithelial tissues [37,41] where the 2D network of polygonal cells is assigned an energy function where A c , L b and P b are cell area, bond length and cell perimeter, respectively. Model parameters specify preferred cell area A 0,c , bond tensions Λ b and cell perimeter stiffness Γ c . To avoid localisation of flow in a narrow shear band, which occurs in the homogeneous vertex model [42], we introduce cell size polydispersity for shear flows, see appendices A and C. In all simulations the network is in solid phase with the normalised preferred perimeter p 0 ≡ −Λ b /(2Γ c √ A 0,c ) 3.41 for all cells, well below the rigidity transition point p * 0 3.81 [38]. Note that our results below may not hold in the fluid phase of the vertex model [43], or in active tension networks with isogonal soft modes [44].
The dynamics of the cellular network is described by overdamped dynamics of vertex positions r α where ν = 1 is a mobility. A T 1 transition occurs when a bond length becomes smaller than a threshold length T 1 , see appendix B. We use two ensembles of isotropic disordered networks with N = 400 and N = 2500 cells, as described in appendix D. An example of a network with N = 400 cells is shown in figure 1(a) left. We perform simple shear strain simulations on these networks at a constant strain rate, illustrated in figure 1(a) right foṙ γ = 10 −4 . The network initially responds elastically: the shear stress Σ in the network grows almost linearly ( figure 1(b)). As the strain is increased, T 1 transitions occur and relax the stress in the network, visible as sudden drops in the stress vs strain curve. In figure 1(a) right we visualise recent T 1 transitions that occurred during a strain increment Δγ = 0.15 by coloring participating cells in red. T 1 transitions appear to be correlated and organised into avalanches of various sizes, corresponding to widely distributed stress drops in figure 1(c). Eventually a steady state is reached in which stress relaxation due to T 1 transitions balances the elastic loading. In figure 1(d) we show the steady state flow curve. It is well described by the Herschel-Bulkley lawγ ∼ (Σ − Σ c ) β [27] with the yield stress Σ c ≈ 0.93 7 and exponent β ≈ 1.3 8 .

Energy cusp at T 1 transitions
Next we characterize the T 1 transitions or elementary plastic events. For this purpose, we identify the bond in a network that will first disappear under strain. We then constrain the length of that bond to a value L * and determine the energy of the network under strain. The original bond length is assigned negative values and the new bond that appears through the T 1 transition is assigned positive values. In figure 2 we show the energy of the network relative to the energy of the unconstrained network, as a function of strain ΔE(L * ; γ) ≡ E(L * ; γ) − E(γ). Originally the system is in a metastable state, corresponding to the local minimum of ΔE. As the shear stress increases the minimum disappears and the T 1 transition occurs. The energy profile shows a cusp at the onset of T 1 where L * = 0, a well known feature of the vertex model energy landscape [38,[46][47][48]. The presence of a cusp in the energy profile allows us to relate the bond length of short bonds to the additional force 9 x at the bond needed to drive a T 1 transition. Namely, expanding ΔE(L * , γ) in bond length around the equilibrium value L reads: and we see that at T 1 transition, corresponding to L * = 0, the energy barrier to the T 1 transition is: Thus, the force on that bond required to trigger the T 1 transition is: Since the effective stiffness of the bond ΔE (L) is expected to be finite at a T 1 transition, we find E b ∼ L 2 and x ∼ L. These relationships were previously reported and supported analytically in a continuum elastic model of a T 1 transition [48]. Here, we find that they are a generic consequence of the cusp in the energy profile 10 . Ultimately, the origin of this cusp lies in the form of the energy function that depends on the bond lengths, cell perimeters and area. These quantities are smooth functions of the vertex positions and the energy function is therefore continuous even when topology changes. However, at a T 1 transition a bond is removed and another one added. This implies that gradients of the energy function are discontinuous at the T 1 transition, leading to a cusp in the energy profile. Consequently, forces can change discontinuously at the transition point. Therefore, we expect to generically find x ∼ L in cellular systems such as epithelial tissues and dry foams for which the dependence of the energy on the vertex position is topology-dependent. By contrast, particle systems in which the energy depends on the particle positions independently of any notion of topology cannot show such a cusp (as long as the interaction potential is smooth). Furthermore, due to the cusp at the T 1 transition the stiffness of the corresponding displacement mode does not vanish, as it would at the plastic event in particle systems 11 . Therefore, we do not expect to 9 Note that exerting a stress increment at the boundary of the system will in general generate a bond force proportional to that increment, so the quantity x we use here characterizes well the distance to a local yield stress. 10 For example, if the energy landscape of vertex model were smooth at T 1 transitions the relationship would be x ∼ √ L and E b ∼ L 3 . 11 In systems with smooth energy function a plastic event corresponds to the usual saddle-node bifurcation. . We find that at the lowest strain rates the effective exponent θ eff. converges to θ ≈ 0.56-0.77, see appendix J for details. As in isotropic networks, at low L r we find P(L r ) ∼ const., as expected for any finite strain rate [33]. find a signature of local plastic events in the eigenvalues of Hessian matrix of energy function, which could explain the lack of soft non-localised modes recently observed in the Voronoi vertex model 12 [49].
We test these predictions in isotropic disordered networks by forcing bond length to attain a very small value L min = 10 −6 and determining the corresponding network energy change and the constraining force magnitude x. Results shown in figures 2(c) and (d) are consistent with our predictions. Therefore, identifying the locations of short bonds allows us to read the map of 'weak spots' in the system, as well as to deduce the distribution P(x).

Stability of the cellular network
After a T 1 transition, the network relaxes to a new metastable state with redistributed shear stresses. We measure the stress redistribution by enforcing a T 1 transition and measuring the change in stress within each cell after the network had relaxed (the cellular stress is defined as in [50]). In an elastic 2D medium we expect the shear stress redistribution to be consistent with that of a force dipole in an elastic medium: Δσ xx ∼ cos(4ϕ)/r 2 , Δσ xy ∼ sin(4ϕ)/r 2 [23]. Figure 3(a) shows the shear stress redistribution, obtained by orienting the disappearing bond direction along the x-axis and averaging over 50 realisations, as detailed in appendix E. We find a clear four-fold symmetry of Δσ xx (inset) as well as inverse quadratic decay of its magnitude, as expected for a force dipole and consistent with simulations of 2D foams [40].
The stress change after a T 1 transition can trigger new T 1 transitions if there are bonds with small x in the network. In the solid phase, the stability of the network with respect to extensive avalanches of T 1 transition imposes P(x) ∼ x θ with θ > 0, otherwise never-ending avalanches would occur [31]. As we have demonstrated x ∼ L in the vertex model. Therefore, bond length distribution should vanish with the same exponent P(L) ∼ L θ , and θ can be extracted from P(L). We measure the cumulative distribution C(L)≡ L 0 P(L )dL in disordered isotropic networks. We find a scaling regime, whose range of validity grows with system size, for which θ ≈ 0.5-0.6 ( figure 3(b)), see appendix H for details. At even smaller L, 12 The voronoi vertex model has the same energy function as the usual vertex model, but its degrees of freedom correspond to cell centers. The network topology is constructed from these centers by performing a Voronoi tessellation at each time-point. the bond lengths distribution departs from this scaling as P(L) ∼ const., as expected due to finite size effects and also observed in elasto-plastic models [28].
Interestingly, these values are consistent with those found in two-dimensional elasto-plastic models [28]. In these coarse-grained models, the material is described as a collection of mesoscopic blocks with a simplified description of plastic events: a block yields when the local yield stress is reached, it accumulates plastic strain and redistributes stress in the material as a force dipole [51,52]. Since both ingredients are present in vertex model as well, as we have seen, it is not surprising that we find a consistent value of θ.
We next studied the evolution of the exponent θ during the transient loading period, between the initially isotropic network and the steady state. Surprisingly, it has been predicted that θ would then non-monotonically depend on strain [33], a result observed in elasto-plastic models [26] but only indirectly observable in amorphous solids where P(x) is very hard to access [34][35][36]. To test directly this prediction, we measure the bond length distribution as a function of strain at a small constant strain rateγ = 10 −4 , close to the quasi-static limit. Note that even in the thermodynamic limit we expect to find singular P(x) and P(L) only in the quasi-static limit of vanishing strain rate (at any finite rate, there are always T 1 transition occurring leading to θ = 0 [33]. However, in a finite system we can still measure an effective exponent θ eff. . We confirm that the evolution of θ eff. with strain is non-monotonic, see figure 3(c).
It is also important to quantify the distribution of bond lengths in steady state flow, see figure 3(d). At high strain rates we find θ eff. → 0 as expected, while in the limit of vanishing strain rates the effective exponent approaches the value θ ≈ 0.56-0.77 (figure 3(d)). Thus, P(L) can be used to locate the distance to the yield stress, at least in this setting where noise is absent.
Finally, note that the fact that the density of bonds vanishes at small argument as P(L) ∼ L θ is a necessary constraint required by network stability [31]. Such a power-law distribution of weak regions is found in amorphous materials under athermal quasi-static shear and in slowly driven glassy materials in various contexts [32]. Generating glasses in which weak regions are depleted even further requires extreme annealing [53]. Cellular networks could also be made more stable, for example by sequentially performing T 1 transitions on the shortest bonds as discussed in reference [54]. Such an annealing process depletes the distribution of short bonds and can create a gap in the distribution P(L) that does not occur in the distributions discussed here. More generally, we expect our results to be relevant for describing disordered developing tissues undergoing shape changes through T 1 transitions, whereas it might not describe very well-annealed tissues.

Collective behavior of T 1 transitions
Quite generally in disordered systems [31], a singularity in the density of weak regions P(x) is synonymous to avalanche-type response where many weak regions-here T 1 's-rearrange in concert. To test if this idea holds in the vertex model, we quantify the correlation between T 1 transitions by using a susceptibility motivated by the four-point susceptibility χ 4 studied in glasses [55,56], elasto-plastic models [57,58] and vertex models [49]. We defineχ 4 as a normalised variance of the number n T 1 (τ ) of T 1 transitions in the time-window τ [57,59]: If T 1 transitions were completely independent, the variance would be equal to the mean at any τ and χ 4 (τ ) = 1. For T 1 transitions organised in avalanchesχ 4 grows and reaches a maximum at τ A , corresponding to the typical avalanche duration. The value at the maximum can be interpreted as a characteristic avalanche size S ≡χ 4 (τ A ). After an avalanche, the stress has relaxed locally and new avalanches are less likely to occur. This effect leads to a decay ofχ 4 (τ ) for τ > τ A . We measureχ 4 (τ ) in steady state shear flow at different strain rates, see figures 4(a) and (b). We find thatχ 4 displays a peak which grows near the transition pointγ → 0, indeed supporting the idea that the dynamics becomes collective at that point. Furthermore, we find that theχ 4 curves at different strain rates collapse when re-scaling the axes asγ 2/3 τ andγ 1/4χ 4 . Therefore, as the strain rate vanishes, the mean avalanche size diverges as S ∼γ −1/4 and the mean avalanche duration as τ A ∼γ −2/3 .

Fly wing epithelia
We have shown that in the vertex model, the bond length distribution P(L) is indicative of the regime in which the material flows: it presents a singular distribution approaching the solid phase, where the dynamics becomes collective and the flow curve is non-linear.
As a first test of the relevance of these ideas to real tissues, we analyze the bond length distribution in wing epithelium of the fruit fly at two stages of development: (i) during pupal wing morphogenesis, imaged in vivo [6] and (ii) the wing disc epithelium, in third instar larva wing disc epithelium imaged ex vivo [60]. In the pupal wing we considered a region defined by the longitudinal veins denoted L4 and L5, and the posterior crossvein (yellow cells in figure 5(a)), imaged at 5 min intervals between 19 and 23 h after puparium formation, collected from 3 experiments [61]. In figure 5(b) we show the analyzed region in the wing disc epithelium, corresponding to the wing disc pouch, imaged at 5 min intervals over about 13 h and collected from 5 experiments. In both figures 5(a) and (b) we indicate in red the cells that have lost or gained a bond as a part of T 1 transition in the last 5 min (the time resolution of experiments).
Note that the wing disc epithelia have been developing for about 100 h before the imaging has started [60], while pupal wings have undergone a significant three-dimensional shape change, called eversion [62], before the pupal morphogenesis. Therefore, the initial state of these tissues at the beginning of the experiments could already contain a significant strain history. Thus, although the amount of strain accumulated during experiments is small (∼ 0.1-0.2 for pupal wings [6,61] and ∼ 0.05-0.1 for wing discs [60]), it is not clear if we are in a large or small strain regime, with respect to the strain where figure 3(c) displays a minimum.
Remarkably, we find that in both tissues the bond length distribution P(L) vanishes at small L with the effective exponent θ eff. ≈ 0.7-0.9 (see figure 5(c) and appendix K for details), similar to those measured in the slowly flowing vertex model at small or large strain. This observation suggests that the scaling relation x ∼ L holds in real tissues as well, in the range of L that we can probe. Clearly for very small bonds, this relationship must eventually break down due to finite size of vertices [63,64].
It would be very interesting to test this scaling relation directly by perturbing the system. It could be achieved by observing tissue response to a localised mechanical perturbation, such as laser ablation of a single bond: if the energy landscape were smooth with no cusps, the stiffness of the corresponding displacement mode would vanish in approach to a T 1 transition as expected near a saddle node bifurcation. As a consequence, a strong locally heterogeneous response would be observed in the experiment at locations of short bonds just before they rearrange, as observed in particulate amorphous solids [65] preceding a plastic event. On the other hand, in presence of a cusp there would be no softening and no strong locally heterogeneous displacement preceding cell rearrangements.
Our observation also raises the possibility that developing tissues can lie in a non-linear regime with δγ/δΣ γ/Σ, where collective effects are important. Unfortunately, such collective effects are very hard to measure in our experimental data, because the strain rate is not stationary, leading to difficulties using the definition ofχ 4 . Thus, it would be important to look for non-linear effects more directly by measuring stress dynamics using laser ablation experiments and comparing them to elastic and plastic flow components [6,[66][67][68]. Alternatively, epithelia obtained as cultured cell monolayers might provide interesting experimental systems that allow for direct rheological experiments. In such systems non-linear flow properties have been observed [69] and bond length distributions P(L) could be measured in different flow regimes.

Discussion
We have shown that in the vertex model in its solid phase, the distribution of bond lengths provides the distribution of local distances to yield stress x. The distribution P(x) reveals properties of the regime in which flow is occurring. This result has two consequences.
First in developing epithelia, we observe a singular distribution of bond lengths, consistent with the one found in the vertex model. This result raises the intriguing possibility that non-linear collective effects may be important in tissues, and suggests further empirical tests. Yet, active forces and cell divisions or extrusions could affect the relationship between geometry and flow properties at shear rates comparable or smaller than cell division and extrusion rates [3,11]. Extending our results to this flow regime will require incorporating these active processes. From a theoretical perspective, how the density of weak spots depends on stress in the presence of noise-even a simple thermal noise-is not well understood in amorphous solids, and is just starting to be investigated [11,70]. In this light, it would be important to study in the future how different kinds of noise affect the flow curve and the distribution P(L) in the vertex model.
Secondly, despite the fact that the energy functional of the vertex model is more complex than that of usual particulate materials (in which interaction can be radially symmetric), the relationship between geometry and the presence of weak spots is much simpler in the vertex model. Outstanding questions in the context of amorphous materials, such as predicting how amorphous solids break by forming shear bands in which most plastic events occur, are hampered by the difficulty of measuring the distribution of weak spot with small x [17][18][19][20][21]. The vertex model may thus be ideal to understand the universal aspects by which amorphous materials break and flow.
• ν = 1 Any change of parameters in a particular simulation is explicitly listed below.

Appendix B. T 1 transition implementation
When a bond length becomes smaller than a threshold value T 1 a T 1 transitions is attempted: old bond and corresponding vertices are destroyed and new ones are created, then forces on the new bond are computed and T 1 transition is allowed if the tension in the new bond is positive (forces are stretching the new). Otherwise, the T 1 transition is canceled and the network is reverted to the original state. The choice of T 1 is specified for particular simulations below. To avoid the possibility of an extrusion we do not allow bond loss by T 1 transition for cells with 3 neighbors. any vertex was below F = 10 −3 ). Once the relaxation is over, the magnitude of the force acting on the vertices of the constrained bond is recorded as x, as well as network energy change ΔE.
where (x α , y α ) are coordinates of vertex α, Δγ is strain increment, and L y simulation box size in y direction. The applied strain rateγ = Δγ/Δt was always constant during a simulation. T 1 transition threshold was T 1 = 10 −2 and simulation time step Δt = 10 −3 .

Appendix H. Isotropic θ
We determined the range of values corresponding to this exponent by fitting cumulative bond length distribution obtained from 50 isotropic networks of size N = 2500, which exhibit a broader range of scaling than N = 400 networks. We performed the power law fit on a range of data [l, 0.2] with varying lower limit l as shown in figure 6(a). We find that for values of lower limit l in the range [0.01, 0.05] the exponent θ varies between 0.5 and 0.60.