The role of non-affine deformations in the elastic behavior of the cellular vertex model

Abstract

The vertex model of epithelia describes the apical surface of a tissue as a tiling of polygonal cells, with a mechanical energy governed by deviations in cell shape from preferred, or target, area, A₀, and perimeter, P₀. The model exhibits a rigidity transition driven by geometric incompatibility as tuned by the target shape index, . For with p_*(6) the perimeter of a regular hexagon of unit area, a cell can simultaneously attain both the preferred area and preferred perimeter. As a result, the tissue is in a mechanically soft compatible state, with zero shear and Young's moduli. For p₀ < p_*(6), it is geometrically impossible for any cell to realize the preferred area and perimeter simultaneously, and the tissue is in an incompatible rigid solid state. Using a mean-field approach, we present a complete analytical calculation of the linear elastic moduli of an ordered vertex model. We analyze a relaxation step that includes non-affine deformations, leading to a softer response than previously reported. The origin of the vanishing shear and Young's moduli in the compatible state is the presence of zero-energy deformations of cell shape. The bulk modulus exhibits a jump discontinuity at the transition and can be lower in the rigid state than in the fluid-like state. The Poisson's ratio can become negative which lowers the bulk and Young's moduli. Our work provides a unified treatment of linear elasticity for the vertex model and demonstrates that this linear response is protocol-dependent.

---

Many biological processes, such as morphogenesis,^1–4 wound healing,^5–8 and cancer metastasis,^9,10 require coordinated motion and shape changes of many cells. An important open question in biology is how the large scale mechanics of biological tissue emerges from the properties of individual cells, which are in turn governed by force-generating proteins within the cytoskeleton and adhesion molecules between cells.^11–13 Many theoretical models have been proposed to describe dense epithelia, single layers of very tightly packed cells.^14–18 Among these, vertex models, originally developed from models of soap films,¹⁹ have proven a powerful starting point for capturing the mechanical properties of epithelia. The vertex model describes the apical surface of a confluent tissue as a polygonal tiling of the plane (Fig. 1a).^18–23 Each polygon represents a cell, each edge a cell–cell junction, and each vertex a multicellular junction. Each cell's mechanics are controlled by multiple bio-mechanical processes that were proposed to be effectively described by a mechanical energy determined by deviations of their area and perimeter from preferred values. These preferred values encode bio-mechanical properties such as cadherin molecules concentration, apical ring-contractiliy and more. Force balance via energy minimization then determines the position of the vertices and thus the shape of cells in the tissue. Topological rearrangements resulting in cell intercalation, cell division and motility have also been incorporated, and the model has been highly successful in capturing a range of biological processes, such as tissue growth,²⁴ wound healing,⁷ and tissue organization.²¹

Fig. 1 The vertex model for epithelia. (a) The apical surface of an epithelium is modeled by a polygonal tiling, with each polygon representing a cell. (b) Vertex model phase diagram in the p₀, r plane. Within the blue region cells are unstable and collapse. Within the green region the tissue is in an incompatible state, with neither preferred perimeter nor area achieved, and the ground state is a regular hexagonal lattice. The tissue acts like a solid in response to shear. Within the yellow region both preferred perimeter and area are achieved and cells have a degenerate ground state and the tissue has zero shear modulus. The cell shapes show example energy minima, where a cell may elongate, increase its pointiness using the angle ϕ, or increasing its shear tilt angle θ as described in ref. 25, in order to increase its perimeter while maintaining unit area.

It has been shown that vertex models exhibit a transition between a fluid-like state and a solid-like state where cells are jammed and unable to rearrange (Fig. 1b). This rigidity transition occurs at constant cell density and is driven by both active processes, such as fluctuations in cell-edge tension and cell motility,^{21–23,26–30} as well as by geometric constraints.^25,31 Recent work by us and others has shown that even in the absence of fluctuations and topological rearrangements, vertex models exhibit a rigidity transition associated with geometrical frustration.^25,31,32 In the rigid or incompatible state cells are unable to achieve the target values of area and perimeter and the system is under finite tension, with a unique gapped ground state. In the soft or compatible state, cells achieve both target area and perimeter and the ground state has zero energy. Due to the underconstrained nature of the vertex model, however, the liquid ground state is degenerate as for a given n-sided polygon there are many shapes that preserve area and perimeter. This allows the system to accommodate small shear deformations by finding a new zero energy shape, resulting in vanishing shear modulus.

The linear elastic response of an ordered hexagonal vertex model to external deformations has been examined through calculations of the shear and bulk moduli.^22,25,33 Staple et al.^22,25 evaluated the elastic moduli and first demonstrated the vanishing of the shear modulus in the compatible state. More recently, we showed that the vanishing of both shear and Young moduli in the soft regime stems from the degeneracy of the compatible ground states, which allows the deformed tissue to spontaneously shear to a new compatible ground state to accommodate the external deformation.²⁵ We additionally discovered that the response is highly singular at the critical point, with breakdown of linear elasticity and anomalous coupling between compression and shear, as quantified by the development of a new elastic constant.²⁵

The above studies only allow for affine deformations of the cells. This approximation can be viewed as appropriate for determining the short-time response of the vertex model to strain. The vertex model has, however, additional degrees of freedom and can relax stress by moving vertices in a non-affine way. Murisic et al.³³ incorporated these effects by considering the hexagonal lattice as the union of two sub-lattices with a microscopic shift between them and found that the shear modulus is 2/3 softer than previously reported. Tong et al.²⁸ used simulations to measure the shear storage modulus and viscosity in both ordered and disordered model tissues.

In this paper, we expand upon previous work by incorporating simple non-affine deformations. Using a mean field model for a hexagonal lattice, we derive analytic expressions for all the linear elastic moduli of the tissue, and verify these results using simulations. We show that, away from the critical point, the elastic constants of a regular VM satisfy the standard relations of two-dimensional elasticity of isotropic solids. Despite this Hookean relationship, the mechanical linear response exhibits robust non-affine contributions that can significantly reduce the elastic constants, as known to happen in amorphous solids.^34–37 For instance, the bulk modulus can be softer in the rigid state than in the soft fluid-like state and jumps discontinuously across the solid to fluid transition. We highlight several novel behaviors of vertex model elasticity, such as negative Poisson's ratio and a softening of the tissue as the ratio of area to perimeter stiffness increases. We verify our analytical results using numerical vertex model simulations of a regular tissue.

The remainder of the paper is organized as follows. In Section 1 we state the vertex model simulation and deformation protocol to extract various elastic constants. In Section 2 we introduce the VM and its mean-field implementation used in the present work, and present a new derivation of the ground states that allows us to quantify the degeneracy of the compatible regime. In Section 3, after highlighting the distinction between the affine and non-affine deformations allowed in our model, we present results for all the elastic constants. We conclude in Section 4 with a brief discussion.

1 Vertex model: simulation and deformation protocol

1.1 The vertex model of epithelia

The vertex model describes cells in a confluent tissue as polygons of area A_α and perimeter P_α (Fig. 1a). The tissue energy is written aswhere α labels individual cells and 〈ij〉 indexes edges connecting vertices i and j. The first term embodies the energy cost of cell area deformations, with K the area elasticity and A_α0 the preferred or target area. The second term represents active contractility and elasticity of the cytoskeleton, with Γ the contractility. The third term represents interfacial energy between neighboring cells, with L_ij the length of edge ij and Λ_ij the associated tension controlled by the interplay of cell–cell adhesion and cortex contractility. The tension can become negative when adhesion overcomes contractile surface forces.

The mechanical force on vertex i with position x_i is given by . The tissue rearranges vertices to locally minimize the energy. This can be described quasi-statically by requiring force balance at each time-step, or dynamically by assuming that vertices relax according to overdamped dynamics where viscous drag balances the mechanical forces: with γ a friction coefficient.

As the network relaxes, edges may shorten and cells may shrink, resulting in topological rearrangements that reconfigure the network. In T1 transitions, also known as cell–cell intercalations, a junction between two cells shrinks to a point and a new edge is formed, causing two originally neighboring cells to lose contact and two previously unconnected cells to form a new interface. T1 transitions allow the tissue to relax shear stresses through cell rearrangements rather than cell elongation. A T2 transition, also known as cell extrusion, occurs as a cell shrinks to zero area and is replaced by a single vertex. The mechanical state of the tissue is controlled by both topological rearrangements driven by active processes and geometric frustration. Both types of processes can drive transitions between rigid and fluid states. Here we neglect topological rearrangements to focus on the role of geometry.

We further simplify the model by assuming that all cells have the same preferred area A_α0 = A₀ and all edges have the same tension Λ_ij = Λ. The interfacial energy can then be written in terms of the cell perimeter, where the factor of arises because the interfacial energy of each edge is shared by two cells. The tissue energy can then be recast in the form

where is the preferred perimeter, and E₀ is a constant term obtained from completing the square. Since we care about the gradient of energy and not the absolute value, we discard E₀ in the following.

Finally, we work in dimensionless units by normalizing the energy with KA₀² and lengths with . The dimensionless tissue energy is then given by

where a_α = A_α/A₀, r = Γ/KA₀ is the rigidity ratio, and is the target shape index of the cell.

1.2 Deformation protocol

To numerically obtain the elastic moduli, we simulate the mechanical response of the vertex model under different deformations using a tissue of 4 hexagonal cells in a periodic box of lengths L_x(0) and L_y(0), and area A(0) = L_x(0)L_y(0) determined by energy minimisation, and implemented in the Surface Evolver software.³⁸ For the incompatible regime, the ground state is a regular hexagonal cell. For the compatible regime, while the ground state is degenerate, we use a hexagon with 120° angles between edges and with the edge lengths determined by energy minimisation. First, we use an intermediate rigidity ratio of r = 0.1, and test the response across a range of preferred values of the shape index, from p₀ = 0 to p₀ = 4.6, covering both the compatible and incompatible regimes.

To calculate the shear modulus, we deform the ground state (Fig. 2 top-left) by applying an initially affine deformation to vertices and the boundaries: x_i(ε) = (1 + ε/2)x_i(0), y_i(ε) = (1 + ε/2)⁻¹y_i(0), and L_x(ε) = (1 + ε/2)L_x(0) and L_y(ε) = (1 + ε/2)⁻¹L_y(0), where ε = 0.001 (Fig. 2 top-middle). We then allow the vertex positions to relax to an energy minima (Fig. 2 top-middle), and record the change in tissue energy δE before and after the deformation. The shear modulus is then numerically estimated by .

Fig. 2 Strain protocols for measuring elastic moduli of the vertex model. (top, middle, bottom) From the ground state, the periodic box lengths and vertex positions are transformed and constrained according to an affine transformation, shown by the arrows. From the constrained state, the system is relaxed according to tissue-scale or box constrained. (Top) The shear modulus is calculated by applying a shear transformation to the box. In the constrained state, every edge has the same tension, producing a net force on the vertices, hence this is not a force-balanced state. After relaxation, forces are balanced through a non-affine transformation on the vertices. During relaxation the box size is fixed. (Middle) The bulk modulus is calculated by applying an isotropic expansion to the box and vertices. During relaxation the box size is fixed. (Bottom) The Young's modulus and Poisson's ratio are calculated by applying a uniaxial strain to the box and vertices. During relaxation the height of the box may change and vertices may move.

In the ground state of the incompatible regime, cell edges are under tension and meet at 120° angles. After the initial affine deformation, the angles change and the tissue is no longer in a force-balanced configuration (Fig. 2 top-middle). As we allow the tissue to relax, it responds with a non-affine deformation; vertices which are of the same y-coordinate alternate between moving left and moving right during relaxation, returning the angles between edges to a stable 120° configuration (Fig. 2 top-right). Such a deformation cannot be described by a single affine transformation, but by two affine transformations applied to different subsets of vertices.³³

To demonstrate the importance of this relaxation step, we report the response to two types of deformation protocols: (i) “constrained” deformations which are obtained where after deformation of the bounding box the cell vertices are not allowed to move to minimise the energy of the tissue, and (ii) “relaxed” deformations where the vertices are allowed to adjust their position to achieve force balance and the global tissue shape remains controlled by the geometry of the deformed box. Note that in the compatible regime the relaxed state can also be achieved by allowing the tissue to change its shape,²⁵ and the resulting linear elastic constants are the same.

For an intermediate rigidity ratio r = 0.1, we find that the shear modulus decreases as p₀ increases and becomes zero at the transition to the compatible regime. In particular, the relaxation step allows cells to decrease their perimeter, and thus energy, resulting in a relaxed shear modulus that is softer than in the constrained case (Fig. 3a). In the compatible regime, the tissue is initially under no tension since the preferred perimeter is achieved. Upon straining the tissue, the perimeter increases and tissue energy increases. The larger the initial perimeter, the higher the change, resulting in constrained shear modulus that increases with p₀. When the tissue is able to relax, the vertices move to reduce the perimeter until the preferred perimeter is achieved again, allowing for the net energy to remain constant, leading to a zero shear modulus.

Fig. 3 Non-affine deformations allow for a softer mechanical response. (a) Shear modulus G, (b) bulk modulus K, (c) Young's modulus Y, and (d) Poisson's ratio ν against target shape index p₀ for a rigidity ratio r = 0.1. The constrained values represent elastic moduli where vertices are constrained by the given deformation. The relaxed values lines represent the moduli allowing for non-affine deformations, where vertices may relax, subject to the boundary conditions. Dots represent simulated values. Lines represent analytic values. Shaded regions show the range of possible values in the constrained case, depending on the initial shape of the cells.

To calculate the bulk modulus, we apply the isotropic transformation , and and where ε = 0.001, such that A(ε) = (1 + ε)A(0). During the relaxation step, we allow the vertices to move, with the box lengths fixed. The bulk modulus is then given by .

In the incompatible regime, force balance requires a constant 120° angle between edges, thus the tissue expands isotropically. We find that the bulk modulus increases as the target shape index p₀ increases, and is equal between the relaxed and constrained cases (Fig. 3b).

In the compatible regime, the deformation initially increases the perimeter. During the relaxation step, the tissue responds in a non-affine way to restore its perimeter to its preferred value and so energy change only arises from the area term and we have a bulk modulus K = 1. Interestingly, this is lower than the bulk modulus in the incompatible regime just before the transition and thus, there is a discontinuity in the bulk modulus as p₀ changes. In contrast, the constrained case is unable to relax the cells perimeters and so has a higher bulk modulus and does not exhibit the discontinuity (Fig. 3b).

Next, we apply a uniaxial deformation to calculate the Young's modulus and Poisson's ratio: x_i(ε) = (1 + ε)x_i(0) and L_x(ε) = (1 + ε)L_x(0) (Fig. 1c). We then allow the vertex positions and box height L_y(ε) to relax to minimise energy. The Young's modulus is given by and the Poisson's ratio by . Note that this definition of the Poison's ratio is equivalent to that in 2D elasticity and therefore its values are limited between −1 < ν < 1. The extreme case ν = 1 corresponds to incompressible solid, analogous to the case of ν_3d = 0.5 for incompressible 3D solids.

Again, the tissue undergoes a similar non-affine relaxation as under shear strain, reducing the shear modulus compared to the constrained case (Fig. 3c). In this case, though, we find that the Young's modulus is non-monotonic. For p₀ close to zero, the Young's modulus increases as p₀ increases. For higher p₀, increasing p₀ further decreases the Young's modulus towards zero at the transition point, after which the Young's modulus is zero. However, in the constrained case the Young's modulus increases after the transition point due to the increased bulk and shear moduli. Interestingly, the Poisson's ratio begins negative for small p₀ and increases towards a value of 1 as p₀ increases, before remaining 1 in the compatible regime (Fig. 3c). In the constrained case, the Poisson's ratio is actually lower than in the relaxed case for small p₀, in particular, in the compatible case the Poisson's ratio decreases as p₀ increases while for a relaxed tissue it remains 1.

This phenomena highlights the counter-intuitive nature of VM mechanics. In classical elasticity ν = 1 corresponds to incompressible solids, commonly considered as very stiff. Here we find that the tissue approaches ν = 1 for higher values of p₀ corresponding to compatible tissue with floppy response. This seeming contradiction is resolved by noting that in this limit cells can accommodate rest area and perimeter simultaneously and therefore upon deformation their area remains intact, just as in incompressible solids.

The simulations highlight the complex mechanical behaviour of the vertex model to applied tissue-level strains, both in its elastic moduli and the vertex-level non-affine deformations while relaxing the energy. The non-affine relaxation step allows the tissue to reduce its elastic moduli in the compatible regime. In particular, while the shear and Young's moduli are zero at the transition point, increasing p₀ increases the moduli in the constrained case, while they remain zero in the relaxed case. However, the simulations do not give an intuitive understanding for why the bulk modulus is discontinuous, or why we can get a negative Poisson's ratio. Thus, in the remainder of the paper, we develop a mean-field theory of the vertex model that can account for non-affine relaxation of the tissue under strain to derive analytic expressions for the elastic moduli and understand the source of the complex phenomena mentioned above.

2 Vertex model: mean-field theory and ground states

2.1 Mean-field theory of vertex model

To understand the numerical results, we construct a mean-field theory by assuming that all cells responds equally. In this case the tissue energy is just E_tissue = NE and one can simply consider the energy E of a single cell, given by

Each cell consists of horizontal edges of length l₁ and diagonal edges of length l₂, with ϕ the angle between horizontal and diagonal edges (Fig. 4a). This parameterization captures the behavior of the tissue observed in our numerical simulations, where it is the angle between edges that changes during relaxation. Although cells have additional degrees of freedom, the description in terms of these three degrees of freedom is sufficient to capture the ground states of the tissue VM, and the response of the tissue under shear and bulk deformations in simulations (Fig. 2). To both examine the ground states and the response to deformation, it is convenient to parametrize each cell in terms of the height h and width w, as shown in Fig. 4a, given by

h = 2l₂ sin ϕ, w = l₁ − l₂ cos ϕ. (5)

Each cell then contributes an area a = wh to the tissue. We stress that the angle ϕ is distinct from the shear tilt angle θ introduced previously in ref. 25, where cells may tilt or untilt in order to change their perimeter. While the shapes obtained from an initial regular hexagon by varying the shear tilt angle θ correspond to affine deformations of the hexagon, those parametrized by ϕ generally correspond to non-affine deformations of the regular hexagon. Inverting eqn (5), we obtainCell area and perimeter can then be written as

a = hw, (8)

p = 2w + hf(ϕ), (9)

whereresulting in an energyThis form makes it evident that the VM energy is underconstrained as area and perimeter do not uniquely determine cell shape.

Fig. 4 Shape parameterization of the vertex model and ground states. (a) Schematic of the vertex model and cell shape parametrization. Cells are defined by the lattice height h, width w, and angle between edges ϕ. (b) The ground state in the solid state is a regular hexagonal lattice, with ϕ = 2π/3. (c) The ground state used in the soft state with ϕ > 2π/3. (d–f) Cell area (d), cell perimeter (e), and edge tension (f) vs. target shape index p₀ for various values of the rigidity ratio r. Dots represent simulated values, lines are the analytical results.

2.2 Ground states

The ground state configurations are obtained by minimizing the energy with respect to the cell width w, height h, and angle ϕ and are solutions of the three coupled equationsAs shown in previous work, we find a transition at p₀ = p_* between two distinct states. For a regular lattice of n-sided polygons p_* is given by the isoperimetric value with . The isoperimetric inequality p ≥ p_*(n) provides a lower bound on the perimeter of a regular n-sided polygon for given area.³⁹ For p₀ > p_*(n) the cell is in a geometrically compatible regime, where both preferred area and perimeter may be achieved, and the tissue has zero shear modulus^22,31,33,40 (Fig. 1b). For p₀ < p_*(n) the cell is in an incompatible regime, where both preferred area and perimeter cannot be simultaneously satisfied, and the tissue behaves like a solid by resisting shear deformation.^21,25,31,32 The corresponding ground state of the tissue is a lattice of identical hexagonal cells (Fig. 1b). As p₀ is further lowered the cell may become unstable and collapse to zero area and perimeter (Fig. 1b). Additionally, in a small range of parameters near the collapsing region more exotic ground states exist, with mixed lattices of square and octagonal, or dodecahedral and triangular cells providing lower energy than hexagonal cells.²²

2.2.1 Compatible state, p₀ > p_*(6). For p₀ > p_*(6) eqn (12) are identically solved by hw = 1 and p = 2w + hf(ϕ) = p₀, and the zero ground state energy vanishes (Fig. 4c–e). We refer to this situation as the compatible state. The ground state configuration is a family of 6-sided polygons parametrized by the angle ϕ, withwhere both roots are acceptable solutions for a given value of ϕ, corresponding to either tall and thin or short and wide cells. It is evident from eqn (15) and (16) that such a solution exists provided p₀² ≥ 8f(ϕ). The function f(ϕ) has a minimum at with corresponding to . At this value of p₀ there is a single zero energy solution that corresponds to a hexagon of unit area. For p₀ > p_* there is degenerate continuum of zero energy solutions corresponding to deformed hexagons of unit area, perimeter p₀ and with ϕ_m determined by p₀² = 8f(ϕ_m) (Fig. 4c). There exist many other parameterizations that can give ground state shapes in the compatible regime, for example, cells becoming tall and thin, cells decreasing the angle ϕ to increase their perimeter, or cells tilting as in ref. 25.

2.2.2 Incompatible state, p₀ < p_*(6). For p₀ < p_* the cell cannot simultaneously realize the target area and perimeter. We refer to this situation as the incompatible state. Eqn (12) requires f′(ϕ) = 0, with solution such that . An intuitive explanation for this fixed angle is that all edges are under identical tension, and so by force balance a junction of three edges must have equally spaced angles. Eqn (13) and (14) then imply that . This gives a perimeter p = 4w and area which means that the cells are regular hexagons in the incompatible state (Fig. 4b).

We can combine eqn (13) and (14) to obtain a cubic equation for the perimeter

with and a = p²/p_*². The cubic equation can be solved perturbatively in the limit of low and high rigidity ratio r.

At low rigidity ratio, i.e., r ≪ 1/p_*², we find

The cell remains close in shape to a hexagon of unit area, with a reduction of the perimeter relative to the value p_* (Fig. 4d and e). The tension, given by decreases monotonically as p₀ increases and vanishes at p₀ = p_*, where the cell reaches the compatible state and is under no tension (Fig. 4f). If p₀ becomes too small, the stable configuration collapses to a point with zero area and perimeter (Fig. 1b).

For high rigidity ratio, i.e., r ≫ 1/p_*², the perimeter and area can be expanded in inverse powers of r, with the result

In this limit the cell shape index is close to the target shape index (Fig. 4d and e). The tension is nonmonotonic in p₀, increasing from almost no tension at p₀ = 0 before decreasing as p₀ approaches p_* (Fig. 4f). As the rigidity ratio increases the tension saturates to the value which has a maximum value of at (Fig. 4d). For p₀ ≤ 0, the cell collapses to a point with zero area and perimeter (Fig. 4b).

Finally, we note that when topological transitions are allowed, tissues may also unjam and undergo a solid-to-liquid phase transition when cell rearrangements cost zero energy. In disordered realizations of the VM, the unjamming transition occurs at p₀ ≈ 3.81, a value close to, but slightly larger than p_*(5).²⁶ Intuitively, for the tissue to rearrange in a T1 transition with zero energy barrier, two hexagonal cells must momentarily lose an edge and become pentagons while still maintaining their preferred perimeter and area. Cell motility can further promote fluidity and lower the transition point.²³

3 Mechanical response of the vertex model

3.1 Deformations protocol

It is evident from eqn (11) that area and perimeter (or equivalently height and width of the box shown in Fig. 4a) do not uniquely specify a polygonal shape. In the compatible regime there is a family of zero energy shapes corresponding to either tilted polygonal shapes obtained by affine deformations or non-affinely deformed polygons parametrized by the angle ϕ. In the incompatible regime, if only affine deformations are allowed, both the ground state and each deformed state are unique for fixed area and perimeter. Allowing non-affine deformations introduces, however, additional degrees of freedom that can lower the energy for a given set of parameters.

In the following we examine the response of a tissue initially in a ground state to an externally imposed strain. The deformation is imposed globally on the tissue by changing the shape of the bounding box. Such a deformation uniformly changes the shape of the cells and generally results in a state where individual vertices are no longer force balanced (Fig. 2 top-middle). We will refer to this state as the “constrained” deformed state. Due to the presence of hidden degrees of freedom the system can, however, lower its energy and relax a state of local force balance. In the compatible regime this relaxation can occur via motion of the vertices that correspond to non-affine deformations (Fig. 2 top-right). In the compatible regime the relaxation can occur either via non-affine deformations with fixed box shape or through a global tilting of the tissue, which entails affine cell deformations, as in Hernandez et al.²⁵ The elastic constants measured in the “relaxed” state of the compatible regime are the same for the two relaxation protocols.

Operationally, constrained deformations are achieved by first fixing either the cell height, width, or both, and then transforming the vertices according to the given deformation, as done in Staple et al.²² We prevent spontaneous tilting of the tissue, which can be used to soften the mechanical response using only affine deformations in the compatible regime.²⁵

We next evaluate the various elastic moduli of the vertex model. As we will see below, a new result of our work is that in the incompatible regime cells can find new deformed states by relaxing through non-affine deformations (Fig. 2), resulting in a softer response than obtained in previous studies²² (Fig. 3).

3.2 Shear modulus

To calculate the shear modulus of the tissue, we apply an area-preserving pure shear deformation, corresponding to w → w(1 + ε/2) and h → h(1 + ε/2)⁻¹, with ε the strain. We allow for a non-affine deformation to relax the tissue by minimising energy with respect to the angle ϕ for each value of strain (Fig. 4). The shear modulus is defined asAs area is preserved under pure shear, we only need to consider the energy cost due to changes in perimeter.

3.2.1 Compatible State, p₀ > p_*. In the compatible case, cells can accommodate shear and maintain their area and perimeter at the target values a = 1 and p = p₀ by changing shape, i.e., by adjusting the angle ϕ to a value other than 2π/3. The perimeter of the deformed cell is given by p(ε,ϕ) = 2w(1 + ε/2) + hf(ϕ)/(1 + ε/2). The cell can maintain p = p₀ by deforming to a new compatible ground state corresponding to an angle ϕ* given by the solution of p(ε,ϕ*) = p₀. Clearly the energy remains zero, demonstrating that the shear deformation cost no energy and

G = 0, (23)

for all rigidity ratios and p₀ > p_* (Fig. 5a). This of course only holds up to a maximum value of strain determined by the angle ϕ_m(p₀). Beyond this value the fluid-like compatible tissue stiffens and acquires a finite shear modulus, as mentioned by ref. 41, and recently by ref. 40.

Fig. 5 Elastic moduli of the vertex model. (a) Shear modulus, (b) bulk modulus, (c) Young's modulus, and (d) Poisson's ratio against target shape index p₀ and rigidity ratio r.

However, if the shape of the cell is constrained then applying a shear transformation could increase the perimeter. To calculate the shear modulus under an affine transformation, we must consider the changes to the cell perimeter from its initial configuration. The shear modulus is defined by

and since for ε = 0 we have a = 1 and p = p₀ this simplifies toSince an affine transformation also changes the angle between the edges, we must consider the effect of the transformation on each edge when calculating the change in perimeter. The length of the two edges shown in Fig. 4a change asandwith first derivativesandThusand

The cell angle ϕ defines a family of solutions with as does the choice between the ± branch in our solution for cell height and, and so we obtain a range of values of the constrained shear modulus, and thus Young's modulus and Poisson's ratio, as ϕ is varied. We find that as p₀ is increased, the minimum shear modulus in the constrained case remains 0 while the maximum possible shear modulus increases (Fig. 3a).

3.2.2 Incompatible case, p₀ < p_*. In the incompatible case cell edges are under uniform tension. By force balance, this implies that the angle between edges remains even under small deformations at the tissue scale. The ground state configuration is a regular hexagon with perimeter . Using eqn (25) and the relations for height and width in terms of the perimeter, and we obtainIn the rigid, incompatible state the shear modulus is a monotonically increasing function of r and vanishes at the transition p₀ = p_* (Fig. 3a and 5a), in agreement with earlier results.³³ The non-affine deformations of the relaxed tissue allow for a softer response of the tissue, with the shear modulus being a factor of 3/2 stiffer when only considering vertices constrained by the affine shear strain,^22,41 as confirmed by simulations (Fig. 3a).

For small rigidity ratio (r ≪ 1/p_*²), we can expand G in powers of r, with the result

The opposite limit of large rigidity ratio (r ≫ 1/p_*²) yields

3.3 Bulk modulus

To calculate the bulk modulus, we change the area a → a(1 + ε) by rescaling the height and width , and allow the angle ϕ to vary to minimize the deformation energy. The bulk modulus is then given bywithTo evaluate this expression we need to consider separately the compatible and incompatible states.

3.3.1 Compatible state, p₀ > p_*. We have previously shown that in the compatible case the angle ϕ can adjust to maintain a fixed cell perimeter under small deformations. Thus and and the bulk modulus is simply

K = 1, (37)

for all r and p₀ > p_* (Fig. 3b and 5b).

By contrast, if we allow for only affine deformations then the perimeter expands isotropically . In this case and giving a bulk modulus equal to

which can be significantly higher than the non-affine result for high r, emphasising the need to consider non-affine displacements (Fig. 2).

3.3.2 Incompatible state, p₀ < p_*. In the incompatible state for p₀ < p_* the angle that minimizes energy remains for small perturbations to cell height and width to ensure tension balance at the cell vertices. The cell then expands isotropically, such that resulting in a bulk modulusshown in Fig. 3b and 5b. We note that as p₀ approaches the critical value p_*(6) from below, the bulk modulus has the value . On the other hand, in the compatible regime K = 1. Thus the bulk modulus exhibits a jump discontinuity at the critical point separating compatible and incompatible states. In contrast, if vertex positions are fixed by uniform dilation without relaxation, the bulk modulus is continuous and higher in the compatible region (Fig. 3b).

In the limit of low rigidity ratio (r ≪ 1/p_*²), we find

The bulk modulus increases with p₀ up to the critical value, at which point it discontinuously jumps to 1 for all p₀ > p_*, independent of r. Interestingly, for the bulk modulus of the incompatible solid is lower than that of the compatible fluid, suggesting that contractility can actually reduce the bulk stiffness of the tissue. Additionally, increasing the rigidity ratio further reduces the bulk modulus for low p₀.

In the limit of high rigidity ratio (r ≫ 1/p_*²), we find

Thus the bulk modulus increases with the rigidity ratio.

We find that for low p₀ the bulk modulus can be significantly lower in the rigid than in the fluid state, which can affect the rate of spreading monolayers.⁴² This is due to the fact that at small p₀ cells have a smaller area while the energy required to deform the cell is proportional to the square of the area change. Therefore it costs more energy to strain a single cell than to strain two cells with half the area, similar to the reduction in effective stiffness obtained when springs are placed in series.

The bulk modulus is also a non-monotonic function of the rigidity ratio. For small r the bulk modulus decreases with r as the size of the cell decreases, reaching a minimum near r = 2/p_*² ≈ 0.144 before increasing linearly in r for high r, due to the growing contribution from the perimeter elasticity.

3.4 Young's modulus and poisson's ratio

Next we calculate the Young's modulus and Poisson's ratio (Fig. 2 bottom row) by stretching the width of the cell w → w(1 + ε) while allowing the cell height h and angle ϕ free to minimize the energy. The Young's modulus is defined asand the Poisson's ratio asTo evaluate Y and ν we use the relationship between the linear elastic constants, and which have been shown to hold away from the critical point.

3.4.1 Compatible state, p₀ > p_*. In the compatible state, the ground state degeneracy allows the cell to achieve the target shape index and area for small strain by reducing cell height and finding values of the angle ϕ different from 2π/3, i.e., by changing its shape, with no energetic cost. As a result for p₀ > p_* we find

Y = 0, ν = 1, (44)

for all r (Fig. 3c, d and 5c, d). When constrained to affine only deformations, the Young's modulus can take a range of values from zero to a maximum value which increases with p₀, due to the increasing shear modulus. Similarly, the Poisson's can range from a maximum of 1 to a minimum value which decreases with p₀.

3.4.2 Incompatible case, p₀ < p_*. It has been shown that away from the critical point the elastic constants of the VM satisfy the familiar relation of linear elasticity of isotropic solids.²⁵ We can therefore use the relations and to evaluate Y and ν in the incompatible regime, with the resultWe find that the Young's modulus is a non-monotonic function of both target shape index and rigidity ratio (Fig. 3c and 5c). The nonmonotonicity with p₀ is most pronounced at intermediate values of r, where at small p₀ the Young's modulus increases, rather than decrease, with increasing p₀.

At both high and low rigidity ratio, the Poisson's ratio remains close to 1 for all p₀, indicating that the cell preserves its area under deformations (Fig. 3d and 5d). At intermediate values of the rigidity ratio and small p₀, the nonmonotonicity of the Young's modulus results in a negative Poisson's ratio, which indicates that a tissue stretched in the x-direction, also expands in the y-direction.

In comparison to the constrained response, we find that the relaxed response gives a softer Young's modulus for all rigidity ratio and target shape index values, since the shear modulus is also softer by a factor of (Fig. 3c). Similarly, the Poisson's ratio is higher in the compatible state for all p₀ < p_* (Fig. 3d).

In the limit of low rigidity ratio (r ≪ 1/p_*²), the Young's modulus and Poisson's ratio are given by

showing that when p₀ is increased towards the critical point from the solid side (p₀ → p_*⁻) the Young's modulus vanishes and the Poisson's ratio increases towards 1.

In the limit of high rigidity ratio (r ≪ 1/p_*²), we obtain approximate expression by expanding in 1/r as

The Young's modulus also has a maximum value of 2 at p₀ = 0.

3.5 Origin of negative Poisson's ratio

We can understand why certain cell parameters give a positive or negative Poisson's ratio by looking at how the energy gradient with respect to cell height changes as we change the cell width. The gradient can be thought of as the effective force acting on the height of the cell, given byIn the ground state, this will be zero. Then, if we vary the cell width but keep the height fixed we can measure the change in force asWhen this value is positive, then as cell width is increased, the effective force acting on cell height increases and so the cell height will decrease as it relaxes to the energy minimum. Consequently, the sign of this value is the same sign as the Poisson's ratio.

The second term comes from the perimeter contribution in the VM and accounts for energy changes due to perimeter elasticity. Since p ≥ p₀, the cell is under tension and the perimeter term aims to shrink the cell. Increasing the width further increases the perimeter and so tension increases, providing more force to shrink the cell. Thus, the perimeter elasticity always acts to shrink the cell and contributes to a positive Poisson's ratio.

The first term, 2hw − 1 = 2a − 1, represents the energy change due to area elasticity. We can write this as 2w(h − 1/2w), which we can think of as an spring like force with stiffness 2w and target height 1/2w. As width increases, the height becomes closer to the target height, reducing the strain. At the same time, the effective stiffness 2w increases, increasing the pressure. Thus there is a trade off between less restoring force on the cell area versus increased effectiveness of changes in cell height. The net effect on whether this increases or decreases the perimeter depends on the size of the cell area: when cell area the area term acts to increase cell height when width is increased, and for the area term reduces the cell height.

For high rigidity ratio, the perimeter term dominates and so an increase in cell width leads to a reduction in cell height. For low rigidity ratio, the area term dominates and the cell area is close to 1 (Fig. 4d), thus an increase in cell width reduces the area pressure and cell height decreases. However, in the intermediate regime for low p₀ cell area is small, meaning the area term acts to increase cell width, and the area contribution and perimeter contributions are of comparable size, resulting in negative Poisson's ratio.

We can calculate the transition to a negative Poisson's ratio exactly at p₀ = 0, which corresponds to the situation where the contribution to cell edge tension from cortical contractility and cell–cell adhesion precisely balance. In this limit the equation for the ground state perimeter, eqn (17), becomes

with solutionfor r < 2/p_*² ≈ 0.144. For r > 2/p_*² the cell is unstable and collapses to zero area. We can calculate whether cell height increases or decreases when width is increased by calculating how the effective force on the height, changes with width. Substituting our formula for area we findwhere we have used . Thus for and p₀ = 0 the tissue has a negative Poisson's ratio.

4 Conclusions

In this paper we studied the linear response of the 2D vertex model by calculating the shear modulus, bulk modulus, Young's modulus and Poisson's ratio, using a mean-field approach that allows for a class of non-affine deformations, which agree well with numerical simulations. We also provide approximate expressions in the limit of high and low rigidity ratio. Our calculations match previous results showing a rigidity transition controlled by purely geometric effects and tuned by the target shape index p₀.

For cells in the incompatible case, p₀ < p_* ≈ 3.772, the tissue has a finite shear modulus which decreases with p₀. For cells in the compatible case, p₀ > p_*, the shear modulus becomes zero for all p₀. However, when the tissue is constrained by the deformation, we find a stiffer mechanical response to shear in the incompatible case, and the compatible case can have a finite shear modulus, which depends on the initial configuration of the cells, and that increases with p₀.

The bulk modulus of the tissue increases with p₀ in the incompatible regime, and then has a jump discontinuity at p₀ = p_*, where it changes from a larger value in the solid state to a value of 1 in the fluid state. In the incompatible case, cells perimeters increase upon isotropic expansion of the tissue, but in the compatible regime cells can change shape to preserve their perimeter under small changes in area. However, this discontinuity is not observed when the tissue is constrained. We also find that in the incompatible regime, the bulk modulus can decrease below 1 for small rigidity ratio. This indicates that cell contractility can reduce the stiffness of the tissue, resulting in a larger bulk modulus in the “soft” phase than in the “solid” phase.

Under uniaxial strain, we find the Young's modulus of the tissue can be non-monotonic with respect to p₀, initially increasing and then decreasing towards zero in the incompatible case. The Poisson's ratio can become negative for small p₀ and intermediate rigidity ratio, as cells can reduce their energy more by increasing their area through an orthogonal expansion, than by reducing their perimeter. Within the compatible regime, the tissue has zero Young's modulus and Poisson's ratio equal to one. However, when only constrained deformations are allowed the tissue can have a finite Young's modulus in the compatible regime, similar to the shear modulus.

Our results highlight the complex linear elastic behaviour that can arise from the simplest version of the vertex model due to its underconstrained nature. For simplicity, we have assumed that cells are regularly arranged and that we only have small strains. This analysis might be most applicable in tissues with a regular crystalline structure, such as the Drosophila pupal wing? However, it would be interesting to extend our calculations to the case of disordered cell networks. Additionally, we highlight the importance of allowing for unconstrained degrees of freedom, in this case non-affine deformations, to relax the system and give a softer mechanical response to strain. The constrained case may be thought of as the short-time response of the tissue to strain, and the relaxed case as the long-time limit.

Finally, we note that the emergence of rigidity observed in this work draws a direct link with the rigidity of mechanical frames and granular models, where mechanical stability occurs at critical coordination number, or at finite strains, and are normally accompanied by a discontinuous jump in the bulk modulus. Nevertheless, the results of our work shows that geometric incompatibility is a crucial ingredient that has to be taken into account for an estimation of the critical coordination number, and is left for a future work.^43,44

In conclusion, our work demonstrates that the vertex model, thought of as a collection of geometric constraints rather than a reference ground state structure, can engender interesting linear mechanical responses. The linear response exhibits a strong non-affine contribution under uniaxial compression and shear, as well as a negative Poisson's ratio. Typically, these two phenomena in crystalline solids require special lattice constructions, whereas in the vertex model exotic mechanical response can be achieved by tuning the relative competition between area and perimeter constraints via r and geometric compatibility via p₀.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

M. C. M. and A. H. were supported by the National Science Foundation Grant No. DMR-2041459. M. M. was supported by the Israel Science Foundation grant No. 1441/19. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. Open Access funding provided by the Max Planck Society.

References

1. A. C. Martin, M. Kaschube and E. F. Wieschaus, Nature, 2009, 457, 495–499.
2. R. Etournay, M. Popović, M. Merkel, A. Nandi, C. Blasse, B. Aigouy, H. Brandl, G. Myers, G. Salbreux and F. Jülicher, et al. , eLife, 2015, 4, e07090.
3. S. J. Streichan, M. F. Lefebvre, N. Noll, E. F. Wieschaus and B. I. Shraiman, eLife, 2018, 7, e27454.
4. E. Maniou, M. F. Staddon, A. R. Marshall, N. D. Greene, A. J. Copp, S. Banerjee and G. L. Galea, Proc. Natl. Acad. Sci. U. S. A., 2021, 118, e2023163118.
5. M. Poujade, E. Grasland-Mongrain, A. Hertzog, J. Jouanneau, P. Chavrier, B. Ladoux, A. Buguin and P. Silberzan, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 15988–15993.
6. A. Brugués, E. Anon, V. Conte, J. H. Veldhuis, M. Gupta, J. Colombelli, J. J. Muñoz, G. W. Brodland, B. Ladoux and X. Trepat, Nat. Phys., 2014, 10, 683–690.
7. R. J. Tetley, M. F. Staddon, D. Heller, A. Hoppe, S. Banerjee and Y. Mao, Nat. Phys., 2019, 15, 1195–1203.
8. V. Ajeti, A. P. Tabatabai, A. J. Fleszar, M. F. Staddon, D. S. Seara, C. Suarez, M. S. Yousafzai, D. Bi, D. R. Kovar and S. Banerjee, et al. , Nat. Phys., 2019, 15, 696–705.
9. P. Friedl and D. Gilmour, Nat. Rev. Mol. Cell Biol., 2009, 10, 445–457.
10. E. N. Arwert, E. Hoste and F. M. Watt, Nat. Rev. Cancer, 2012, 12, 170–180.
11. G. Salbreux, G. Charras and E. Paluch, Trends Cell Biol., 2012, 22, 536–545.
12. M. Murrell, P. W. Oakes, M. Lenz and M. L. Gardel, Nat. Rev. Mol. Cell Biol., 2015, 16, 486–498.
13. B. Ladoux and R.-M. Mège, Nat. Rev. Mol. Cell Biol., 2017, 18, 743–757.
14. F. Graner and J. A. Glazier, Phys. Rev. Lett., 1992, 69, 2013.
15. M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143.
16. S. Banerjee and M. C. Marchetti, Cell Migrations: Causes Functions, 2019, 45–66.
17. R. Alert and X. Trepat, Ann. Rev. Condens. Matter Phys., 2020, 11, 77–101.
18. S. Alt, P. Ganguly and G. Salbreux, Philos. Trans. R. Soc., B, 2017, 372, 20150520.
19. T. Nagai and H. Honda, Philos. Mag. B, 2001, 81, 699–719.
20. A. G. Fletcher, M. Osterfield, R. E. Baker and S. Y. Shvartsman, Biophys. J., 2014, 106, 2291–2304.
21. R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton and F. Jülicher, Curr. Biol., 2007, 17, 2095–2104.
22. D. B. Staple, R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton and F. Jülicher, Eur. Phys. J. E: Soft Matter Biol. Phys., 2010, 33, 117–127.
23. D. Bi, X. Yang, M. C. Marchetti and M. L. Manning, Phys. Rev. X, 2016, 6, 021011.
24. L. Hufnagel, A. A. Teleman, H. Rouault, S. M. Cohen and B. I. Shraiman, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 3835–3840.
25. A. Hernandez, M. F. Staddon, M. J. Bowick, M. C. Marchetti and M. Moshe, Phys. Rev. E, 2022, 105, 064611.
26. D. Bi, J. Lopez, J. M. Schwarz and M. L. Manning, Nat. Phys., 2015, 11, 1074–1079.
27. D. M. Sussman, M. Paoluzzi, M. C. Marchetti and M. L. Manning, EPL, 2018, 121, 36001.
28. S. Tong, N. K. Singh, R. Sknepnek and A. Kosmrlj, PLoS Comput. Biol., 2022, 18(5), e1010135.
29. C. Duclut, J. Paijmans, M. M. Inamdar, C. D. Modes and F. Jülicher, Cells Dev., 2021, 203746.
30. C. Duclut, J. Paijmans, M. M. Inamdar, C. D. Modes and F. Jülicher, Eur. Phys. J. E: Soft Matter Biol. Phys., 2022, 45, 1–10.
31. M. Moshe, M. J. Bowick and M. C. Marchetti, Phys. Rev. Lett., 2018, 120, 268105.
32. M. Merkel, K. Baumgarten, B. P. Tighe and M. L. Manning, Proc. Natl. Acad. Sci. U. S. A., 2019, 116, 6560–6568.
33. N. Murisic, V. Hakim, I. G. Kevrekidis, S. Y. Shvartsman and B. Audoly, Biophys. J., 2015, 109, 154–163.
34. M. L. Falk and J. S. Langer, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1998, 57, 7192.
35. B. Utter and R. Behringer, Phys. Rev. Lett., 2008, 100, 208302.
36. W. G. Ellenbroek, Z. Zeravcic, W. van Saarloos and M. van Hecke, EPL, 2009, 87, 34004.
37. A. Zaccone and E. Scossa-Romano, Phys. Rev. B, 2011, 83, 184205.
38. K. A. Brakke, Exper. Math., 1992, 1, 141–165.
39. R. Osserman, Bull. Am. Math. Soc., 1978, 84, 1182–1238.
40. J. Huang, J. O. Cochran, S. M. Fielding, M. C. Marchetti and D. Bi, Phys. Rev. Lett., 2022, 128, 178001.
41. R. Farhadifar, Dynamics of cell packing and polar order in developing epithelia, PhD thesis, Technische Universität Dresden, 2009.
42. M. F. Staddon, M. P. Murrell and S. Banerjee, Soft Matter, 2022, 18, 7877–7886.
43. O. K. Damavandi, V. F. Hagh, C. D. Santangelo and M. L. Manning, Phys. Rev. E, 2022, 105, 025003.
44. O. K. Damavandi, V. F. Hagh, C. D. Santangelo and M. L. Manning, Phys. Rev. E, 2022, 105, 025004.

---