Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Advertisement

  • Loading metrics

Open Access

Peer-reviewed

Research Article

Linear viscoelastic properties of the vertex model for epithelial tissues

  • Sijie Tong,

    Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing

    Affiliation Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, New Jersey, United States of America

    ⨯
  • Navreeta K. Singh,

    Roles Data curation, Formal analysis, Investigation, Visualization, Writing – review & editing

    Affiliation Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, New Jersey, United States of America

    ⨯
  • Rastko Sknepnek ,

    Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Writing – original draft, Writing – review & editing

    r.sknepnek@dundee.ac.uk (RS); andrej@princeton.edu (AK)

    Affiliations School of Science and Engineering, University of Dundee, Dundee, United Kingdom, School of Life Sciences, University of Dundee, Dundee, United Kingdom

    ⨯
  • Andrej Košmrlj

    Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing

    r.sknepnek@dundee.ac.uk (RS); andrej@princeton.edu (AK)

    Affiliations Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, New Jersey, United States of America, Princeton Institute of Materials, Princeton University, Princeton, New Jersey, United States of America

    ⨯
  • Article
  • Authors
  • Metrics
  • Comments
  • Media Coverage
  • Peer Review

Abstract

Epithelial tissues act as barriers and, therefore, must repair themselves, respond to environmental changes and grow without compromising their integrity. Consequently, they exhibit complex viscoelastic rheological behavior where constituent cells actively tune their mechanical properties to change the overall response of the tissue, e.g., from solid-like to fluid-like. Mesoscopic mechanical properties of epithelia are commonly modeled with the vertex model. While previous studies have predominantly focused on the rheological properties of the vertex model at long time scales, we systematically studied the full dynamic range by applying small oscillatory shear and bulk deformations in both solid-like and fluid-like phases for regular hexagonal and disordered cell configurations. We found that the shear and bulk responses in the fluid and solid phases can be described by standard spring-dashpot viscoelastic models. Furthermore, the solid-fluid transition can be tuned by applying pre-deformation to the system. Our study provides insights into the mechanisms by which epithelia can regulate their rich rheological behavior.

Author summary

Epithelial tissues line organs and cavities in the body, and serve as barriers that separate organisms from their environment. Epithelia are robust yet adaptable; they have the ability to change their own viscoelastic behavior in response to internal or external stimuli by actively tuning the mechanical properties of the constituent cells and interactions between them. The mesoscopic mechanics of epithelia are commonly described with the vertex model. Here we present a detailed study of the linear rheological properties of the vertex model for both regular hexagonal and disordered cell configurations over a wide range of driving frequencies. The linear viscoelastic responses of the vertex model are mapped to standard spring-dashpot models. Our work, therefore, shows that the vertex model is a suitable base model to study the rich rheological behavior of epithelial tissues.

Citation: Tong S, Singh NK, Sknepnek R, Košmrlj A (2022) Linear viscoelastic properties of the vertex model for epithelial tissues. PLoS Comput Biol 18(5): e1010135. https://doi.org/10.1371/journal.pcbi.1010135

Editor: Tobias Bollenbach, Universitat zu Koln, GERMANY

Received: January 27, 2022; Accepted: April 25, 2022; Published: May 19, 2022

Copyright: © 2022 Tong et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: The simulation and analysis codes and initial configurations are available on GitHub: https://github.com/sknepneklab/RheoVM.

Funding: This research was primarily supported by NSF through the Princeton University’s Materials Research Science and Engineering Center DMR-2011750 (to AK) and by the Project X Innovation Research Grant from the Princeton School of Engineering and Applied Science awarded to AK. This research was also supported by the UK Biotechnology and Biological Sciences Research Council (Award BB/N009789/1) awarded to RS. AK and RS initiated this project during the KITP program "Symmetry, Thermodynamics and Topology in Active Matter" (ACTIVE20), which was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

The development and maintenance of tissues requires close coordination of mechanical and biochemical signaling [1–3]. There is, for instance, mounting evidence for the key role played by tissue material properties and their regulation during embryonic development [4]. Tissues must be able to adjust their mechanical properties in response to internal and external stimuli. In particular, epithelial tissues, which line all cavities in the body and demarcate organs, must sustain substantial mechanical stresses while also supporting numerous biological processes such as selective diffusion and absorption/secretion of molecules [5]. In homeostasis, epithelia must maintain their shape and resist deformation while remaining flexible. The tissue must also be able to regenerate and repair itself, often with fast turnover, e.g., in gut epithelia [6]. Furthermore, in morphogenesis, the epithelial tissue must take up a specific shape and function [7], but this shape is lost during metastasis when cancer cells invade surrounding healthy tissues [8]. All of these processes require that cells be able to move, often over distances much larger than the cell size. During cell migration, however, the epithelial tissue must maintain its integrity. It is, therefore, not surprising that epithelia exhibit rich viscoelastic behavior [9]. Unlike passive viscoelastic materials, an epithelial tissue can actively tune its rheological response, making the study of its rheology not only important for understanding biological functions but also an interesting problem from the perspective of the physics of active matter systems [10].

Collective cell migration has been extensively studied in biology [11] and biophysics [12]. In vitro studies of confluent cell monolayers [13–17] focused on the physical aspects of force generation and transmission and showed that cell migration is an inherently collective phenomenon. Some aspects of collective cell migration are remarkably similar to the slow dynamics of structural glasses [18–24]. This suggests that many of the observed behaviors share common underlying mechanisms and can be understood, at least at mesoscales (i.e., distances beyond several cell diameters), using physics of dense active systems [25]. A particularly intriguing observation is that tuning cell density [18, 26, 27], strength of cell-cell and cell-substrate interactions [28], or cell shape parameters [21, 29] can stop collective migration. In other words, the epithelium undergoes a fluid to solid transition. Signatures of such behavior have been reported in several in vitro [19, 30] and developmental systems [31–33]. This suggests that important aspects of morphogenetic development may rely on epithelial tissue’s ability to undergo phase transitions [4].

How an epithelial tissue responds to external and internal mechanical stresses depends on its rheological (i.e., material) properties. While there have been numerous studies focusing on the rheology of a single cell [34–36], much less is known about tissue rheology, particularly during development. In order to develop a comprehensive understanding of epithelial tissue mechanics, such insight is key. Though single cell measurements are valuable, the mechanics of an epithelial tissue can be drastically different from that of its constituent cells. The stiffness of cell monolayers, for example, is orders of magnitude higher than the stiffness of constituent cells, while the time dependent mechanical behaviors of monolayers in response to deformation vary depending on the magnitude of loading [37]. Embryonic cell aggregates have been shown to behave elastically (i.e., solid-like) at short timescales, but they flow like fluids at long timescales, which facilitates both the robustness needed to maintain integrity and the flexibility to morph during development [9]. Experiments have characterized the mechanical behaviors of epithelial tissues at various loading conditions, which led to a phenomenological description that models the relaxation properties of epithelial monolayers based on fractional calculus [38]. Notably, a recent particle-based model that includes cell division and apoptosis provided a plausible microscopic model for nonlinear rheological response [39]. Particle-based models are, however, unable to capture geometric aspects such as cell shape. It is, therefore, necessary to investigate rheological response in geometric models.

The vertex model [40–42] and more recent, closely related Voronoi models [21, 43, 44] have played an important role in modeling mechanics of epithelial tissues since they account for the shapes of individual cells and provide a link to cellular processes, such as cell-cell adhesion, cell motility, and mitosis [42]. These geometric models are also able to capture the solid to fluid transition and demonstrate rich and unusual nonlinear mechanical behavior [45, 46]. While the mechanical properties of the vertex and Voronoi models have been extensively studied, most works to date focused on the long-time behavior. These include studies of the quasistatic shear modulus [20], effective diffusion constant of cells related to the tissue viscosity [21], correlations between a structural property called “softness” and the likelihood of topological rearrangements of cells [47], and steady state flow profiles around a sphere dragged through the tissue [48]. The rheological properties of the vertex model that cover a broad range of timescales, however, have not yet been systematically explored. In this paper, we model the response of a model epithelial tissue adhered to a substrate by studying the response of the regular hexagonal and disordered cell configurations in the vertex model (see Fig 1) to applied oscillatory shear and bulk deformations of small amplitude, i.e., in the linear response regime. We measured the response stresses and used them to compute the storage and loss moduli in both the solid and fluid phases. We show that the dynamical response of the vertex model can be fitted to standard spring-dashpot viscoelastic models over seven decades in the driving frequency and that the solid-fluid transition can be tuned by applying pre-deformation to the system. Thus we argue that the vertex model makes a suitable basis for studies of dynamics of epithelial tissues beyond the quasistatic limit.

thumbnail
Download:
Fig 1. Epithelial tissue is represented as a polygon tiling of the plane subject to periodic boundary conditions.

We studied the rheology of the vertex model for both (a) regular hexagonal and (b) disordered tilings. Colors represent the number of neighbors of each cell; 4-white, 5-red, 6-gray, 7-blue, 8-yellow.

https://doi.org/10.1371/journal.pcbi.1010135.g001

Model and methods

Vertex model

In the vertex model, the state of an epithelial tissue is approximated as a polygonal tiling of the plane (Fig 1). The degrees of freedom are vertices, i.e., meeting points of three or more cell-cell junctions. In the simplest formulation, junctions are assumed to be straight lines. The energy of the vertex model is a quadratic function of cell areas and perimeters [41], i.e., (1) where KC and ΓC are the area and perimeter elastic moduli, and AC and AC0 are the actual and preferred areas of cell C, respectively. Similarly, PC and PC0 are the actual and preferred perimeters of the same cell. In this work, we assumed KC, ΓC, AC0, and PC0 to be identical for all cells (i.e., KC ≡ K, ΓC ≡ Γ, AC0 ≡ A0, PC0 ≡ P0). Further, we fixed the values of K and A0, and measured the energy in units of , stresses in units of KA0, and lengths in units of . Since the ratio between the perimeter and area elastic moduli does not qualitatively change the behavior of the vertex model [20, 41], we fixed that ratio to Γ/(KA0) ≈ 0.289 for all simulations. The only variable parameter in simulations was the preferred cell perimeter P0, which sets the dimensionless cell-shape parameter, defined as the ratio .

The cell-shape parameter, p0, plays a central role in determining whether the system behaves as a fluid or solid [20]. Bi, et al. [20] argued that the rigidity transition occurs at p0 = pc ≈ 3.812 for a disordered polygonal tiling, while Merkel, et al. [49] reported p0 = pc ≈ 3.92. For a regular hexagonal tiling, the transition point is at [50]. In the fluid phase, the energy barrier for neighbor exchanges vanishes and cells can flow past each other [51]. As p0 is reduced below pc, the energy barrier becomes finite, neighbor exchanges cease and the system becomes solid. While the transition point for regular hexagonal tilings can be understood in terms of the mechanical stability and the excess perimeter [45, 52], the mechanism that leads to a larger value for random tilings is more subtle and only partly understood [53]. For example, recent studies [49, 54] have shown that the rigidity transition of random tilings depends on the procedure used to generate the tilings. The presence of vertices with coordination greater or equal to four, as well as the presence of cells with five or less neighbors, increases the critical value of the cell-shape parameter p0 [54].

Simulation setup

We first studied the rheology of regular hexagonal tilings (Fig 1a) subject to periodic boundary conditions. The shape of the simulation box was chosen to be as close to a square as allowed by the geometry of a hexagon, and the area of the box was such that it accommodated N cells of area AC that matched the preferred areas A0. Most simulations started with hexagonal tiling with Nx = 15 cells in the horizontal direction (i.e., N = 240 cells in total, Fig 1a). Simulations of larger system sizes (Nx = 37, 51, i.e., N = 1406, 2652 total cells, respectively) were performed for a subset of values of p0 to explore the finite size effects. No quantitative differences between the system with N = 240 cells and larger systems were observed. All simulations were performed with an in-house code [55] and snapshots of cell configurations were visualized with ParaView [56].

For the solid phase with p0 ≲ 3.722, the ground state of the energy in Eq (1) is the honeycomb lattice [50], and it was directly used to investigate rheological properties. Note that there was some residual hydrostatic stress due to the mismatch of actual cell perimeters PC from their preferred values P0, which could be eliminated by appropriate rescaling of the simulation box. This hydrostatic stress, however, does not qualitatively affect the rheological behavior of the system (see Sec D in S1 File for further discussion). For the fluid phase with p0 ≳ 3.722, the hexagonal tiling corresponds to a saddle point of the energy in Eq (1) [50]. A small random perturbation was applied to each vertex, i.e., each vertex was displaced from its original position in the hexagonal tiling by a vector δri = δxiex + δyiey, where δxi and δyi were Gaussian random variables with zero mean and standard deviation ; the system was then relaxed using the FIRE algorithm [57] to reach a local energy minimum with the relative accuracy of 10−12. Note that the energy landscape in the fluid phase has many local minima and a large number of soft modes (see Sec H in S1 File). We repeated simulations to investigate rheological properties for multiple configurations corresponding to different local energy minima.

The study of the rheological properties of the vertex model for a regular hexagonal tiling is appealing since one can make comparisons to analytical treatments. The regular hexagonal tiling is, however, a rather crude approximation of real epithelial tissues, which are typically irregular [58]. Therefore, we also investigated the rheology of the vertex model of disordered tilings (Fig 1b). The disordered tilings were created as follows (see Fig A in the Sec A in S1 File for a schematic illustration). We used the random sequential addition algorithm [59] to place N = 200 seed points inside a square box of size L = 15 without overlaps. We then created periodic images of the seed points and used SciPy to build the periodic Voronoi tessellation [60]. The preferred area of each cell was set to A0 = L2/N. The energy of the system given in Eq (1) was then relaxed using the FIRE algorithm to reach a local minimum. During the energy minimization, T1 transitions (exchanges of cell neighbors) were allowed but were not common. We generated an ensemble of 10 different random initial configurations using different values of the random number generator seed and repeated rheology simulations for each of those configurations to probe the rheology for a range of values of p0.

Dynamics and probing the rheology

In order to probe the dynamic response of the vertex model, we need to specify the microscopic equations of motion for vertices. Assuming the low Reynolds number limit, which is applicable to most cellular systems due to their slow speed, inertial effects can be neglected [61]. The equations of motion are then a force balance between friction with the substrate and elastic forces Fi due to deformations of cell shapes, i.e., (2) Here, we assume that friction between the tissue and the substrate arises from binding and unbinding of adhesion molecules. In particular, on time scales much longer than the characteristic unbinding time, the tissue–substrate adhesive bonds undergo stick-and-slip processes leading to a form of viscous friction [62–64]. In the above Eq (2), ri is the position vector of vertex i in a laboratory frame of reference, is the mechanical force on vertex i due to deformation of cells surrounding it, γ is the friction coefficient, and dot denotes the time derivative. Therefore, each vertex experiences dissipative drag proportional to its instantaneous velocity. We fixed the value of γ in simulations, which sets the unit of time as γ/(KA0). Furthermore, we neglected thermal fluctuations and hence omit the stochastic term in Eq (2). This is a reasonable assumption since typical energy scales in tissues significantly exceed the thermal energy, kBT, at room temperature T, where kB is the Boltzmann constant. It is, however, worth noting that there are other sources of stochasticity in epithelia (e.g., fluctuations of the number of force-generating molecular motors) which are important for tissue scale behaviors [65]. Here, we did not consider such effects but note that they could be directly included in the model as additional forces in Eq (2).

We applied an oscillatory affine deformation to investigate the rheological behavior of the vertex model. The affine deformation can be described by a deformation gradient tensor defined as , where the mapping x = x(X0, t) maps the reference configuration X0 to a spatial configuration x at time t. Note that the affine deformation was applied to all vertices inside the bounding box as well as to the periodic images of vertices, which ensured deformation of cells across periodic boundaries. Specifically, for a vertex i with a position vector ri(t) in a deformed configuration, we also consider its periodic images at positions ri(t) + m1a1(t) + m2a2(t), where a1,2(t) are the unit cell vectors of the periodic box and m1,2 are integers. For cells that cross the periodic boundary, junction lengths are computed using the minimum image convention. Note that the unit cell vectors of the periodic box transform as due to the affine deformation.

The deformation gradient of simple shear is and of biaxial deformation is , where ϵ(t) = ϵ0 sin(ω0t) and ω0 is the frequency of the oscillatory deformation. In all simulations, we used a small magnitude of deformation, i.e., ϵ0 = 10−7, so that we probed the linear response and the measured moduli were independent of the magnitude of the deformation. In every time step after the affine deformation was applied, the system evolved according to the overdamped dynamics in Eq (2). During the oscillatory deformations, T1 transitions were allowed but were not common. Equations of motion were integrated using the first-order Euler method [66] with the time step Δt ≈ 0.00866γ/(KA0) when the frequency of oscillatory deformation ω0γ/(KA0) < 29.02, but with a smaller time step Δt ≈ 0.000866γ/(KA0) when ω0γ/(KA0) > 29.02 so that there were enough sampling points (at least 25) over one period of oscillatory deformation.

The response stress tensor, , for each cell C was computed using the formalism introduced in Refs. [67–69] as (3) where the summation is over all junctions e belonging to cell C. Here, is the hydrostatic pressure inside a cell, is the unit tensor, and is the tension along the junction e with le being a vector joining the two vertices on it [67–69]. The average stress tensor , with wC = AC/∑C AC, was used as a measure for the response of the system. Measurements of the response stresses for each cell [see Eq (3)] and the entire system were taken 25 times within each cycle of oscillatory deformation.

To ensure that we were probing the steady state, we performed the following analysis. For example, in the case of shear deformation, the shear stress signal was divided into blocks of length T = 3T0, each containing 3 cycles of the time period T0 = 2π/ω0 of the driving shear deformation. Within each block n, we performed the Fourier transform of τ(t) and obtained as (4) where n is a positive integer. Similar Fourier transform analysis was performed for the strain, ϵ(t), of which the Fourier transform is denoted as . The length of the simulation was chosen such that it contained a sufficient number of blocks in order for the to reach a steady state value . The obtained steady state value of was used to calculate the dynamic shear modulus at a given frequency ω0 of applied shear strain. We ensured that simulations ran long enough to reach a steady state. An analogous procedure was applied to the hydrostatic stress, , in the case of the bulk deformation. Please refer to Sec C in S1 File for a representative example of the steady state analysis.

Results

Response to a shear deformation

The hexagonal ground state in the solid phase and states corresponding to local energy minima in the fluid phase were used to investigate the rheological behavior by applying an oscillatory affine shear deformation to the substrate (Fig 2a and 2b). Due to the binding and unbinding of adhesive molecules, deformation of the tissue follows the deformation of the substrate on short timescales, and then tissue can relax on longer timescales. Thus, at each time step, we first applied the affine shear deformation to the simulation box and all vertices, which was followed by internal relaxation of the vertices according to Eq (2). The affine simple shear deformation can be described by a deformation gradient tensor, , where ϵ(t) = ϵ0 sin(ω0 t). Sufficiently small amplitude ϵ0 = 10−7 ≪ 1 was used to probe the linear response properties.

thumbnail
Download:
Fig 2. Storage and loss shear moduli in the solid (top row) and fluid phase (bottom row) for hexagonal tilings.

(a-b) An overlay of the representative reference (grey) and sheared (yellow) configurations in (a) the solid and (b) the fluid phase. The magnitude of the shear is highly exaggerated for demonstration purposes. (c-d) Representative storage (G′) and loss (G″) shear moduli as functions of the shearing frequency, ω0, for different values of the cell-shape parameter, p0. Dashed curves are the fits based on (c) the Standard Linear Solid (SLS) model in the solid phase [see Eq (5)] and (d) the Burgers model in the fluid phase [see Eq (8)]. (e-f) The collapse of the moduli curves for different values of p0 for (e) the solid phase and (f) the fluid phase. The insets show the representation of (e) the SLS model and (f) the Burgers model in terms of the springs and dashpots. The majority of the data corresponds to the system of nearly square shape with Nx = 15 cells in the horizontal direction, and we also show examples of larger systems with Nx = 37 and Nx = 51 cells in the horizontal direction.

https://doi.org/10.1371/journal.pcbi.1010135.g002

We measured the response stresses as described in the Model and methods section above. The dynamic shear modulus was then calculated at a given frequency ω0 of applied shear strain, where and are the Fourier transforms of the response shear stress and the applied strain ϵ(t), respectively (see Model and methods). We ensured that simulations ran long enough to reach a steady state (see Model and methods and Sec C in S1 File). The real part of the dynamic shear modulus, G′ = Re(G*), is the storage shear modulus and the imaginary part, G″ = Im(G*), is the loss shear modulus. The storage shear modulus corresponds to the in-phase response and measures the elastic (i.e., reversible) response of the system, while the loss shear modulus corresponds to the out-of-phase response and measures the system’s irreversible dissipation [70] (see also Sec B in S1 File). For systems under an oscillatory simple shear, storage and loss shear moduli were obtained for different values of p0 and different system sizes in the solid and the fluid phases for a broad range of driving frequencies ω0 spanning over seven orders of magnitude, as shown in Fig 2c and 2d. Most simulations were performed for systems with nearly square shapes with Nx = 15 cells in the horizontal direction. We repeated several simulations for systems with Nx = 37 and Nx = 51, which showed that the finite size effects are negligible (Fig 2c–2f).

In the solid phase there are two different regimes (see Fig 2c). At low frequencies, ω0, the storage shear modulus G′ has a constant value, while the loss shear modulus scales as G″ ∝ ω0. At high frequencies, the storage shear modulus G′ has a higher constant value, while the loss shear modulus scales as . Such rheological behavior is characteristic of the Standard Linear Solid (SLS) model [70]. Storage and loss shear moduli for the SLS model are [70], respectively, where we used the representation of the SLS model (Fig 2e, inset) that consists of a spring with elastic constant E2 connected in parallel with a Maxwell element, which comprises a spring with elastic constant E1 and a dashpot with viscosity η1 connected in series. The above expressions in Eq (5) were used to fit the storage and loss shear moduli obtained from simulations. The fitted curves, represented with dashed lines in Fig 2c, show an excellent match with the simulation data, indicating that the SLS model is indeed appropriate to describe the shear rheology in the solid phase. This was also confirmed in Fig 2e, where we collapsed the storage and loss shear moduli for different values of the shape parameter, p0, by rescaling the moduli and frequencies with the fitted values of spring and dashpot constants. Note that the SLS response in the solid phase is consistent with recent experiments on suspended MDCK monolayers [71].

As the value of the p0 increases, we observe that the storage shear modulus reduces at all frequencies and that the loss shear modulus reduces at high frequencies. Furthermore the crossover between the two regimes shifts towards lower frequencies (Fig 2c). This is because the elastic constants E1 and E2 decrease linearly with increasing p0 and they become zero exactly at the solid-fluid transition with p0 = pc ≈ 3.722 (Fig 3a). The dashpot constant η1 is nearly independent of p0 (Fig 3b) and scales with the friction parameter γ, which is the only source of dissipation in the vertex model. The crossover between the two regimes for both the storage and loss shear moduli corresponds to a characteristic timescale, η1/E1, which diverges as ∼ γ(KA0)−1(pc − p0)−1 as p0 approaches the solid-fluid transition (Fig 3c) due to the vanishing elastic constant (Fig 3a). Note that the values of the elastic constants E1 and E2 can be estimated analytically. In the quasistatic limit (ω0 → 0), the external driving is sufficiently slow that the system can relax internally. In this limit, Murisic, et al. [72] showed that the storage shear modulus is (6) where α(p0, Γ/KA0) is a scaling factor chosen such that the hydrostatic stress vanishes once the system box size is rescaled from L to αL (see Sec D in S1 File). In the high frequency limit (ω0 → ∞), on the other hand, the system follows the externally imposed affine deformation and has no time for internal relaxation. Thus, by considering the energy cost for a hexagonal tiling under affine deformation, we obtained the storage shear modulus (see Sec G in S1 File) (7) The above Eqs (6) and (7) were used to extract the values of elastic constants E1 and E2, which showed excellent agreement with the fitted values from simulations (Fig 3a).

thumbnail
Download:
Fig 3.

(a-b) Fitted values of spring-dashpot models for hexagonal tilings under simple shear. (a) Elastic constants as a function of target cell-shape parameter, p0. In the solid phase (i.e., for p0 < pc ≈ 3.722), fitted values of the spring constants show excellent match with the analytical predictions obtained from Eqs (6) and (7) (dashed lines). Inset shows the spring constants near the critical point. (b) Dashpot viscosity constants as a function of the target cell-shape parameter, p0. (c-d) Characteristic timescales in (c) the solid and (d) fluid phase for hexagonal tilings obtained from the fitted values of the elastic constant and the dashpot viscosity. The normalization factor t* = γ/(KA0) sets the unit of time. For the fluid phase (i.e., for p0 > pc ≈ 3.722), errorbars correspond to the standard deviation for simulations that were repeated for configurations that correspond to different local energy minima.

https://doi.org/10.1371/journal.pcbi.1010135.g003

In the fluid phase, the storage and loss shear moduli show a markedly different behavior (Fig 2d). There are three different regimes with two crossover frequencies, which correspond to two characteristic timescales. At low frequencies, ω0, the storage shear modulus and the loss shear modulus G″(ω0) ∝ ω0. The storage modulus approaches 0 for ω0 → 0, which indicates that the system is indeed a fluid. At high frequencies the storage shear modulus has a constant value, while the loss shear modulus scales as . To capture this behavior we used the Burgers model, which consists of two Maxwell models connected in parallel (Fig 2f, inset), to fit the shear moduli measured in the simulations. The storage and loss shear moduli for a Burgers model are [70], respectively, where p1 = η1/E1 + η2/E2, p2 = η1η2/(E1 E2), q1 = η1 + η2, q2 = η1η2(E1 + E2)/(E1E2). The dashed curves in Fig 2d show fits of the storage and loss shear moduli for a range of values of p0, which show good agreement with simulations. Unlike for the solid phase, it is not possible to collapse the data for storage and loss shear moduli onto single universal curves because the fluid phase is characterized by two independent timescales η1/E1 and η2/E2. Thus we show two different collapses for the storage and loss shear moduli in the low frequency range (Fig 2f) and in the high frequency range (Fig E in the Sec E in S1 File).

As the value of the p0 decreases, we observe that both the storage and loss shear moduli reduce at intermediate and high frequencies, but they increase at low frequencies (Fig 2d). We also observe that the first crossover shifts towards lower frequencies, while the second crossover remains at approximately the same frequency. This is because the elastic constants E1 and E2 decrease linearly toward zero as p0 approaches the solid-fluid transition at pc ≈ 3.722 (Fig 3a). The dashpot constant η2 also decreases linearly toward zero, while the dashpot constant η1 increases but remains finite as p0 approaches the solid-fluid transition (Fig 3b). As a consequence, one of the characteristic timescales η1/E1 ∼ γ(KA0)−1(p0 − pc)−1 diverges, while the second timescale η2/E2 ∼ γ(KA0)−1 remains finite as p0 approaches the solid-fluid transition (Fig 3d). The diverging characteristic timescale captures the macroscopic behavior of the system, while the second timescale captures the microscopic details of the vertex model. Note that at the solid-fluid transition there is a discontinuous jump in the values of the dashpot constant η1 (see Fig 3b). This is because at p0 = pc the storage and loss shear moduli are identically equal to zero (G′(ω0) = G″(ω0) ≡ 0) due to the vanishing elastic constants (E1 = E2 = 0), while the dashpot constants can have arbitrary values [see Eqs (5) and (8)]. Finally, we note that the values of the spring and dashpot constants are somewhat sensitive to the local energy minimum configuration used to probe the response in the fluid phase. The errorbars in Fig 3 show standard deviation for different configurations that were obtained by using the same magnitude of the initial perturbation (see Model and methods). In Fig F in the Sec F in S1 File, we show how the values of the spring and dashpot constants are affected when configurations were obtained by using different magnitudes of the initial perturbation.

Cells in real epithelial tissues are, however, unlikely to have a perfect hexagonal shape and form a honeycomb tiling. The observed tilings are disordered, often with a rather specific distribution of the number of neighbor cells conserved across several species [73]. To mimic the geometry of real tissues, we constructed an ensemble of uncorrelated disordered tilings of polygons corresponding to local energy minima at different values of p0 (see Fig 1b and Model and methods). We then probed the shear rheology of each such configuration following the same procedure as for the hexagonal tilings. Fig 4 shows the average storage and shear moduli. We found that the critical value of p0 for the solid-fluid transition was at pc ≈ 3.93, which is consistent with refs. [49, 54], but somewhat higher than what was reported in [20].

thumbnail
Download:
Fig 4. Average storage and loss shear moduli in the solid and fluid phase for disordered tilings.

(a,c) Average storage (G′) and loss (G″) shear moduli as functions of the shearing frequency, ω0, for different values of the cell-shape parameter, p0, (a) deep in the solid phase and (c) deep in the fluid phase. The error bars represent the standard error of the mean. (b,d) The collapse of the moduli curves for different values of p0 for (b) the solid phase and (d) the fluid phase. The insets show the representation of (b) the Standard Linear Solid (SLS) model and (d) the Burgers model in terms of the springs and dashpots. (e,f) Average storage (G′) and loss (G″) shear moduli as functions of the shearing frequency, ω0, for intermediate values of the cell-shape parameter, p0, in (e) the solid phase and (f) the fluid phase. Dashed curves are the fits based on (a,e) the SLS model in the solid phase [see Eq (5)] and (c,f) the Burgers model in the fluid phase [see Eq (8)]. A representative example of random cell configurations used to produce these plots is shown in Fig 1b.

https://doi.org/10.1371/journal.pcbi.1010135.g004

For p0 < pc values corresponding to the system being deep in the solid phase (Fig 4a), the storage and loss moduli are described accurately by the Standard Linear Solid (SLS) model, the same as for the hexagonal tiling in the solid phase (Fig 2c). As p0 increases, however, a second shoulder develops in the loss moduli (see p0 = 3.71 and p0 = 3.77 curves in Fig 4a), which indicates the presence of multiple time scales. The fits to the SLS model shown with the dashed lines also begin to deviate from the measured moduli. The scaling collapse is only possible for values of p0 deep in the solid phase (Fig 4b). This supports the observation that SLS is no longer able to capture the rheology in the solid phase as p0 approaches the critical point.

In the opposite limit, i.e., when the value of p0 > pc is deep in the fluid phase (Fig 4c), the storage and loss moduli can be modeled with the Burgers model, the same as for the local energy minima states relaxed from the hexagonal tiling in the fluid phase (Fig 2d). The fits represented by the dashed lines, however, deviate from the measured moduli as p0 decreases (see p0 = 3.97 and p0 = 3.99 curves in Fig 4c). Fig 4d shows the rescaling of the moduli and frequencies by the fitted spring and dashpot constants.

Near the critical value (i.e., for p0 = 3.93 and p0 = 3.95), the ensemble of random tilings contains the solid and the fluid configurations (see Fig I in the Sec I in S1 File), which was determined based on the presence or absence of non-trivial zero modes. We separated the solid and the fluid configurations and calculated average storage and loss shear moduli on each set. Fig 4e and 4f show the average storage and loss moduli for values of p0 close to the critical value in the solid phase (Fig 4e) and the fluid phase (Fig 4f). The dashed curves are the fits based on the SLS model in the solid phase and the Burgers model in the fluid phase, which do not fully capture the behavior of the measured moduli curves due to the presence of multiple time scales. As the value of p0 approaches the critical value, the spread of the moduli increases, especially for low frequencies, which is captured by the size of error bars. This can also be seen in Fig I in the Sec I in S1 File, which presents the raw data of storage and loss shear moduli at different values of p0.

In Fig 5, we summarize the fitted values of spring-dashpot models. The values of spring constants decrease as the system approaches the solid-fluid transition, while the values of dashpot constants diverge near the transition. In the intermediate regime (shaded regions in Fig 5), the SLS model (in the solid phase) and the Burgers model (in the fluid phase) cannot accurately fit the measured moduli due to the presence of additional timescales.

thumbnail
Download:
Fig 5. Fitted values of spring-dashpot models for disordered tilings under simple shear.

(a) Elastic constants as a function of the target cell-shape parameter, p0. (b) Dashpot viscosity constants as a function of the target cell-shape parameter, p0. The shaded regions indicate the intermediate regime between the solid and fluid phases.

https://doi.org/10.1371/journal.pcbi.1010135.g005

Response to bulk deformations

We further studied the bulk rheological properties of the hexagonal tilings by applying an oscillatory biaxial deformation to the substrate (Fig 6a and 6b) described by the deformation gradient , where ϵ(t) = ϵ0 sin(ω0t). We applied a sufficiently small amplitude ϵ0 = 10−7 ≪ 1 to probe the linear response properties characterized by the average normal stress . As in the simple shear test, we then computed the dynamic bulk modulus as from which we obtained the storage bulk modulus B′ = Re(B*) and the loss bulk modulus B″ = Im(B*) (see Fig 6c and 6d).

thumbnail
Download:
Fig 6. Loss and storage bulk moduli in the solid (top row) and fluid phase (bottom row) for hexagonal tilings.

(a-b) An overlay of the representative reference (grey) and biaxially deformed (yellow) configurations in (a) the solid and (b) the fluid phase. The magnitude of the bulk deformation is highly exaggerated for demonstration purposes. (c-d) Representative storage (B′) and loss (B″) bulk moduli as functions of the deformation frequency, ω0, for different values of the cell-shape parameter, p0. For the solid phase in (c), the loss bulk modulus B″ ≡ 0. For the fluid phase in (d), dashed curves are the fits based on the Standard Linear Solid (SLS) model [see Eq (5)]. (e-f) The collapse of the moduli curves for different values of p0 for (e) the solid phase and (f) the fluid phase. The insets show the representation of (e) the spring model and (f) the SLS model in terms of the springs and dashpots. In panel (e), Btheory corresponds to the analytical prediction in Eq (9) for the storage bulk modulus in the solid phase.

https://doi.org/10.1371/journal.pcbi.1010135.g006

In the solid phase, the storage bulk modulus is independent of the driving frequency and the loss bulk modulus is zero. This is because in the solid phase, the hexagonal tiling is stable to biaxial deformation and there is no relative motion of vertices with respect to the substrate, which is the sole source of dissipation. Thus the response of the system can be captured by a single spring Esolid (Fig 6e, inset). The measured value of the storage bulk modulus matches the analytical prediction, (9) by Staple, et al. [50], where the hexagonal tiling is assumed to undergo affine deformation under biaxial deformation. Storage bulk moduli, normalized by Btheory, for different values of p0 all collapse to 1 (Fig 6e).

In the fluid phase, the bulk response behavior of the system can be described by the SLS model (Fig 6f, inset). While it might appear counter-intuitive to model a fluid with the SLS model, this is a direct consequence of the fact that in the fluid state, the bulk modulus is finite but the shear modulus vanishes, i.e., the fluid flows in response to shear but resists bulk deformation. The fitted storage and loss bulk moduli for the SLS model [see Eq (5)] show an excellent match with the simulation data (Fig 6d). This was also confirmed in Fig 6f, where we collapsed the storage and loss bulk moduli for different values of p0.

The fitted values of elastic spring and dashpot viscosity constants for different values of p0 are plotted in Fig 7. In the fluid phase, the storage bulk modulus in the high frequency limit B′(ω0 → ∞) = E1 + E2 [see Eq (5)] continuously increases from the value for the solid phase Btheory [see Eq (9)] as the system transitions from solid to fluid (Fig 7a). The storage bulk modulus in the quasistatic limit B′(ω0 → 0) = E2 [see Eq (5)] emerges at the transition point with a finite value and increases as p0 increases from pc (Fig 7a). Fig 7b shows that the dashpot constant η1 diverges as the p0 decreases toward pc. Thus, the characteristic timescale η1/E1 also diverges (Fig 7c), but for a different reason than for the shear deformation, where the spring constant E1 is vanishing (see Fig 3). Finally, we note that, unlike for the response to shear, the values of the spring and dashpot constants for bulk deformation are not sensitive to the local energy minimum configuration used to probe the response in the fluid phase, which is reflected by the very small errorbars in Fig 7. This is because the bulk moduli are dominated by the changes in cell areas.

thumbnail
Download:
Fig 7. Fitted values of spring-dashpot models for the system under bulk deformation as a function of the target cell-shape parameter, p0.

(a) Elastic constants as a function of the target cell-shape parameter, p0. In the solid phase (p0 < pc ≈ 3.722), the bulk storage modulus Esolid agrees with the analytical prediction Btheory in Eq (9) (dashed line). At the solid-fluid transition point (p0 = pc ≈ 3.722), it continuously changes to the high frequency limit of the bulk storage modulus, i.e., B′(ω0 → ∞) = E1 + E2, of the fluid phase. The low frequency limit of the bulk storage modulus is B′(ω0 → 0) = E2 in the fluid phase. (b) Dashpot viscosity constant as a function of the target cell-shape parameter, p0. (c) Characteristic timescales in the fluid phase obtained from the fitted values of the elastic constant and the dashpot viscosity. The normalization factor t* = γ/(KA0) sets the unit of time. For the fluid phase (p0 > pc ≈ 3.722), errorbars correspond to the standard deviation for simulations that were repeated for configurations that correspond to different local energy minima.

https://doi.org/10.1371/journal.pcbi.1010135.g007

The same procedures were applied to the ensemble of disordered tilings to probe the bulk rheology. Fig 8 shows the average storage and loss bulk moduli for different values of p0. When the system is deep in the solid phase (Fig 8a) or deep in the fluid phase (Fig 8c), the bulk rheology can be described by the SLS model, which is confirmed by the fits (dashed curves) and the collapse in Fig 8b and 8d. The fitted values of spring and dashpot constants for different values of p0 are shown in Fig 9. As p0 approaches the value of solid-fluid transition, the fits based on the SLS model deviate from the measured moduli curves (Fig 8a and 8c). At intermediate frequencies the storage moduli have a lower slope than predicted by the SLS model and the peak in the loss moduli is flattened and a second peaks starts to develop (p0 = 3.71, 3.77, and 3.80 in Fig 8a and p0 = 3.99 in Fig 8c). Fig 8e shows the storage and loss moduli for values of p0 near the solid-fluid transition, and the collapsed data is shown in Fig 8f. As the value of p0 approaches the critical value, the spread of the moduli increases, especially for low frequencies, which is seen in Fig J in the Sec J in S1 File, that presents the raw data of storage and loss bulk moduli at different values of p0.

thumbnail
Download:
Fig 8. Average storage and loss bulk moduli in the solid and fluid phase for disordered tilings.

(a,c) Average storage (B′) and loss (B″) bulk moduli as functions of the deformation frequency, ω0, for different values of the cell-shape parameter, p0, (a) deep in the solid phase and (c) deep in the fluid phase. The error bars represent the standard error of the mean. (b,d) The collapse of the moduli curves for different values of p0 for (b) the solid phase and (d) the fluid phase. The insets show the representation of the Standard Linear Solid (SLS) model in terms of the springs and dashpots. (e) Average storage (B′) and loss (B″) bulk moduli as functions of the deformation frequency, ω0, for intermediate values of the cell-shape parameter, p0. (f) The collapse of the moduli curves for for intermediate values of the cell-shape parameter, p0. Dashed curves in (a,c,e) are the fits based on the SLS model [see Eq (5)].

https://doi.org/10.1371/journal.pcbi.1010135.g008

thumbnail
Download:
Fig 9. Fitted values of spring-dashpot models for disordered tilings under bulk deformation.

(a) Elastic constants as a function of the target cell-shape parameter, p0. (b) Dashpot viscosity constant as a function of the target cell-shape parameter, p0. The shaded regions indicate the intermediate regime between the solid and fluid phases.

https://doi.org/10.1371/journal.pcbi.1010135.g009

Response to a shear deformation of a uniaxially pre-deformed system

The solid-fluid transition for the regular hexagonal tiling occurs when p0 ≈ 3.722, above which the hexagonal tiling is unstable. This is consistent with the vanishing of the affine shear modulus in Eq (7) at the transition point. If the regular hexagonal tiling is pre-compressed or pre-stretched uniaxially by a factor a, which is described by the deformation gradient , then the high frequency limit of the linear storage shear modulus that is dominated by affine deformation becomes (see Sec G in S1 File), (10) By setting the affine shear modulus to 0, we obtained the solid-fluid transition boundary in the a − p0 plane as (11) The above analytical prediction for the phase boundary (Fig 10a, blue line) shows an excellent agreement with the stability analysis in terms of the eigenvalues of the Hessian matrix of the energy function [74] (Fig 10a, red dots). A given configuration is stable if all eigenvalues of the Hessian matrix are positive and the loss of mechanical stability occurs when the lowest eigenvalue becomes 0. For a given p0, the value of the lowest eigenvalue reduces with decreasing a, i.e., as the magnitude of compression is increased. Thus, the compression (stretching) shifts the solid-fluid transition towards the lower (higher) values of p0 (see Fig 10a).

thumbnail
Download:
Fig 10. Tuning the solid to fluid transition by applying uniaxial pre-deformation.

(a) The solid-fluid transition boundary in the a − p0 plane, where a measures the amount of uniaxial pre-deformation described by the deformation gradient . Blue line shows the analytical prediction from Eq (11), which matches the stability analysis with the Hessian matrix (red dots). (b,c) The fitted values of the (b) spring and (c) dashpot constants for the SLS model in the solid phase [see Eq (5)] and the Burgers model in the fluid phase [see Eq (8)] when the system is under uniaxial compression (a = 0.95), no pre-deformation (a = 1.00), and under uniaxial tension (a = 1.05).

https://doi.org/10.1371/journal.pcbi.1010135.g010

We also probed the response to oscillatory shear applied to uniaxially pre-compressed and pre-stretched systems. This analysis was done on the uniaxially deformed hexagonal tiling in the solid phase as well as a system in the fluid phase obtained by relaxing the unstable, uniaxially deformed hexagonal tiling after an initial random perturbation (see Model and methods). The response to the shear deformation is qualitatively similar and can still be described by the SLS model in the solid phase and the Burgers model in the fluid phase. Fig 10b and 10c shows fitted values of the parameters for spring-dashpot models when the system is under uniaxial compression (a = 0.95), no pre-deformation (a = 1.00, i.e., same as Fig 3a and 3b), and uniaxial tension (a = 1.05). In both the solid and fluid phases, all spring elastic constants decrease to 0 as p0 approaches the critical value predicted by Eq (11). The dashpot constant η1 remains constant in the solid phase. Once the system enters the fluid phase as p0 increases, a new dashpot constant η2 emerges and increases from 0, while the value of the dashpot constant η1 decreases. As in the simple shear case, we note that the dashpot constant η1 has a discontinuous jump at the solid-fluid transition (see Fig 3c) and that the values of the spring and dashpot constants are somewhat sensitive to the local energy minimum configuration used to probe the response in the fluid phase. The errorbars in Fig 3 show standard deviation for configurations that were obtained by using different random initial perturbation (see Model and methods). Finally, we note that besides the uniaxial pre-deformation, the solid-fluid transition point can be tuned by other modes of pre-deformation (see Fig G and Sec G in S1 File).

Discussion and conclusions

We have performed a detailed analysis of the rheological properties of the vertex model subject to small-amplitude oscillatory deformations over seven orders of magnitude in the driving frequency. Our analysis shows that the vertex model exhibits non-trivial viscoelastic behavior that can be tuned by a single dimensionless geometric parameter—the shape parameter, p0. In order to characterize the response, we constructed constitutive rheological models that use combinations of linear springs and dashpots connected in series and in parallel. These models allowed us to match the shear response of the vertex model to that of the Standard Linear Solid model in the solid phase and the Burgers model in the fluid phase. In the low-frequency, i.e., quasistatic regime, our results are fully consistent with many previous studies [20, 21, 50, 72]. Our work, however, provides insights into the time-dependent response of the vertex model over a broad range of driving frequencies, which is important if one is to develop full understanding of the rheological properties of the vertex model and how they inform our understanding of epithelial tissue rheology.

While the SLS and the Burgers model accurately describe rheology of the vertex model of disordered tilings deep in the solid and liquid phases, respectively, these models deviate from the data for p0 values in the vicinity of the solid-fluid transition. This is because close to the transition points additional relevant time scales start to emerge. As shown in Fig K in the Sec K in S1 File, adding additional Maxwell elements in parallel to the spring-dashpot models increases accuracy of the fits. This is to be expected since each Maxwell element introduces a new time scale. The physical interpretation of these additional time scales has clearly to do with the local arrangements of the cells for a particular disordered configuration but tying it to a specific cell pattern is, however, not easy. In addition, we also found that for disordered tilings the loss shear modulus crosses over from the linear scaling in frequency at low ω0 to the ∼ ωα with α ≈ 0.73 at intermediate frequencies (Fig L in the Sec L in S1 File). The crossover moves to lower frequencies as the system size increases. This behavior suggests a large (potentially infinite) number of relevant timescales.

It is important to note that we considered only friction between cells and the substrate and neglected any internal dissipation within the tissue. Therefore, the dissipation is solely due to relative motion of the cells with respect to the substrate as a result of non-affine relaxation of the tissue. The approach used in this study, therefore, would not be suitable for modeling the rheological response of epithelia not supported by a solid substrate, e.g., for early stage embryos or suspended epithelia in the experiments of Harris, et al. [37]. Furthermore, dissipative processes in epithelia are far more complex than simple viscous friction and are not fully understood. It has, for example, recently been argued that internal viscoelastic remodelling of the cortex can lead to interesting collective tissue behaviors [75]. We have addressed some of these questions in a separate recent work [76] using a semi-analytic approach based on the normal-mode expansion.

In this work we kept the ratio Γ/(KA0) fixed, since this ratio does not qualitatively change the behavior of the vertex model [20, 41]. In Ref. [76], however, we further explored the effect of parameter Γ/(KA0), where we showed that the perimeter stiffness Γ and the area stiffness K affect the low-frequency and high-frequency rheological behavior, respectively.

We also showed that the critical value for the solid-fluid transition can be tuned by applying pre-deformation. Interestingly, under uniaxial and biaxial (i.e., isotropic) pre-compression the solid to fluid transition shifts to lower values of p0, leading to the non-intuitive prediction that one can fluidize the system by compressing it. This is, however, unsurprising, since the transition is driven by a geometric parameter that is inversely proportional to the square root of the cell’s native area. Compressing the system reduces its area and, hence, effectively increases p0. It is, however, important to note that this is just a property of the vertex model and it does not necessarily imply that actual epithelial tissue would behave in the same way. Cells are able to adjust their mechanical properties in response to applied stresses, and it would be overly simplistic to assume that compression would directly lead to changes in the preferred area. In fact, experiments on human bronchial epithelial cells show that applying apical-to-basal compression, which effectively expands the tissue laterally (i.e., corresponds to stretching in our model), fluidizes the tissue [19].

Furthermore, the transition from solid phase to fluid phase is accompanied by the emergence of a large number of soft modes. As we have noted, it has recently been argued that these soft modes lead to a nonlinear response distinct from that obtained in classical models of elasticity [45]. Approximately half of the eigenmodes are zero modes (see Fig H in the Sec H in S1 File). While the analysis of soft modes in the vertex model is an interesting problem [53], it is beyond the scope of this work. Other models in this class have intriguing non-trivial mechanical properties, such as the existence of topologically protected modes [77–81].

We briefly comment on the values of vertex model parameters and timescales that are relevant for experimental systems. While obtaining accurate in vivo measurements of elastic coefficients of epithelial tissues is notoriously difficult, it is possible to make order of magnitude estimates. For example, recent experiments on human corneal epithelial cells estimated K/γ ≈ 0.5 μm−2h−1 and A0 ≈ 500 μm2 [25]. This would correspond to the timescale γ/(KA0) ≈ 15 s, or the relevant frequencies in the ∼ 100−102 Hz range. Characteristic values of stress have been estimated to be KA0 ∼ 10nN/μm for a number of different epithelia [71, 82].

It is also important to note that our work focuses on the rheological behavior in the tangent moduli approximation to a general stress-strain curves, where applied shear and bulk deformations are infinitesimal with respect to a potentially significantly predeformed state (e.g., see Fig 10). This is clearly an idealized case and, in reality, one would be interested in the response to higher values of the applied deformation (≳ 1 − 10%). It would be possible to study such finite deformations using our simulation protocol. Interpreting the results would, however, be more challenging, e.g., due to shear-driven rigidity transition in model tissues and the presence of plastic events [44, 46].

The vertex model provides a rather coarse description of real epithelial tissues and omits many important aspects such as cell polarization, chemical signalling, cells’ ability to actively adjust their properties in response to their environment, etc. The framework presented in this study would, however, be able to address the linear response in the presence of such effects, provided that the vertex model is suitably augmented [83–85].

Regardless of whether cells in an epithelial tissue are arrested or able to move, the rheological response of the tissue is viscoelastic with multiple timescales [38]. This response arises as a result of the complex material properties of individual cells combined with four basic cellular behaviors: movement, shape change, division, and differentiation. The tissue not only has a non-trivial rheological response but is also able to tune it. There is growing evidence that this ability of biological systems to tune their rheology, and in particular, transition between solid-like and fluid-like behaviors, plays a key role during morphogenesis [4]. How such cellular processes are regulated and coordinated to form complex morphological structures is only partly understood. It is, however, clear that the process involves mechano-chemical feedback between mechanical stresses and the expression of genes that control the force-generating molecular machinery in the cell. Any models that aim to describe morphological processes, therefore, need to include coupling between biochemical processes and mechanical responses. The base mechanical model, however, must be able to capture the underlying viscoelastic nature of tissues. Our work provides evidence that the vertex model, a model commonly used to study the mechanics of epithelial tissues, has interesting non-trivial rheological behavior. This, combined with its ability to capture both fluid- and solid-like behavior by tuning a single geometric parameter shows it to be an excellent base model to build more complex descriptions of real tissues.

Supporting information

S1 File. Supporting information file.

https://doi.org/10.1371/journal.pcbi.1010135.s001

(PDF)

Acknowledgments

We would like to acknowledge useful discussions with Ricard Alert, Moumita Das, and Mikko Haataja.

References

  1. 1. Lecuit T, Lenne PF, Munro E. Force generation, transmission, and integration during cell and tissue morphogenesis. Annual Review of Cell and Developmental Biology. 2011;27:157–184. pmid:21740231
  2. 2. Heisenberg CP, Bellaïche Y. Forces in tissue morphogenesis and patterning. Cell. 2013;153(5):948–962. pmid:23706734
  3. 3. Hannezo E, Heisenberg CP. Mechanochemical feedback loops in development and disease. Cell. 2019;178(1):12–25. pmid:31251912
  4. 4. Petridou NI, Heisenberg CP. Tissue rheology in embryonic organization. The EMBO Journal. 2019;38(20):e102497. pmid:31512749
  5. 5. Ross MH, Pawlina W. Histology. Lippincott Williams & Wilkins; 2006.
    • 6. Krndija D, El Marjou F, Guirao B, Richon S, Leroy O, Bellaïche Y, et al. Active cell migration is critical for steady-state epithelial turnover in the gut. Science. 2019;365(6454):705–710. pmid:31416964
    • 7. Wolpert L, Tickle C, Arias AM. Principles of Development. Oxford University Press, USA; 2015.
      • 8. Weinberg RA. The Biology of Cancer. Garland Science; 2013.
        • 9. Forgacs G, Foty RA, Shafrir Y, Steinberg MS. Viscoelastic properties of living embryonic tissues: a quantitative study. Biophysical Journal. 1998;74(5):2227–2234. pmid:9591650
        • 10. Marchetti M, Joanny J, Ramaswamy S, Liverpool T, Prost J, Rao M, et al. Hydrodynamics of soft active matter. Reviews of Modern Physics. 2013;85(3):1143.
        • 11. Friedl P, Gilmour D. Collective cell migration in morphogenesis, regeneration and cancer. Nature Reviews Molecular Cell Biology. 2009;10(7):445–457. pmid:19546857
        • 12. Alert R, Trepat X. Physical models of collective cell migration. Annual Review of Condensed Matter Physics. 2020;11:77–101.
        • 13. Poujade M, Grasland-Mongrain E, Hertzog A, Jouanneau J, Chavrier P, Ladoux B, et al. Collective migration of an epithelial monolayer in response to a model wound. Proceedings of the National Academy of Sciences. 2007;104(41):15988–15993. pmid:17905871
        • 14. Trepat X, Wasserman MR, Angelini TE, Millet E, Weitz DA, Butler JP, et al. Physical forces during collective cell migration. Nature Physics. 2009;5(6):426–430.
        • 15. Tambe DT, Hardin CC, Angelini TE, Rajendran K, Park CY, Serra-Picamal X, et al. Collective cell guidance by cooperative intercellular forces. Nature Materials. 2011;10(6):469–475. pmid:21602808
        • 16. Brugués A, Anon E, Conte V, Veldhuis JH, Gupta M, Colombelli J, et al. Forces driving epithelial wound healing. Nature Physics. 2014;10:683–690. pmid:27340423
        • 17. Etournay R, Popović M, Merkel M, Nandi A, Blasse C, Aigouy B, et al. Interplay of cell dynamics and epithelial tension during morphogenesis of the Drosophila pupal wing. eLife. 2015;4:e07090. pmid:26102528
        • 18. Angelini TE, Hannezo E, Trepat X, Marquez M, Fredberg JJ, Weitz DA. Glass-like dynamics of collective cell migration. Proceedings of the National Academy of Sciences. 2011;108(12):4714–4719. pmid:21321233
        • 19. Park JA, Kim JH, Bi D, Mitchel JA, Qazvini NT, Tantisira K, et al. Unjamming and cell shape in the asthmatic airway epithelium. Nature Materials. 2015;14(10):1040–1048. pmid:26237129
        • 20. Bi D, Lopez J, Schwarz JM, Manning ML. A density-independent rigidity transition in biological tissues. Nature Physics. 2015;11(12):1074–1079.
        • 21. Bi D, Yang X, Marchetti MC, Manning ML. Motility-driven glass and jamming transitions in biological tissues. Physical Review X. 2016;6(2):021011. pmid:28966874
        • 22. Atia L, Bi D, Sharma Y, Mitchel JA, Gweon B, Koehler SA, et al. Geometric constraints during epithelial jamming. Nature Physics. 2018;14(6):613–620. pmid:30151030
        • 23. Sussman DM, Paoluzzi M, Marchetti MC, Manning ML. Anomalous glassy dynamics in simple models of dense biological tissue. EPL (Europhysics Letters). 2018;121(3):36001.
        • 24. Czajkowski M, Sussman DM, Marchetti MC, Manning ML. Glassy dynamics in models of confluent tissue with mitosis and apoptosis. Soft Matter. 2019;15(44):9133–9149. pmid:31674622
        • 25. Henkes S, Kostanjevec K, Collinson JM, Sknepnek R, Bertin E. Dense active matter model of motion patterns in confluent cell monolayers. Nature Communications. 2020;11(1):1–9. pmid:32179745
        • 26. Szabo B, Szöllösi G, Gönci B, Jurányi Z, Selmeczi D, Vicsek T. Phase transition in the collective migration of tissue cells: experiment and model. Physical Review E. 2006;74(6):061908. pmid:17280097
        • 27. Sadati M, Qazvini NT, Krishnan R, Park CY, Fredberg JJ. Collective migration and cell jamming. Differentiation. 2013;86(3):121–125. pmid:23791490
        • 28. Garcia S, Hannezo E, Elgeti J, Joanny JF, Silberzan P, Gov NS. Physics of active jamming during collective cellular motion in a monolayer. Proceedings of the National Academy of Sciences. 2015;112(50):15314–15319.
        • 29. Merkel M, Manning ML. A geometrically controlled rigidity transition in a model for confluent 3D tissues. New Journal of Physics. 2018;20(2):022002.
        • 30. Mitchel JA, Das A, O’Sullivan MJ, Stancil IT, DeCamp SJ, Koehler S, et al. In primary airway epithelial cells, the unjamming transition is distinct from the epithelial-to-mesenchymal transition. Nature Communications. 2020;11(1):1–14. pmid:33028821
        • 31. Bénazéraf B, Francois P, Baker RE, Denans N, Little CD, Pourquié O. A random cell motility gradient downstream of FGF controls elongation of an amniote embryo. Nature. 2010;466(7303):248–252. pmid:20613841
        • 32. Lawton AK, Nandi A, Stulberg MJ, Dray N, Sneddon MW, Pontius W, et al. Regulated tissue fluidity steers zebrafish body elongation. Development. 2013;140(3):573–582. pmid:23293289
        • 33. Mongera A, Rowghanian P, Gustafson HJ, Shelton E, Kealhofer DA, Carn EK, et al. A fluid-to-solid jamming transition underlies vertebrate body axis elongation. Nature. 2018;561(7723):401–405. pmid:30185907
        • 34. Desprat N, Guiroy A, Asnacios A. Microplates-based rheometer for a single living cell. Review of Scientific Instruments. 2006;77(5):055111.
        • 35. Salbreux G, Charras G, Paluch E. Actin cortex mechanics and cellular morphogenesis. Trends in Cell Biology. 2012;22(10):536–545. pmid:22871642
        • 36. Berthoumieux H, Maître JL, Heisenberg CP, Paluch EK, Jülicher F, Salbreux G. Active elastic thin shell theory for cellular deformations. New Journal of Physics. 2014;16(6):065005.
        • 37. Harris AR, Peter L, Bellis J, Baum B, Kabla AJ, Charras GT. Characterizing the mechanics of cultured cell monolayers. Proceedings of the National Academy of Sciences. 2012;109(41):16449–16454. pmid:22991459
        • 38. Bonfanti A, Fouchard J, Khalilgharibi N, Charras G, Kabla A. A unified rheological model for cells and cellularised materials. Royal Society Open Science. 2020;7(1):190920. pmid:32218933
        • 39. Matoz-Fernandez D, Agoritsas E, Barrat JL, Bertin E, Martens K. Nonlinear rheology in a model biological tissue. Physical Review Letters. 2017;118(15):158105. pmid:28452505
        • 40. Nagai T, Honda H. A dynamic cell model for the formation of epithelial tissues. Philosophical Magazine B. 2001;81(7):699–719.
        • 41. Farhadifar R, Röper JC, Aigouy B, Eaton S, Jülicher F. The influence of cell mechanics, cell-cell interactions, and proliferation on epithelial packing. Current Biology. 2007;17(24):2095–2104. pmid:18082406
        • 42. Fletcher AG, Osterfield M, Baker RE, Shvartsman SY. Vertex models of epithelial morphogenesis. Biophysical Journal. 2014;106:2291–2304. pmid:24896108
        • 43. Barton DL, Henkes S, Weijer CJ, Sknepnek R. Active vertex model for cell-resolution description of epithelial tissue mechanics. PLoS Computational Biology. 2017;13(6):e1005569. pmid:28665934
        • 44. Huang J, Cochran JO, Fielding SM, Marchetti MC, Bi D. Shear-driven solidification and nonlinear elasticity in epithelial tissues. Physical Review Letters. 2022; 128(17):178001.
        • 45. Moshe M, Bowick MJ, Marchetti MC. Geometric frustration and solid-solid transitions in model 2D tissue. Physical Review Letters. 2018;120(26):268105. pmid:30004729
        • 46. Popović M, Druelle V, Dye N, Jülicher F, Wyart M. Inferring the flow properties of epithelial tissues from their geometry. New Journal of Physics. 2021;23:033004.
        • 47. Tah I, Sharp T, Liu A, Sussman DM. Quantifying the link between local structure and cellular rearrangements using information in models of biological tissues. Soft Matter. 2021;17(45):10242–10253. pmid:33463648
        • 48. Sanematsu PC, Erdemci-Tandogan G, Patel H, Retzlaff E, Amack JD, Manning L. 3D viscoelastic drag forces drive changes to cell shapes during organogenesis in the zebrafish embryo. Cells & Development. 2021; 168: 203718 pmid:34273601
        • 49. Merkel M, Baumgarten K, Tighe BP, Manning ML. A minimal-length approach unifies rigidity in underconstrained materials. Proceedings of the National Academy of Sciences. 2019;116(14):6560–6568. pmid:30894489
        • 50. Staple D, Farhadifar R, Röper JC, Aigouy B, Eaton S, Jülicher F. Mechanics and remodelling of cell packings in epithelia. The European Physical Journal E. 2010;33(2):117–127. pmid:21082210
        • 51. Bi D, Lopez JH, Schwarz J, Manning ML. Energy barriers and cell migration in densely packed tissues. Soft Matter. 2014;10(12):1885–1890. pmid:24652538
        • 52. Hernandez A, Staddon MF, Bowick MJ, Marchetti MC, Moshe M. Geometric rigidity and anomalous elasticity of cellular tissue vertex model. arXiv preprint arXiv:210910407. 2021.
          • 53. Yan L, Bi D. Multicellular rosettes drive fluid-solid transition in epithelial tissues. Physical Review X. 2019;9(1):011029.
          • 54. Wang X, Merkel M, Sutter LB, Erdemci-Tandogan G, Manning ML, Kasza KE. Anisotropy links cell shapes to tissue flow during convergent extension. Proceedings of the National Academy of Sciences. 2020;117(24):13541–13551. pmid:32467168
          • 55. Sknepnek R, Tong S, Košmrlj A. Rheological Vertex Model (RheoVM); 2022. https://github.com/sknepneklab/RheoVM.
            • 56. Ayachit U. The paraview guide: a parallel visualization application. Kitware, Inc.; 2015.
              • 57. Bitzek E, Koskinen P, Gähler F, Moseler M, Gumbsch P. Structural relaxation made simple. Physical Review Letters. 2006;97(17):170201. pmid:17155444
              • 58. Hočevar A, Ziherl P. Degenerate polygonal tilings in simple animal tissues. Physical Review E. 2009;80(1):011904. pmid:19658726
              • 59. Torquato S, Haslach H Jr. Random heterogeneous materials: microstructure and macroscopic properties. Applied Mechanics Reviews. 2002;55(4):B62–B63.
              • 60. Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D, et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods. 2020;17:261–272. pmid:32015543
              • 61. Purcell EM. Life at low Reynolds number. American Journal of Physics. 1977;45(1):3–11.
              • 62. Walcott S, Sun SX. A mechanical model of actin stress fiber formation and substrate elasticity sensing in adherent cells. Proceedings of the National Academy of Sciences. 2010;107(17):7757–7762. pmid:20385838
              • 63. Sens P. Rigidity sensing by stochastic sliding friction. EPL (Europhysics Letters). 2013;104(3):38003.
              • 64. Schwarz US, Safran SA. Physics of adherent cells. Reviews of Modern Physics. 2013;85(3):1327.
              • 65. Curran S, Strandkvist C, Bathmann J, de Gennes M, Kabla A, Salbreux G, et al. Myosin II controls junction fluctuations to guide epithelial tissue ordering. Developmental Cell. 2017;43(4):480–492. pmid:29107560
              • 66. Leimkuhler B, Matthews C. Molecular Dynamics. Springer; 2016.
                • 67. Chiou KK, Hufnagel L, Shraiman BI. Mechanical stress inference for two dimensional cell arrays. PLoS Computational Biology. 2012;8(5):e1002512. pmid:22615550
                • 68. Yang X, Bi D, Czajkowski M, Merkel M, Manning ML, Marchetti MC. Correlating cell shape and cellular stress in motile confluent tissues. Proceedings of the National Academy of Sciences. 2017;114(48):12663–12668. pmid:29138312
                • 69. Nestor-Bergmann A, Goddard G, Woolner S, Jensen OE. Relating cell shape and mechanical stress in a spatially disordered epithelium using a vertex-based model. Mathematical Medicine and Biology: A Journal of the IMA. 2018;35(Supplement_1):i1–i27.
                • 70. Larson RG. The Structure and Rheology of Complex Fluids. vol. 150. Oxford University Press New York; 1999.
                  • 71. Wyatt TP, Fouchard J, Lisica A, Khalilgharibi N, Baum B, Recho P, et al. Actomyosin controls planarity and folding of epithelia in response to compression. Nature Materials. 2020;19(1):109–117. pmid:31451778
                  • 72. Murisic N, Hakim V, Kevrekidis IG, Shvartsman SY, Audoly B. From discrete to continuum models of three-dimensional deformations in epithelial sheets. Biophysical Journal. 2015;109(1):154–163. pmid:26153712
                  • 73. Gibson MC, Patel AB, Nagpal R, Perrimon N. The emergence of geometric order in proliferating metazoan epithelia. Nature. 2006;442(7106):1038–1041. pmid:16900102
                  • 74. Sussman DM, Merkel M. No unjamming transition in a Voronoi model of biological tissue. Soft Matter. 2018;14(17):3397–3403. pmid:29667689
                  • 75. Noll N, Mani M, Heemskerk I, Streichan SJ, Shraiman BI. Active tension network model suggests an exotic mechanical state realized in epithelial tissues. Nature Physics. 2017;13(12):1221–1226. pmid:30687408
                  • 76. Tong S, Sknepnek R, Košmrlj A. Normal mode analysis of the linear viscoelastic response of dissipative systems: application to vertex model. arXiv preprint arXiv:220203261. 2022.
                    • 77. Kane C, Lubensky T. Topological boundary modes in isostatic lattices. Nature Physics. 2014;10(1):39–45.
                    • 78. Lubensky T, Kane C, Mao X, Souslov A, Sun K. Phonons and elasticity in critically coordinated lattices. Reports on Progress in Physics. 2015;78(7):073901. pmid:26115553
                    • 79. Paulose J, Chen BGg, Vitelli V. Topological modes bound to dislocations in mechanical metamaterials. Nature Physics. 2015;11(2):153–156.
                    • 80. Rocklin DZ, Zhou S, Sun K, Mao X. Transformable topological mechanical metamaterials. Nature Communications. 2017;8(1):1–9. pmid:28112155
                    • 81. Mao X, Lubensky TC. Maxwell lattices and topological mechanics. Annual Review of Condensed Matter Physics. 2018;9:413–433.
                    • 82. Vincent R, Bazellières E, Pérez-González C, Uroz M, Serra-Picamal X, Trepat X. Active tensile modulus of an epithelial monolayer. Physical Review Letters. 2015;115(24):248103. pmid:26705659
                    • 83. Lan H, Wang Q, Fernandez-Gonzalez R, Feng JJ. A biomechanical model for cell polarization and intercalation during Drosophila germband extension. Physical Biology. 2015;12(5):056011. pmid:26356256
                    • 84. Blanchard GB, Fletcher AG, Schumacher LJ. The devil is in the mesoscale: mechanical and behavioural heterogeneity in collective cell movement. In: Seminars in Cell & Developmental Biology. vol. 93. Elsevier; 2019. p. 46–54.
                      • 85. Sknepnek R, Djafer-Cherif I, Chuai M, Weijer CJ, Henkes S. Generating active T1 transitions through mechanochemical feedback. arXiv preprint arXiv:210612394. 2021.