Mechanical characterisation of disordered and anisotropic cellular monolayers

Mechanical characterisation of disordered and anisotropic cellular monolayers

Alexander Nestor-Bergmann [ Wellcome Trust Centre for Cell-Matrix Research, School of Medical Sciences, Faculty of Biology, Medicine & Health,
Manchester Academic Health Science Centre, University of Manchester, Oxford Road, Manchester M13 9PT UK
   Emma Johns Wellcome Trust Centre for Cell-Matrix Research, School of Medical Sciences, Faculty of Biology, Medicine & Health,
Manchester Academic Health Science Centre, University of Manchester, Oxford Road, Manchester M13 9PT UK
   Sarah Woolner Wellcome Trust Centre for Cell-Matrix Research, School of Medical Sciences, Faculty of Biology, Medicine & Health,
Manchester Academic Health Science Centre, University of Manchester, Oxford Road, Manchester M13 9PT UK
   Oliver E. Jensen School of Mathematics, University of Manchester, Manchester M13 9PL, UK
July 7, 2019

We consider a cellular monolayer, described using a vetex-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 neighbour 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.

Cell stress, cell shape, vertex model, tissue mechanics, upscaling.
preprint: PREPRINT

Also at ]School of Mathematics, University of Manchester, Manchester M13 9PL, UK

I Introduction

Epithelial tissues have significant roles in embryonic development, tissue homeostasis and disease development (Guillot and Lecuit, 2013). Recent work has revealed that many critical functions in biological tissues are dependent on the accurate organisation of the shapes and packing geometry of the constituent cells (Lecuit and Lenne, 2007). Disturbances in this organisation have been associated with problems during embryonic development and diseases in adult life (Carney and Jacobs, 1984; Jiao et al., 2011). Furthermore, there is evidence that mechanical forces may directly trigger biochemical responses to regulate morphogenetic processes (Schluck et al., 2013; Heisenberg and Bellaïche, 2013). However, due to difficulties in quantifying stresses in tissues, the mechanisms by which tissue behaviour emerges from these multiscale feedback processes remain poorly understood.

Continuum descriptions can provide useful insights to tissue-level behaviour. For example, elastic-viscoplastic continuum models can capture the solid- and liquid-like response of tissues to small and large deformations over differing timescales (Preziosi et al., 2010; Weickenmeier and Jabareen, 2014). However, they are in many cases not built from explicit physical description of cells. Furthermore, the interactions of multiple cells can lead to rich emergent behaviour at the tissue scale, such as yielding and remodelling, that is not easily accessed through conventional continuum frameworks (Pathmanathan et al., 2009; Ishihara et al., 2017a).

Discrete vertex-based models of epithelia have proven to be a useful tool in linking mechanics to tissue morphology (Farhadifar et al., 2007; Nagai and Honda, 2001; Fletcher et al., 2014; Alt et al., 2017; Fletcher et al., 2017). These models have more recently been developed to characterise the mechanical properties of tissues (Murisic et al., 2015; Merzouki et al., 2016) and to infer local and global stresses (Brodland et al., 2014; Sugimura and Ishihara, 2013; Nestor-Bergmann et al., 2017a). 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 (Nestor-Bergmann et al., 2017a; Yang et al., 2017). Relationships between discrete models and traditional continuum approaches have been found for spatially periodic cell networks (Murisic et al., 2015; Staple et al., 2010a), while equivalent relationships for disordered tissues have only been partially established for isotropic disordered cell networks (Nestor-Bergmann et al., 2017a) or for analogous physical systems such as two-dimensional (2D) dry foams Kruyt (2007).

Many mechanical models of biological tissues assume that the material is isotropic. However, recent observations in the Drosophila melanogaster embryo (Gao et al., 2016) have provided evidence that biological tissues may exhibit orientational, as well as positional, structure. Likewise, models of 3D foams have explored how orientational structure can be introduced through uniaxial deformations (Evans et al., 2017). The deformation induces a net stress in the material, leading to a series of irreversible deformations (such as neighbour exchanges). Experimental observations of epithelial tissues have revealed similar patterns of orientational order following stretching (Sugimura and Ishihara, 2013; Nestor-Bergmann et al., 2017b).

To explore how deformation induces anisotropy, in this paper we use a variant of a well-studied vertex-based model to quantify the mechanical behaviour 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 timescale of minutes Doubrovinski et al. (2017), evidence of shorter but distinct stress-relaxation timescales in the cytosol and cortex Bambardekar et al. (2015) and sub-cellular observations of dissipation at cell contacts Clément et al. (2017). 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 behaviour 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 characterising the perturbation stress of a pre-stressed 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 (Farhadifar et al., 2007; Fletcher et al., 2014; Bi et al., 2015), using the notational framework presented in Nestor-Bergmann et al. (2017a). Details of our experimental protocol follow a summary of the model.

ii.1 The vertex-based model

The epithelial monolayer, , is represented as a spatially disordered planar network of vertices, labelled , connected by straight edges and bounding convex polygonal cells, labelled . The cells are assumed to have identical mechanical properties described in terms of a preferred area , a preferred perimeter , a bulk stiffness , a cortical stiffness , a bulk viscosity and a cortical viscosity . Scaling all distances on , the vector from the coordinate origin to vertex is given by . Each vertex is shared by three cells and edges are shared by two cells (excluding cells at the boundary of ). Quantities specific to cell are labelled by a Greek subscript and defined relative to its centroid (Figure 1A).

Figure 1: A. Vertex model representation of tissue geometry. The centroid of cell is located at relative to the fixed origin, . The position of cell vertices, , are given relative to the centroid of a cell. Each vertex has three vectors, , 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 behaviour. Following (Farhadifar et al., 2007), region I is a ‘soft’ network with no shear resistance; the network becomes solid-like in region II. For hexagons, (see (16)) has a single positive root in region II and two positive roots in region . The network collapses in Region III. Circular markers indicate locations of parameters used for simulations in Figure 4B below.

Cell has vertices labelled anticlockwise by relative to . We define as the vector from the cell centroid to vertex , such that . Anticlockwise tangents are defined by (taking modulo ), unit vectors along a cell edge by and outward normals to edges by (where is a unit vector pointing out of the plane). The length, , of an edge belonging to cell between vertices and , the cell perimeter, , and the cell area, , are given by


We note for later reference that and where ; furthermore Nestor-Bergmann et al. (2017a),


is a symmetric tensor characterising the shape of cell , satisfying .

The dimensionless mechanical energy of an individual cell, (scaled on ), is assumed to be


Here the dimensionless parameter represents the stiffness of the cell’s cortex relative to its bulk; the preferred cell perimeter, is often expressed in terms of a second dimensionless parameter . Major features of the -parameter space are shown in Figure 1B. A rigidity transition characteristic of a glassy material takes place along for regular hexagons (Farhadifar et al., 2007; Staple et al., 2010b), where (a regular -gon has an exact perimeter area relationship, ). For a disordered monolayer, the transition occurs along (Bi et al., 2015). The transition between the fluid regime (I) and the solid regime (II) is indicated in Figure 1B. We avoid region III below, where the network collapses.

We label derivatives and of (3) as a pressure and a tension respectively:


Setting , we can then interpret as the elastic restoring force generated by cell when vertex undergoes a small displacement.

In a departure from many previous models we introduce the cell’s (dimensionless) energy dissipation rate as


where a dot denotes a time derivative. This accounts for viscous dissipation associated with shape changes of individual cells. It follows from (1) that , , so that and . The parameters and can be related to their dimensional counterparts via and through a choice of timescale that we do not specify immediately. It follows from (5) that and . Following Fozard et al. Fozard et al. (2009), 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 through ,


This is achieved by minimizing the Lagrangian for some Lagrange multiplier . The first variation of with respect to the velocity of each vertex in must vanish, i.e.


where for each the sum is over the three cells adjacent to ; likewise yields (6). Acting on (7) with yields , which with (6) implies and hence (7) gives the net force balance on each vertex as () where


can be interpreted as the viscoelastic restoring force due to cell alone following a small displacement of its th vertex.

It is computationally convenient (particularly when modelling viscous interaction with a substrate) simply to impose a drag on each vertex, leading to an explicit set of ODEs (of the form ) 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 Alt et al. (2017). 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 of the form (using (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 (Nestor-Bergmann et al., 2017a).

ii.2 Tissue-level stress

The stress tensor for a single cell (9) may be rewritten as


where is deviatoric, satisfying . The isotropic component of the stress is given in terms of the effective cell pressure, defined as


Positive (negative) values of indicate that the cell is under net tension (compression).

Assuming the monolayer forms a simply connected region of tissue, the tissue-level stress, , satisfies Nestor-Bergmann et al. (2017a)


where the sum is over all cells in and . Correspondingly, the isotropic component of tissue level stress is expressed in terms of


For a monolayer under isotropic external loading, the deviatoric component of the global stress must vanish. Thus, once in equilibrium, the system satisfies , where is the peripheral pressure, assumed uniform. An isolated monolayer under conditions of zero external loading must satisfy , i.e.


allowing cells to be grouped into those that are under net tension () and net compression (). The deviatoric component of the tissue-level stress,


has eigenvalues , where . These quantify the tissue-level shear stress, and provide a measure of the spatial anisotropy of the monolayer Kraynik et al. (2003).

The expression for tissue-level stress using (13, 15) extends results derived in (Nestor-Bergmann et al., 2017a; Sugimura and Ishihara, 2013; Xu et al., 2016) to account for viscous resistance to area and perimeter changes, although it does not include additional dissipative stresses associated with neighbour exchanges or extrusion of very small cells. An alternative derivation of directly from and is provided in Appendix A, where we show that for a small-amplitude homogeneous strain ; larger deformations, leading to neighbour 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 .

ii.3 Simulation methodology

Simulations were generated using the methodology outlined in Nestor-Bergmann et al. (2017a). Initial distributions of cell centres 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 initialisation of the disordered geometry, a series of isotropic expansions or contractions were 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 is the area of the hexagon. We define , to satisfy , where the hexagon is stress free, identifying as a length scale. During relaxation, T1 transitions (neighbour exchanges) were performed on edges with length less than . 3-sided cells with area less than were removed via a T2 transition (extrusion). Fletcher et al. (2014) and Spencer et al. (2017) outline refined treatments of these deformations.

Our focus is primarily on mechanical properties of a monolayer across region II in Figure 1. For comparison with experiments we adopt parameters fitted previously to our experimental system Nestor-Bergmann et al. (2017a), namely , acknowledging the imperfection of the fit and the inherent challenges of parameter estimation in this system Kursawe et al. (2017).

ii.4 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 PDMS membrane. The tissue layer is three cell layers thick; the basal cells attached to the membrane while the apical cells were imaged with confocal microscopy. 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 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 global effective pressure of the unstretched monolayer was presumed to be zero, allowing the preferred area parameter to be calculated following Nestor-Bergmann et al. (2017b) using fitted parameters from Nestor-Bergmann et al. (2017a). was assumed to be unchanged when the epithelium was stretched, allowing 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 Nestor-Bergmann et al. (2017a).

Iii Perturbing tissue structure using deformations

iii.1 Stretching embryonic epithelium

Figure 2A shows apical cell boundaries, rendered as polygons, within a Xenopus embryonic epithelium in the unstretched configuration. Cell shapes are used to infer the relative isotropic stress as described in Sec. II.4. Cells estimated to be under net tension (red, with ) and under net compression (blue, with ) appear in roughly equal proportions, with their orientations distributed approximately isotropically (Figure 2C). The stress field is strongly heterogeneous at the single-cell level: cells under tension (compression) generate a primarily contractile (expansive) stress along their long (short) axis. Figure 2B shows the same epithelium after an application of an instantaneous 35% uniaxial stretch (in the horizontal direction) to the PDMS 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 (red) tend to align their principal axis of stress with the direction of the global stretch, whereas cells under net compression (blue) tend to align their principal axis of stress vertically (Figure 2D). Stretch therefore induces anisotropy and a degree of order to the monolayer, consistent with previous observations (Sugimura and Ishihara, 2013; Nestor-Bergmann et al., 2017b; Wyatt et al., 2015).

Figure 2: A. Epithelial apical layer of a Xenopus laevis animal cap, showing 801 cells rendered as polygons superimposed on the original microscopy image. for each cell was calculated assuming and . Line segments indicate the principal axis of shape and stress for each cell. Red (blue) cells have and exert a net inward (outward) force along each line segment. B. The apical layer in A following a 35% instantaneous uniaxial stretch (horizontal) of the membrane beneath the basal cells, resulting in an (95% confidence interval) uniaxial stretch of the apical cells. C, D. Histograms showing orientation of the principal axis of stress for cells under tension (red) and compression (blue), for apical layers given in A (corresponding to C) and B (corresponding to D). Purple colour is overlay of red and blue. Bin size was selected using Freedman–Diaconis rule.

iii.2 Simulated tissue stretching

We mimic these observations using simulations, seeking to characterize the mechanical properties of the deformed monolayer. Figure 3 shows the result of performing a area-preserving stretch on a simulated monolayer, with the mapping of vertices


where , followed by relaxation to equilibrium via T1 transitions. In the undeformed state, the cell orientations are approximately isotropic (Figure 3C), and roughly equal proportions of cells are under net tension and net compression (Figure 3A). Stretching increases the proportion of cells under net tension (Figure 3B), inducing a striking degree of orientational order. As in Figure 2D, 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 (Figure 3D). 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 (leading to linear pressure and tension relations in (4)) are needed to make the model quantitatively accurate.

Figure 3: A. A simulation of a representative monolayer realisation satisfying with 800 cells, for (see Figure 1 for location in parameter space). Cell colours and line segments follow the scheme in Figure 2. B. The monolayer in A following a 20% area-preserving uniaxial stretch and subsequent relaxation. C&D. Histograms showing orientation of the principal axis of stress for cells under tension (red) and compression (blue), for monolayers given in A (corresponding to C) and B (corresponding to D). Purple colour is overlay of red and blue. Bin size was selected using Freedman–Diaconis rule.

In addition, the area-weighted mean effective pressure , measuring the isotropic component of stress in the monolayer (11), increases with the degree of stretch (Figure 4A). is not significantly affected by the number of steps in which stretching is done, even when the monolayer is relaxed after every step (Figure 4A). In the linear regime (), (17) is a pure shear deformation, leading to negligible increase in for small stretches; however a nonlinear response emerges at larger amplitudes. Figure 4B demonstrates that the tissue shear stress, (see (15)), induced by individual stretches of increasingly large amplitude from an unstressed isotropic inital condition (with subsequent relaxation via T1 transitions), increases with stretch magnitude. This holds for both the monolayer given in Figure 3 and an example closer to the region I/II boundary (Figure 1B) where . The latter has a significantly smaller slope, indicating much lower resistance to shear in this region of parameter space.

Figure 4: A. Testing the effect of incremental stretch on tissue stress. The tissue shown in Figure 3A was subjected to a 30% area-preserving uniaxial stretch in a varying number of steps (indicated by multicoloured lines; straight lines are drawn between datapoints). The total stretch was divided into equally spaced increments and the tissue was relaxed after every stretch. The tissue starts at and ends at approximately , regardless of how many steps were used. Translucent shading indicates 95% confidence intervals over the 5 simulations. B. Shear stress, , versus magnitude of stretch with (solid; the monolayer given in Figure 3) and (dashed). See Figure 1B for locations in parameter space. Each datapoint represents an instantaneous stretch performed on the same initial isotropic monolayer satisfying . The monolayers were relaxed to equilibrium following stretch.

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. pre-stretched) configuration has some pre-stress . A further small-amplitude homogeneous strain, , changes the global stress as . We evaluate the perturbation stresses arising from unidirectional strains and 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 pre-stretch, for the tissue shown in Figure 3A, when subjected to additional 1% strains in the - and - directions. When the tissue is subject to weak initial stretch, subsequent perturbation stresses are almost equal (), indicating that the tissue is mechanically isotropic. However, under increased pre-strain, we see that , indicating increased anisotropy. The figure demonstrates how the tissue can be preferentially stiffened in one direction by imposition of pre-stretch. Correspondingly, of the membrane elements shown in Figure 3B, a greater net membrane length is oriented in the direction (resisting further stretch) in comparison to the direction.

Figure 5: Perturbation stress response to small amplitude deformations in a pre-stretched monolayer. The stretched monolayers used were the same as those for the solid line in Figure 4B, with and the magnitude of pre-stretch is indicated on the -axis. The equilibrium pre-stretched monolayers were subjected to small deformations in the - () and - () directions, with . The component of the perturbation stress tensor in the direction of stretch is indicated on the -axis, with giving the -directed stress following (green, lower line) and giving the -directed stress following (red, upper line). Solid lines indicate values directly evaluated using (12) and dashed lines are predicted values using (28).

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 and appearing in the linearised relation connecting perturbation strain and perturbation stress for some imposed strain . This small-amplitude relation does not account directly for stress relaxation via neighbour exchanges, although these can influence the pre-stress of the base state.

iv.1 The perturbation stress of a single cell

We begin by considering a single cell with stress given by (9). From a stationary state with pre-stress


we impose a small-amplitude, homogenous, symmetric, time-varying strain , such that position vectors transform as . Linearising in the small strain amplitude yields the mappings (Appendix A)


where is a fourth-order tensor (41) satisfying


so that . Dynamic area and perimeter changes are coupled to via and .

Writing the stress following the deformation as , 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 is the perturbation effective pressure


and the traceless contribution characterising perturbation shear is


iv.2 The perturbation stress of the tissue

Applying to the entire monolayer, the global stress transforms as so that (12) becomes


neglecting terms quadratic in . Linearising 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 pre-stretched monolayers being subjected to small-amplitude strains are tested in Figure 5, showing good agreement with direct stress computations. Thus, for a given , the stiffness matrix, , can be evaluated directly from the terms proportional to in (28) and its viscous analogue from terms proportional to . We now illustrate this in the special case of an initially unstressed disordered monolayer.

iv.3 Elastic moduli for a disordered isotropic monolayer

When the base state is a disordered isotropic monolayer at zero stress (satisfying (14) with ), we can derive bulk moduli by imposing a small isotropic expansion with , where . The deformation satisfies


Under an isotropic load, the deviatoric components of the perturbation stress vanish (). Using (30), the bulk perturbation effective pressure (29) is


The bulk and cortical viscosities (appearing in the coefficient of in (31)) contribute in a similar manner to as the bulk and cortical stiffnesses, except via a nonlinear dependence on . Recalling that , the bulk elastic modulus, , can be derived from (31) with as


in agreement with Nestor-Bergmann et al. (2017a). Figure 6A demonstrates how varies across parameter space, with the tissue becoming less resistant to isotropic deformation towards the region II/III boundary. The tissue stress arises through area-weighting of cellular stress (12), leading to a nonlinear area dependence of bulk modulus on cell area in (32); thus when cells are substantially smaller than their target area (near the II/III boundary) the bulk modulus falls accordingly, approaching near-zero values.

Figure 6: A. Heat map across discrete intervals of parameter space showing the value of the analytic bulk elastic modulus (calculated using (32)) of a disordered isotropic monolayer with 800 cells. B. Equivalent plot of the analytic shear modulus (calculated using (36)). C. Computationally estimated shear modulus across parameter space. The elastic shear modulus, , was estimated from the global perturbation stress as , following a 1% simple shear deformation () on the simulated monolayers used in B. Each datapoint is taken as an average from 5 realisations of a monolayer with 800 cells, with . D. Percentage difference () between the exact modulus (36) and the approximation (37), across parameter space; the difference is below 2% in the hatched region.

For the shear moduli, we impose a small simple shear deformation with , where and are the Cartesian coordinate basis, and seek . This simple deformation satisfies . To evaluate and , we write , where satisfies and . Then


Similarly, noting that and , we have


Thus from (28), we have


To evaluate the elastic shear modulus of the disordered monolayer, we set with . The shear modulus is given by


Equation (36) recovers previous predictions for the shear modulus of periodic monolayers, where all cells are perfect hexagons (, all edges have equal length, and the terms with sums over and vanish) (Nestor-Bergmann et al., 2017a; Murisic et al., 2015); however it extends these results by allowing the direct evaluation of the shear modulus for a disordered monolayer. Figure 6B demonstrates how , as predicted by (36), varies across parameter space. Interestingly, the tissue becomes less resistant to shear as it becomes increasingly resistant to isotropic deformations. For comparison, Figure 6C shows the computationally simulated shear modulus, directly evaluated from the global perturbation stress tensor as , following a 1% simple shear deformation () 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 and in (36) 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 falls to zero as the tension in each cell approaches zero. The percentage difference between and in example simulated isotropic tissues with 800 cells is shown in Figure 6D, showing close agreement ( relative difference) almost everywhere across region II. Discrepancies arise only close to the region I/II boundary, where and both approach zero. It is notable that the dynamic shear resistance term in (35) is discarded under the approximation that leads to (37).

iv.4 Composite isotropic and shear deformations

Finally, we recall that the perturbation stress tensor has eigenvalues , where and are the eigenvalues of the isotropic and deviatoric (shear) components of (28) respectively. Figure 7 shows how varies across parameter space for isotropic monolayers subjected to static strain (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.

Figure 7: Heat map across discrete intervals of parameter space showing the value of . Dashed line represents contour where , where tissues show dominant resistance to area change in region A and to shear in region B. A log-colour scale was used to help display the differences across a large range of values. The monolayers used for all heat maps were the same as those used in Figure 6, where each datapoint is taken as an average from 5 realisations of a monolayer with 800 cells, with .

V Discussion

Previous reports have found that internal patterning in tissues can be linked to the mechanical properties of the material (Donev et al., 2005; Majmudar and Behringer, 2005; Honnell et al., 1990). We find that stretching induces ordering within the tissue, with cells being elongated on average in the direction of stretch, consistent with previous observations (Sugimura and Ishihara, 2013; Wyatt et al., 2015). Inferring relative stresses using the vertex-based model provides additional insight: distinguishing cells that are under net tension (with positive isotropic stress ) from those under net compression, we find strong alignment of the former with the direction of stretch (Figure 2C) 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 (Figure 3). Increased spatial organisation of cells is associated with anisotropic mechanical properties, which we characterised by deriving an explicit tissue-level stress/strain/strain-rate relationship (25, 26, 28) describing the response of a pre-stressed tissue to small-amplitude homogeneous deformations.

The stress tensor we employed builds on the formulation derived by Nestor-Bergmann et al. (2017a) and others (Ishihara and Sugimura, 2012; Guirao et al., 2015), neglecting non-planarity (Hannezo et al., 2014; Bielmeier et al., 2016), curved cell edges (Brodland et al., 2014; Ishimoto and Morishita, 2014) 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, for example by ensuring no net change in internal energy under a homogeneous deformation (Appendix A). (This model is in the spirit of, but differs from, that of Okuda et al. Okuda et al. (2015), who proposed a drag force depending on an average of nearby vertex velocities.) The linearized stress/strain relationship (25, 26, 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 Ishihara et al. (2017a); simulations of large-amplitude plastic tissue deformations under external loading using a discrete cell model are provided by Pathmanathan et al. (2009).

Under the present vertex model, the perturbation stress of a pre-stressed 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 , which describes how anisotropic tissues resist deformation through reversible elastic deformations. This formulation extends previous approaches to upscaling the vertex-based model in spatially periodic (Murisic et al., 2015) and disordered isotropic networks (Nestor-Bergmann et al., 2017a). Exact expressions for elastic moduli for a given monolayer realisation are provided as explicit sums (for isotropic monolayers, these are (32) and (36)). Further work is required to derive a priori predictions for the behaviour of these quantities over multiple monolayer realisations. However a crude simplification of the estimate of the elastic shear modulus (37) 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. Merzouki et al. (2016), 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 non-monotonic behaviour close to the region II/III boundary (revealed most clearly in the stress ratios plotted in Figure 7). 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 (Figure 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 of epithelial tissues in various organ systems.

In summary, this study demonstrates how loading organises 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.

ANB was supported by a BBSRC studentship. OEJ acknowledges EPSRC grant EP/K037145/1. SW was supported by a Wellcome Trust/Royal Society Sir Henry Dale Fellowship [098390/Z/12/Z].

Appendix A A homogeneous small-amplitude strain

We derive the mappings of key geometric quantities under a small-amplitude homogenous and symmetric strain, , which transforms position vectors as , and then apply these mappings to the mechanical energy .

We assume that all quantities, , follow a mapping of the form , defined relative to the same cell and therefore temporarily drop the subscript. Tangents are defined as , giving , and hence . The length of an edge is given by . To linear order,


demonstrating that . The cell perimeter is , so that


Thus . It follows that . Similarly, . From the product rule, . Now, and using (38) gives . Thus, . Note that , ensuring that unit vectors are rotated but not stretched.

We now consider the deformation , writing


to linear order. Thus we see that , where in component form is