Mechanical characterization of disordered and anisotropic cellular monolayers

We consider a cellular monolayer, described using a vertex-based model, for which cells form a spatially disordered array of convex polygons that tile the plane. Equilibrium cell configurations are assumed to minimize a global energy defined in terms of cell areas and perimeters; energy is dissipated via dynamic area and length changes, as well as cell neighbor exchanges. The model captures our observations of an epithelium from a Xenopus embryo showing that uniaxial stretching induces spatial ordering, with cells under net tension (compression) tending to align with (against) the direction of stretch, but with the stress remaining heterogeneous at the single-cell level. We use the vertex model to derive the linearized relation between tissue-level stress, strain, and strain rate about a deformed base state, which can be used to characterize the tissue’s anisotropic mechanical properties; expressions for viscoelastic tissue moduli are given as direct sums over cells. When the base state is isotropic, the model predicts that tissue properties can be tuned to a regime with high elastic shear resistance but low resistance to area changes, or vice versa.


I. INTRODUCTION
Epithelial tissues have significant roles in embryonic development, tissue homeostasis, and disease development [1]. Recent work has revealed that many critical functions in biological tissues are dependent on the accurate organization of the shapes and packing geometry of the constituent cells [2]. Disturbances in this organization have been associated with problems during embryonic development and diseases in adult life [3,4]. Furthermore, there is evidence that mechanical forces may directly trigger biochemical responses that regulate morphogenetic processes [5,6]. However, due to difficulties in quantifying stresses in tissues, the mechanisms by which tissue behavior emerges from these multiscale feedback processes remain poorly understood.
Continuum descriptions can provide useful insights to tissue-level behavior. For example, elastic-viscoplastic continuum models can capture the solid-and liquidlike response of tissues to small and large deformations over differing time scales [7,8]. However, they are in many cases not built from an explicit physical description of cells. Furthermore, the interactions of multiple cells can lead to rich emergent behavior at the * oliver.jensen@manchester.ac.uk tissue scale, such as yielding and remodeling, that is not easily accessed through conventional continuum frameworks [9,10].
Discrete vertex-based models of epithelia have been a useful tool in linking mechanics to tissue morphology [11][12][13][14][15]. These models have more recently been developed to characterize the mechanical properties of tissues [16,17] and to infer local and global stresses [18][19][20]. This work has predicted interesting long-range mesoscopic mechanical patterning arising purely from the mechanical properties and short-range mechanical interactions of cells within the tissue, which are not seen in traditional continuum descriptions [20,21]. Relationships between discrete models and traditional continuum approaches have been found for spatially periodic cell networks [16,22], while equivalent relationships for disordered tissues have only been partially established for isotropic disordered cell networks [20] or for analogous physical systems such as two-dimensional (2D) dry foams [23].
Many mechanical models of biological tissues assume that the material is isotropic. However, recent observations in the Drosophila melanogaster embryo [24] have provided evidence that biological tissues may exhibit an orientational, as well as positional, structure. Likewise, models of 3D foams have explored how an orientational structure can be introduced through uniaxial deformations [25]. The deformation induces a net stress in the material, leading to a series of irreversible deformations (such as neighbor exchanges). Experimental observations of epithelial tissues have revealed similar patterns of orientational order following stretching [19,26,27].
To explore how deformation induces anisotropy, in this paper we use a variant of a well-studied vertex-based model to quantify the mechanical behavior of a disordered cellular monolayer under an external load. We ignore cell division or motility but take account of dissipation arising from changes in cell area and perimeter, motivated by observations of damping in the cytosol of cells in the Drosophila embryo on a time scale of minutes [28], evidence of shorter but distinct stressrelaxation time scales in the cytosol and cortex [29], subcellular observations of dissipation at cell contacts [30], and evidence of viscoelastic stress relaxation in a freely suspended cultured monolayer [27]. We provide observations of stretched Xenopus embryonic epithelium demonstrating that a uniaxial stretch enforces order in the tissue, in which cells under net tension tend to align their principal axis of stress with the stretch direction and those under compression align perpendicularly. This behavior is captured in simulations, which are further used to quantify the tissue's anisotropy using the deviatoric (shear) component of the global stress. We then derive a linearized stress/strain/strain-rate relationship characterizing the perturbation stress of a prestressed tissue subjected to a small homogeneous deformation. This allows viscoelastic moduli to be computed for an anisotropic disordered cellular monolayer. Finally, we show that the mechanical parameters of an isotropic tissue can be tuned to elicit high shear resistance but low resistance to area changes, or vice versa.

II. METHODS
We use a modification of a popular vertex-based model to describe a planar epithelium [11,13,31], using the notational framework presented in Ref. [20]. Details of our experimental protocol follow a summary of the model.

A. The vertex-based model
The epithelial monolayer M is represented as a spatially disordered planar network of N v vertices, labeled j = 1, . . . ,N v , connected by straight edges and bounding N c convex polygonal cells, labeled α = 1, . . . ,N c . The cells are assumed to have identical mechanical properties described in terms of a preferred area A * 0 , a preferred perimeter L * 0 , a bulk stiffness K * , a cortical stiffness * , a bulk viscosity γ * , and a cortical viscosity μ * . Scaling all distances on A * 0 , the vector from the coordinate origin to vertex j is given by R j . Each vertex is shared by three cells and edges are shared by two cells (excluding cells at the boundary of M). Quantities specific to cell α are labeled by a Greek subscript and defined relative to its centroid R α [ Fig. 1(a)]. Cell α has Z α vertices labeled anticlockwise by i = 0,1,2, . . . ,Z α − 1 relative to R α . We define R i α as the vector from the cell centroid to vertex i, such that Z α −1 i=0 R i α = 0. Anticlockwise tangents are defined by t i α = R i+1 α − R i α (taking i + 1 modulo Z α ), unit vectors along a cell edge byt i α , and outward normals to edges by n i α = t i α ×ẑ (whereẑ is a unit vector pointing out of the plane). The length l i α of an edge belonging to cell α between vertices i and i + 1, the cell perimeter L α , and the cell area A α are given by We note for later reference that [20], Q α is a symmetric tensor characterizing the shape of cell α, satisfying Tr(Q α ) = 1. The dimensionless mechanical energy of an individual cell U α (scaled on K * A * 2 0 ) is assumed to be Here, the dimensionless parameter = * /(K * A * 0 ) represents the stiffness of the cell's cortex relative to its bulk; the preferred cell perimeter L 0 = L * 0 / A * 0 is often expressed in terms of a second dimensionless parameter = −2 L 0 . Major features of the ( , ) parameter space are shown in Fig. 1(b). 1. (a) Vertex model representation of tissue geometry. The centroid of cell α is located at R α relative to the fixed origin O. The position of cell vertices R i α are given relative to the centroid of a cell. Each vertex has three vectors, R i α ,R i α ,R i α , pointing to the same vertex from cells, α,α ,α . Cell properties, such as area and tangents along edges, are also defined relative to the cell centroid. (b) ( , ) parameter space, showing regimes in which tissues exhibit distinct behavior. Following Ref. [11], region I is a "soft" network with no shear resistance; the network becomes solidlike in region II. For hexagons, P eff 6 = 0 [see (16)] has a single positive root in region II a and two positive roots in region II b . The network collapses in region III. Circular markers indicate the locations of parameters used for simulations in Fig. 4 A rigidity transition characteristic of a glassy material takes place along L 0 = μ 6 for regular hexagons [11,22], where μ Z = 2[Z tan(π/Z)] 1/2 (a regular Z-gon has an exact perimeter area relationship L = μ Z √ A). For a disordered monolayer, the transition occurs along L 0 ≈ μ 5 [31]. The transition between the fluid regime (I) and the solid regime (II) is indicated in Fig. 1(b). We avoid region III below, where the network collapses.
We label derivatives ∂U α /∂A α and ∂U α /∂L α of (3) as a pressure and a tension, respectively, we can then interpret −f i α as the elastic restoring force generated by cell α when vertex i undergoes a small displacement.
In a departure from many previous models we introduce the cell's (dimensionless) energy dissipation rate as where a dot above a variable denotes a time derivative. This accounts for viscous dissipation associated with the shape changes of individual cells, and is expressed in terms of the two geometric variables characterizing cell shape, A α and L α , that appear in the strain energy U α ; this model does not describe frictional interactions with a substrate. It follows from (1) thaṫ The parameters γ and μ can be related to their dimensional counterparts via γ = γ * /K * T * and μ = μ * L * 2 0 /K * T * A * 0 through a choice of time scale T * that we do not specify immediately. It follows from (5) that α =  [32], who treated the analogous 1D problem, we minimize the total dissipation rate across the monolayer, = α α , subject to a constraint ensuring the dissipation of total mechanical energy U = α U α through , This is achieved by minimizing the Lagrangian L = + ζ (U + ) for some Lagrange multiplier ζ . The first variation of L with respect to the velocity of each vertex in M must vanish, i.e., where for each j the sum is over the three cells adjacent to R j ; likewise, L ζ = 0 yields (6). Acting on (7) with jṘ j · yields (1 + ζ )2 + ζU = 0, which with (6) implies ζ = −2 and hence (7) gives the net force balance on each vertex as 2F j = 0 (j = 1,2, . . . ,N v ), where F i α can be interpreted as the viscoelastic restoring force due to cell α alone following a small displacement of its ith vertex.
It is computationally convenient (particularly when modeling a viscous interaction with a substrate) simply to impose a drag on each vertex, leading to an explicit set of ordinary differential equations (ODEs, of the formṘ j ∝ − {R i α =R j } f i α ) that can be used to step the system forward in time; in contrast, (8) couples time derivatives in a more complex manner, falling into a class of models reviewed in Ref. [14]. The tradeoff is a formulation that combines elastic responses to area and perimeter changes [via (4)] with their natural viscous counterparts, leading to an expression for the cell stress tensor (2) and (8)] This retains the property that the principal axes of cell stress and cell shape (as defined by a shape tensor based on vertex locations) are aligned [20].

B. Tissue-level stress
The stress tensor for a single cell (9) may be rewritten as where J α = 1 2 I − Q α is deviatoric, satisfying Tr(J α ) = 0. The isotropic component of the stress is given in terms of the effective cell pressure, defined as Positive (negative) values of P eff α indicate that the cell is under net tension (compression).
Assuming the monolayer forms a simply connected region of tissue, the tissue-level stress σ M satisfies [20] where the sum is over all cells in M and A M = α A α . Correspondingly, the isotropic component of tissue-level stress is expressed in terms of the effective tissue pressure For a monolayer under isotropic external loading, the deviatoric component of the global stress must vanish. Thus, once in equilibrium, the system satisfies P eff = P ext , where P ext is the peripheral pressure, assumed uniform. An isolated monolayer under conditions of zero external loading must satisfy P ext = 0, i.e., allowing cells to be grouped into those that are under net tension (P eff α > 0) and net compression (P eff α < 0). The deviatoric component of the tissue-level stress, has eigenvalues ±ξ , where ξ = det (σ M ) − Tr(σ M ) 2 . These quantify the tissue-level shear stress, and provide a measure of the spatial anisotropy of the monolayer [33]. The expression for tissue-level stress σ M = −P eff I +σ M using (13) and (15) extends results derived in Refs. [19,20,34] to account for viscous resistance to area and perimeter changes, although it does not include additional dissipative stresses associated with neighbor exchanges or extrusion of very small cells. An alternative derivation of σ M directly from U and is provided in Appendix A, where we show that σ M :Ė = −(U + ) = 0 for a small-amplitude homogeneous strain E; larger deformations, leading to neighbor exchanges, are therefore needed in order to change the system's internal energy. The equilibria reported below are, however, unaffected by the choice of dissipation in that they are always local minima of U (R 1 , . . . ,R N v ).

C. Simulation methodology
Simulations were generated using the methodology outlined in Ref. [20]. Initial distributions of cell centers were generated using a Matérn type II random sampling process. Cell edges and vertices were formed by constructing a Voronoi tessellation about the seed points, imposing periodic boundary conditions in a square domain. The system was then relaxed towards the nearest energy minimum. Following the initialization of the disordered geometry, a series of isotropic expansions or contractions was imposed until (14) was satisfied within a prescribed tolerance. Stretching deformations were imposed by mapping all vertices and the domain boundary by an affine transformation and then allowing the system to relax.
The effective pressure (11) of a regular hexagon at equilibrium is given by where A is the area of the hexagon. We define A * 6 ( , ) to satisfy P eff 6 (A * 6 ) = 0, where the hexagon is stress free, identifying A * 6 as a length scale. During relaxation, T1 transitions (neighbor exchanges) were performed on edges with a length less than 0.1 A * 6 . Three-sided cells with an area less than 0.3A * 6 were removed via a T2 transition (extrusion). Fletcher et al. [13] and Spencer et al. [35] outline refined treatments of these deformations. Our focus is primarily on the mechanical properties of a monolayer across region II in Fig. 1. For comparison with experiments we adopt parameters fitted previously to our experimental system [20], namely, ( , ) = (−0.259,0.172), acknowledging the imperfection of the fit and the inherent challenges of parameter estimation in this system [36].

D. Experimental methods
Our Xenopus embryonic animal cap preparation, stretch assay, and imaging protocol are described in Appendix B. Briefly, explants from stage 10 Xenopus laevis embryos were cultured on a fibronectin-coated poly(dimethylsiloxane) (PDMS) membrane. The tissue layer is three cell layers thick; the basal cells were attached to the membrane while the apical cells were imaged with confocal microscopy. The apical cells were not in direct contact with the PDMS membrane, shielding them from influences such as substrate-mediated integrin activation or focal adhesion formation [37]. Uniaxial in-plane stretching of the rectangular PDMS membrane (with a strain imposed on two opposite lateral boundaries and no stress imposed on the other two) deformed the tissue layer, with strains transmitted from the membrane to the apical layer via the basal cells. Using green fluorescent protein (GFP)-α-tubulin cell cortex and cherry-histone nuclear markers, the images of the apical cells were manually segmented and cell boundaries were approximated as polygons using a PYTHON script. The effective tissue pressure P eff of the unstretched monolayer was presumed to be zero, allowing the preferred area parameter A * 0 to be calculated following Ref. [26] using fitted parameters ( , ) = (−0.259,0.172) from Ref. [20]. A * 0 was assumed to be unchanged when the epithelium was stretched, allowing P eff α of individual cells to be estimated. The principal axis of stress for individual cells was identified using the shape tensor based on each cell's vertex locations, as described in Ref. [20].

III. PERTURBING TISSUE STRUCTURE USING DEFORMATIONS
A. Stretching embryonic epithelium membrane; the finite thickness of the tissue means that the apical layer experiences a lower level of strain. The isotropic distribution of apical cells in the undeformed state is disrupted by stretching: Cells under net tension (darker) tend to align their principal axis of stress with the direction of the global stretch, whereas cells under net compression (lighter) tend to align their principal axis of stress vertically [ Fig. 2(d)]. Stretch therefore induces anisotropy and a degree of order to the monolayer, consistent with previous observations [19,26,27,38,39].

B. Simulated tissue stretching
We mimic these observations using simulations, seeking to characterize the mechanical properties of the deformed monolayer. We describe the passive response to stretch, ignoring additional mechanosensitive effects of the kind described at the single-cell level by Xu et al. [40]. Figure 3 shows the result of performing a 20% area-preserving stretch on a simulated monolayer, with the mapping of vertices where λ = 0.2, followed by relaxation to equilibrium via T1 transitions. In the undeformed state, the cell orientations are approximately isotropic [ Fig. 3(c)], and roughly equal proportions of cells are under net tension and net compression [ Fig. 3(a)]. Stretching increases the proportion of cells under net tension [ Fig. 3(b)], inducing a striking degree of orientational order. As in Fig. 2(d), cells under net tension tend to align their principal axis of stress with the direction of stretch, whereas cells under net compression tend to align their principal axis of stress perpendicularly [ Fig. 3(d)]. The ordering is more extreme than in experiments; it is likely that for this large deformation, refinements of the functional form of the mechanical energy U [leading to linear pressure and tension relations in (4)] are needed to make the model quantitatively accurate.
In addition, the effective tissue pressure P eff , measuring the isotropic component of stress in the monolayer (11), increases with the degree of stretch [ Fig. 4(a)]. P eff is not significantly affected by the number of steps in which stretching is done, even when the monolayer is relaxed after every step [ Fig. 4(a)]. In the linear regime (λ 1), (17) is a pure shear deformation, leading to a negligible increase in P eff for small stretches; however, a nonlinear response emerges at larger amplitudes. Figure 4(b) demonstrates that the tissue shear stress ξ [see (15)], induced by individual stretches of increasingly large amplitude from an unstressed isotropic initial condition (with subsequent relaxation via T1 transitions), increases with stretch magnitude. This holds for both the monolayer given in Fig. 3 and an example closer to the region I/II boundary [ Fig. 1(b)] where ( , ) = (−0.569,0.145). The latter has a significantly smaller slope, indicating much lower resistance to shear in this region of parameter space. For reference, the anisotropy in these examples measured by a more traditional order parameter is illustrated in Appendix C. The present model does not account for stress relaxation that may arise under large strains via cell division, as illustrated in simulations by Ref. [39].
Given that stretching induces a strong degree of anisotropy in the tissue, we now assess how this ordering influences further tissue deformation. An initially isotropic tissue is uniaxially stretched via (17) and relaxed to equilibrium, as above. This deformed (i.e., prestretched) configuration has some prestress σ 0 . A further small-amplitude homogeneous strain I → I + E changes the global stress as σ 0 → σ 0 + σ M . We evaluate the perturbation stresses σ M arising from unidirectional strains E X = diag(λ,0) and E Y = diag(0,λ) directly from (12), noting that each deformation combines expansion and shear, e.g., It follows that both the shear and bulk elasticity of the tissue will contribute to the induced perturbation stress. Figure 5 plots two components of the perturbation stress against the magnitude of prestretch, for the tissue shown in Fig. 3(a), when subjected to additional 1% strains in the x and y directions. When the tissue is subject to a weak FIG. 4. (a) The effect of incremental stretch on effective tissue pressure. The tissue shown in Fig. 3(a) was subjected to a 30% area-preserving uniaxial stretch in a varying number of steps (straight lines are drawn between data points). The total stretch was divided into equally spaced increments and the tissue was relaxed after every stretch. The tissue starts at P eff = 0 and ends at approximately P eff = 0.57, regardless of how many steps were used. Translucent shading indicates 95% confidence intervals over the five simulations. (b) Shear stress ξ vs magnitude of stretch with ( , ) = (−0.259,0.172) (solid; the monolayer given in Fig. 3) and ( , ) = (−0.569,0.145) (dashed). See Fig. 1(b) for locations in parameter space. Each data point represents an instantaneous stretch performed on the same initial isotropic monolayer satisfying P ext = 0. The monolayers were relaxed to equilibrium following stretch. initial stretch, subsequent perturbation stresses are almost , indicating that the tissue is mechanically isotropic. However, under increased prestrain, we see that σ M xx (E X ) > σ M yy (E Y ), indicating increased anisotropy. The figure demonstrates how the tissue can be preferentially stiffened in one direction by an imposition of prestretch. Correspondingly, of the membrane elements shown in Fig. 3(b), a greater net membrane length is oriented in the x direction (resisting further stretch) in comparison to the y direction.

IV. DERIVING THE STIFFNESS TENSOR FOR ANISOTROPIC TISSUES
We have demonstrated that anisotropic deformations induce ordering in a tissue, leading to an anisotropic response to external loading reminiscent of an orthotropic material. We now evaluate directly the stiffness and damping tensors C(σ 0 ) and D(σ 0 ) appearing in the linearized relation connecting perturbation strain and perturbation stress σ = −(C : E + D : E) for some imposed strain E. This small-amplitude relation does not account directly for stress relaxation via neighbor exchanges, although these can influence the prestress σ 0 of the base state.

A. The perturbation stress of a single cell
We begin by considering a single cell with stress given by (9). From a stationary state with prestress we impose a small-amplitude, homogenous, symmetric, timevarying strain E, such that position vectors transform as Linearizing in the small strain amplitude yields the mappings (Appendix A) where B α is a fourth-order tensor (A4) satisfying so that Tr(B α : E) = Q α : E. Dynamic area and perimeter changes are coupled to E(t) viaȦ α = A α Tr(Ė) andL α = L α Q α :Ė.
Writing the stress following the deformation as σ α = σ o α + σ α , we use (20) and (9) to give, for a monolayer at equilibrium, to first order. Including time-dependent terms, the perturbation stress is therefore We can separate the isotropic and deviatoric contributions as where P eff α is the perturbation effective pressure (25) and the traceless contribution characterizing perturbation shear is 052409-7

B. The perturbation stress of the tissue
Applying E to the entire monolayer, the global stress transforms as σ M → σ M + σ M so that (12) becomes neglecting terms quadratic in E. Linearizing the left-hand side of (27), the global perturbation stress is given by Thus the effective perturbation tissue pressure is The predictions arising from (28) in prestretched monolayers being subjected to small-amplitude strains are tested in Fig. 5, showing good agreement with direct stress computations. Thus, for a given E, the stiffness matrix C can be evaluated directly from the terms proportional to E in (28) and its viscous analog D from terms proportional toĖ, using the notation {A ⊗ B} ij kl = A ij B kl . We now illustrate this in the special case of an initially unstressed disordered monolayer.

C. Elastic moduli for a disordered isotropic monolayer
When the base state is a disordered isotropic monolayer at zero stress [satisfying (14) withσ M = 0], we can derive bulk moduli by imposing a small isotropic expansion with E = λI, where λ 1. The deformation satisfies Under an isotropic load, the deviatoric components of the perturbation stress vanish ( N c α=1 T α L α J α = 0). Using (31), the bulk perturbation effective pressure (29) is The bulk and cortical viscosities [appearing in the coefficient ofλ in (32)] contribute in a similar manner to P eff as the bulk and cortical stiffnesses, except via a nonlinear dependence on L α . Recalling that A M = A M Tr(E), the bulk elastic modulus K e can be derived from (32) withλ = 0 as in agreement with Ref. [20]. Figure 6(a) demonstrates how K e varies across parameter space, with the tissue becoming less resistant to isotropic deformation towards the region II/III boundary. The tissue stress arises through the area weighting of cellular stress (12), leading to a nonlinear area dependence of the bulk modulus on the cell area in (33); thus when cells are substantially smaller than their target area (near the II/III boundary), the bulk modulus falls accordingly, approaching near-zero values.
For the shear moduli, we impose a small simple-shear deformation with E = κe x e y , where |κ| 1 and e x = (1,0), e y = (0,1) are the Cartesian coordinate bases, and seek σ xy . This simple deformation satisfies Tr(E) = 0. To evaluate Q α : E and B α : E, we writet i α = cos θ i α e x + sin θ i α e y , where θ i α satisfiest i α · e x = cos θ i α andt i α · e y = sin θ i α . Then, Similarly, noting that E ·t i α = κ sin θ i α e x andt i α · E ·t i α = (κ/2) sin (2θ i α ), we have Thus from (28) we have To evaluate the elastic shear modulus of the disordered monolayer, we set σ M xy = −κG e withκ = 0. The shear modulus is given by Equation (37) recovers previous predictions for the shear modulus of periodic monolayers, where all cells are perfect hexagons (L 2 α = 8 √ 3A α , all edges have equal length, θ i α = 2πi/6, and the terms with sums over cos and sin vanish) [16,20]; however, it extends these results by allowing the direct evaluation of the shear modulus for a disordered monolayer. Figure 6(b) demonstrates how G e , as predicted by (37), varies across parameter space. Interestingly, the tissue becomes less resistant to shear as it becomes increasingly resistant to isotropic deformations. For comparison, Fig. 6(c) shows the computationally simulated shear modulus, directly evaluated from the global perturbation stress tensor as σ M 12 = −κG e , following a 1% simple-shear deformation (κ = 0.01) on the simulated monolayers. There is excellent agreement between the analytic and simulated results.
For a sufficiently large disordered but isotropic monolayer, we might assume that the terms with sums over cos and sin in (37) vanish when summed over all cells [the degree of anisotropy can be assessed with (15)]. The elastic shear modulus for an isotropic monolayer is then approximated by showing how G iso e falls to zero as the tension in each cell approaches zero. The percentage difference between G e and G iso e in example simulated isotropic tissues with 800 cells is shown in Fig. 6(d), showing close agreement (< 2% relative difference) almost everywhere across region II. Discrepancies arise only close to the region I/II boundary, where G e and G iso e both approach zero. It is notable that the dynamic shear resistance term μĖ in (36) is discarded under the approximation that leads to (38).

D. Composite isotropic and shear deformations
Finally, we recall that the perturbation stress tensor has eigenvalues σ M ± = − P eff ± ξ , where (− P eff , − P eff ) and ( ξ, − ξ ) are the eigenvalues of the isotropic and deviatoric (shear) components of (28) respectively. Figure 7 shows how P eff / ξ varies across parameter space for isotropic monolayers subjected to static strain E X [see (18)] of amplitude 0.01, which induces both an isotropic and deviatoric stress response. The parameter space partitions into a region that is more resistant to area change (region A) and one that is more resistant to shear (region B). Figure 7 highlights how the isotropic stress, resisting area change, falls dramatically near the region II/III boundary [see Fig. 1(b)].

V. DISCUSSION
Previous reports have found that internal patterning in tissues can be linked to the mechanical properties of the material [41][42][43]. We find that stretching induces ordering within the tissue, with cells being elongated on average in the direction of stretch, consistent with previous observations [19,27,38,39]. Inferring relative stresses using the vertex-based model provides additional insight: Distinguishing cells that are under net tension (with positive isotropic stress P eff α > 0) from those under net compression, we find strong alignment of the former with the direction of stretch [ Fig. 2(c)] and of the latter with the perpendicular direction, in response to the imposed compressive stress, while retaining heterogeneity at the single-cell level. This feature emerges strongly in direct simulations also (Fig. 3). Increased spatial organization of cells is associated with anisotropic mechanical properties, which we characterized by deriving an explicit tissue-level stress/strain/strain-rate relationship (25), (26), and (28) describing the response of a prestressed tissue to small-amplitude homogeneous deformations.
The stress tensor we employed builds on the formulation derived by Nestor-Bergmann et al. [20] and others [44,45], neglecting nonplanarity [46,47], curved cell edges [18,48], and further refinements but including internal dissipation due to dynamic area and length changes of individual cells in a way that naturally complements the assumed strain energy. Our formulation ensures no net change in internal energy under a homogeneous deformation (Appendix A) and is suited to describing the viscoelastic properties of freely suspended monolayers, as described by Ref. [27]. The model is in the spirit of, but differs from, that of Okuda et al. [49], who proposed a drag force depending on an average of nearby vertex velocities. The linearized stress/strain relationship (25), (26), and (28) does not include additional dissipative effects of substrate drag or irreversible cell rearrangements. A framework for including additional plastic stresses and strains has been proposed within a coarse-grained model [10]; simulations of large-amplitude plastic tissue deformations under external loading using a discrete cell model are provided by Ref. [9].
Under the present vertex model, the perturbation stress of a prestressed tissue is given by the area-weighted sum of perturbation stresses of the individual cells (28). This leads to an expression for the fourth-order stiffness tensor C (30a), which describes how anisotropic tissues resist deformation through reversible elastic deformations, and its viscous analog D (30b). This formulation extends previous approaches to upscaling the vertex-based model in spatially periodic [16] and disordered isotropic networks [20]. Exact expressions for elastic moduli for a given monolayer realization are provided as explicit sums [for isotropic monolayers, these are (33) and (37)]. Further work is required to derive a priori predictions for the behavior of these quantities over multiple monolayer realizations. However, a crude simplification of the estimate of the elastic shear modulus (38) for a disordered isotropic monolayer is accurate across the bulk of region II of parameter space, but not close to the phase transition along the region I/II boundary where the shear modulus approaches zero. Our predictions can be compared with those of Merzouki et al. [17], who used simulations to impose stress and measure strain of a periodic hexagonal monolayer in order to infer bulk elastic parameters. Our results are broadly consistent with theirs, including evidence of nonmonotonic behavior close to the region II/III boundary (revealed most clearly in the stress ratios plotted in Fig. 7). Likewise, Xu et al. [39] report reduced stress under uniaxial strain for increasing and , mirroring the bulk elastic response in Fig. 6(d). By comparing the relative size of isotropic and shear stresses induced by a strictly uniaxial deformation, the vertex model can also be used to partition parameter space into regions of higher shear modulus and lower bulk modulus, and vice versa (Fig. 7). This may be relevant to phases during development when tissues undergo extreme shape changes and may have a bearing on the different mechanical environments of epithelial tissues in various organ systems.
In summary, this study demonstrates how loading organizes the cell-scale stress field in a stretched monolayer and how mechanical viscoelastic moduli of disordered or anisotropic cellular monolayers can be determined as explicit sums over cells. Further steps towards deriving well-grounded homogenized descriptions of such media will require assessment of the statistical distributions of different cell classes over the plane.

APPENDIX A: A HOMOGENEOUS SMALL-AMPLITUDE STRAIN
We derive the mappings of key geometric quantities under a small-amplitude homogeneous and symmetric strain E, which transforms position vectors as R → R + E · R, and then apply these mappings to the mechanical energy U .
We assume that all quantities X follow a mapping of the form X → X + X, defined relative to the same cell, and therefore temporarily drop the α subscript. Tangents are defined as The length of an edge is given by l i = (t i · t i ) 1 2 . To linear order, demonstrating that l i = l i (t i · E ·t i ). The cell perimeter is Note thatt i · ( t i ) = 0, ensuring that unit vectors are rotated but not stretched.
We now consider the deformation LQ → LQ + (LQ), writing to linear order. Thus we see that (LQ) = LB : E, where B in component form is To evaluate area changes, we begin by considering the area of a subtriangle, = {R α ,R i α ,R i+1 α }, comprising the cell centroid and two vertices from an edge (Fig. 8). We have used the α subscript notation for clarity in defining the centroid, but drop it again from this point. The area of this triangle is to linear order. Since E is symmetric, we can make use of the spectral theorem to write E = λ 1 e 1 e 1 + λ 2 e 2 e 2 , where (e 1 ,e 2 ) form an orthonormal basis of eigenvectors. Inserting into (A5) we have Defining θ i and θ i+1 such that λ 1 R i × e 1 = −λ 1 |R i | sin θ iẑ , and e 1 · R i+1 = |R i+1 | cos θ i+1 (see Fig. 8) gives Summing over subtriangles gives A = Tr(E)A. Having shown that under the deformation R → R + E · R, lengths and areas of individual cells transform according to L α → L α (1 + Q α : E), A α → A α (1 + Tr(E)), it is straightforward to determine the associated change in recoverable mechanical energy (3) as α . Here, we have decomposed the stress in (9) into its elastic and viscous components σ α = σ (e) α + σ (v) α . The sign change between σ and arises from the difference between the stresses exerted on, or by, a cell or tissue. Similarly defining (v) ≡ − A α σ (v) α , and noting that for this small deformationȦ α = A α Tr(Ė), it follows that where = α α is the dissipation rate (5). In the absence of neighbor exchanges, which would contribute additional stresses and deformations (treated in a coarse-grained approximation by Ishihara et al. [10]), the total rate of change of internal energy of the system is thereforeĖ = :Ė = ( (e) + (v) ) :Ė =U + = 0 [by (6)]. Positing the thermodynamic relation E = T S + U , where S is an entropy change at temperature T , it follows that for the imposed deformation E, TṠ = −U = 0.
Animal cap tissue was dissected from the embryo at stage 10 of development (early gastrula stage) following a protocol previously described by Joshi and Davidson [53], and cultured on a 20 mm × 20 mm × 1 mm elastomeric PDMS membrane coated with fibronectin (incubated at 4 • C overnight with 1 ml of 10 μg/ml fibronectin). The fibronectin was removed and each membrane was subsequently washed three times with 1× phosphate buffered saline (PBS) followed by two washes with Danilchik's for Amy explant culture media (DFA; 53 mM NaCl 2 , 5 mM Na 2 CO 3 , 4.5 mM potassium gluconate, 32 mM sodium gluconate, 1 mM CaCl 2 , 1 mM MgSO 4 ) prior to the introduction of each animal cap. Animal cap explants were excised using forceps and hair knives to make neat squares of tissue. Each explant was transferred to a PDMS membrane filled with DFA and a coverslip with vacuum grease at each end was placed over the top to ensure the explant adhered to the fibronectin-coated membrane. Each membrane was then incubated at 18 • C for at least 2 h prior to imaging.
Each PDMS membrane was attached to a stretch apparatus (custom made by Deben UK Limited) fixed securely to the stage of a Leica TCS SP5 acousto-optical beam splitter (AOBS) upright confocal and a 0.5 or 8.6 mm uniaxial stretch was applied for control (unstretched) and stretched samples, respectively. Images were collected on the Leica TCS SP5 AOBS upright confocal using a 20×/0.50 HCX Apo U-V-I [W (dipping lens)] objective and 2× (or 1×) confocal zoom. The confocal settings were as follows: pinhole 1 Airy unit, scan speed 400 Hz bidirectional, format 512 × 512 (or 1024 × 1024). Images were collected using hybrid detectors with the following detection mirror settings: HyD2 fluorescein isothiocyanate (FITC) 500-550 nm; HyD4 Texas red 590-690 nm using the 488 nm (20%) and 543 nm (100%) laser lines, respectively. Images were collected sequentially. The distance between each optical stack was maintained at 4.99 μm and the time interval between each capture was restricted to 20 s, with each sample being imaged for up to 2.5 h. Only the maximum intensity projections of these 3D stacks are shown in the results.

APPENDIX C: ORDER PARAMETER
A more traditional measure of spatial disorder is provided by the order parameter Q = cos 2θ , where θ measures the angle of the principal axis of the shape tensor of each cell ( Z α −1 i=0 R i α ⊗ R i α ) with respect to the direction of stretch and the average is taken over all the cells in the monolayer. Figure 9 illustrates the evolution of Q for the stretch realizations illustrated in Fig. 4(b). It is notable that while the two examples are quite similar by this geometric measure, the shear stress [ Fig. 4(b)] is substantially lower for parameters closer to the region I/II boundary in parameter space [ Fig. 1(b)].
FIG. 9. The order parameter Q, for the realizations plotted in Fig. 4(b).