Linear Viscoelastic Properties of the Vertex Model for Epithelial Tissues
Sijie Tong, Navreeta K. Singh, Rastko Sknepnek, Andrej Kosmrlj
LLinear Viscoelastic Properties of the Vertex Model for Epithelial Tissues
Sijie Tong, Navreeta K. Singh, Rastko Sknepnek,
2, 3, ∗ and Andrej Koˇsmrlj
1, 4, † Department of Mechanical and Aerospace Engineering,Princeton University, Princeton, New Jersey 08544, USA School of Science and Engineering, University of Dundee, Dundee DD1 4HN, United Kingdom School of Life Sciences, University of Dundee, Dundee DD1 5EH, United Kingdom Princeton Institute for the Science and Technology of Materials (PRISM),Princeton University, Princeton, New Jersey 08544, USA
Epithelial tissues act as barriers and, therefore, must repair themselves, respond to environmentalchanges and grow without compromising their integrity. Consequently, they exhibit complex vis-coelastic rheological behavior where constituent cells actively tune their mechanical properties tochange the overall response of the tissue, e.g., from fluid-like to solid-like. Mesoscopic mechanicalproperties of epithelia are commonly modeled with the vertex model (VM). We systematically stud-ied rheological properties of the VM by applying small oscillatory shear and bulk deformations andmeasuring the response stresses. We found that the shear and bulk responses in the fluid and solidphases can be described by standard spring-dashpot viscoelastic models. Further, we found thatthe solid-fluid transition can be tuned by applying pre-deformation to the system. Thus, the VM’srich rheological behavior makes it suitable for describing mesoscale dynamics of epithelial tissuesbeyond the quasistatic limit.
INTRODUCTION
The development and maintenance of tissues requiresclose coordination of biochemical and mechanical signal-ing [1–3]. There is, for instance, mounting evidence forthe key role played by tissue material properties andtheir regulation during embryonic development [4]. Tis-sues must be able to adjust their mechanical propertiesin response to internal and external stimuli. In par-ticular, epithelial tissues, which line all cavities in thebody and demarcate organs, must sustain substantialmechanical stresses while also supporting numerous bi-ological processes such as selective diffusion and absorp-tion/secretion [5]. In homeostasis, the tissue must main-tain its shape and resist deformation while remainingflexible. The tissue must also be able to regenerate andrepair itself, often with fast turnover, e.g., in gut ep-ithelia [6]. Furthermore, in morphogenesis, the tissuemust take up a specific shape and function [7]. Duringmetastasis, however, the shape is lost and cancer cellsinvade surrounding healthy tissue [8]. All of these pro-cesses require that cells be able to move, often over dis-tances much larger than the cell size. During cell mi-gration, however, the tissue must maintain its integrity.It is, therefore, not surprising that tissues exhibit richviscoelastic behavior [9]. Unlike passive viscoelastic ma-terials, a tissue can actively tune its rheological response,making the study of its rheology not only important forunderstanding biological functions but also an interest-ing problem from the perspective of the physics of activematter systems [10].Collective cell migration has been extensively studiedin biology [11] and biophysics [12]. In vitro studies of ∗ [email protected] † [email protected] confluent cell monolayers [13–17] focused on the phys-ical aspects of force generation and transmission andshowed that cell migration is an inherently collectivephenomenon. Some aspects of collective cell migrationare remarkably similar to slow dynamics of structuralglasses [18–24]. This suggests that many of the observedbehaviors share common underlying mechanisms and canbe understood, at least at mesoscales (i.e., distances be-yond several cell diameters), using physics of dense ac-tive systems [25]. A particularly intriguing observation isthat tuning cell density [18, 26, 27], strength of cell-celland cell-substrate interactions [28], or cell shape parame-ters [21, 29] can cause the collective migration to stop. Inother words, the tissue undergoes fluid to solid transition.Signatures of such behavior have been reported in severalin vitro [19, 30] and developmental systems [31–33]. Thissuggests that important aspects of morphogenetic devel-opment might rely on tissue’s ability to undergo phasetransitions [4].How a tissue responds to external and internal mechan-ical stresses will depend on its rheological (i.e., material)properties. While there have been numerous studies fo-cusing on the rheology of a single cell [34–36], much lessis known about tissue rheology, particularly during de-velopment. In order to develop a comprehensive under-standing of tissue mechanics, such insight is key. Thoughsingle cell measurements are valuable, the mechanics ofa tissue can be drastically different from that of its con-stituent cells. The stiffness of cell monolayers, for ex-ample, is orders of magnitude higher than the stiffnessof constituent cells, while the time dependent mechani-cal behaviors of monolayers in response to deformationvary depending on the magnitude of loading [37]. Em-bryonic cell aggregates have been shown to behave elas-tically (i.e., solid-like) at short time scales but flow likefluids at long time scales, which facilitates both the ro-bustness needed to maintain integrity and the flexibil- a r X i v : . [ c ond - m a t . s o f t ] F e b ity to morph during development [9]. Experiments havecharacterized the mechanical behaviors of tissues at var-ious loading conditions, which led to a phenomenologi-cal description that models the relaxation properties ofepithelial monolayers based on fractional calculus [38].Notably, a recent particle-based model that includes celldivision and apoptosis provided a plausible microscopicmodel for nonlinear rheological response [39]. Particle-based models are, however, unable to capture geometricaspects such as cell shape. It is, therefore, necessary toinvestigate rheological response in geometric models.The vertex model (VM) [40–42] and more recent,closely related self-propelled Voronoi models (SPV) [21,43] have played an important role in modeling mechanicsof epithelial tissues since they account for the shapes ofindividual 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 solidto fluid transition and demonstrate rich and unusual non-linear mechanical behavior [44, 45]. While the mechani-cal properties of the VM and SPV have been extensivelystudied, most works to date focused on the long-timebehavior, e.g., by studying the quasistatic shear modu-lus [20], effective diffusion constant of cells that is relatedto the viscosity of tissues [21], and correlations betweena structural property called “softness” and the likelihoodof topological rearrangements of cells [46]. The rheolog-ical properties of the VM that cover a broad range oftime scales, however, have not yet been systematicallyexplored. In this paper, we present a detailed study ofthe rheology of the VM by studying its response to ap-plied oscillatory shear and bulk deformations of smallamplitude, i.e., in the linear response regime. We mea-sured the response stresses and used them to computethe storage and loss moduli in both the solid and fluidphases. We show that the dynamical response of theVM can be fitted to standard spring-dashpot viscoelasticmodels over seven decades in the driving frequency andthat the solid-fluid transition can be tuned by applyingpre-deformation to the system. Thus we argue that theVM makes a suitable basis for studies of tissue dynamicsbeyond the quasistatic limit. RESULTS
Vertex model dynamics.
In the VM, the state of anepithelial tissue is approximated as a polygonal tiling ofthe plane. The degrees of freedom are vertices, i.e., meet-ing points of three or more cell-cell junctions. In the sim-plest formulation, junctions are assumed to be straightlines. The energy of the VM is a quadratic function ofcell areas and perimeters [41], i.e., E = (cid:88) C (cid:20) K C A C − A C ) + Γ C P C − P C ) (cid:21) , (1)where K C and Γ C are the area and perimeter elastic mod-uli. A C and A C are the actual and preferred areas of cell C , respectively. Similarly, P C and P C are, respectively,the actual and preferred perimeters of the same cell. Inthis work, we assumed K C , Γ C , A C , and P C to be iden-tical for all cells (i.e., K C ≡ K, Γ C ≡ Γ , A C ≡ A , P C ≡ P ). Furthermore, we fixed the values of K and A , andmeasured the energy in units of KA , stresses in unitsof KA , and lengths in units of A / . Since the ratiobetween the perimeter and area elastic moduli does notqualitatively change the behavior of the VM [20, 41], wefixed that ratio to Γ / ( KA ) ≈ . P , which sets the dimensionless cell-shape parameter,defined as the ratio p = P / √ A .The cell-shape parameter, p , plays a central role indetermining whether the system behaves as a fluid orsolid [20]. Bi, et al . [20] argued that the rigidity transitionoccurs at p = p c ≈ .
812 for a random polygonal tiling.The transition point is, however, at p c = (cid:112) √ ≈ . p is reduced below p c , the energy barrier becomes finite, neighbor exchangescease and the system becomes solid. While the transitionpoint for hexagonal tilings can be understood in terms ofthe mechanical stability and the excess perimeter [44],the mechanism that leads to a larger value for randomtilings is more subtle and not fully understood [49].For definiteness, in all simulations we always startedfrom a regular hexagonal tiling in a nearly square boxsubject to periodic boundary conditions, where the initialcell areas A C matched the preferred areas A (see Meth-ods). For the solid phase with p (cid:46) . P C from the preferredvalues P ), which could be eliminated by the appropriaterelaxation of the simulation box. This hydrostatic stress,however, does not qualitatively affect the rheological be-havior of the system (see Supplementary Information,Sec. IV for further discussion). For the fluid phase with p (cid:38) . . × − √ A (see Meth-ods) and then the system was relaxed using the FIRE al-gorithm [50] to reach a local energy minimum. Note thatthe energy landscape in the fluid phase has many localminima and a large number of soft modes (see SupportingInformation, Sec. VIII). Thus we repeated simulations toinvestigate rheological properties for multiple configura-tions corresponding to different local energy minima (seeMethods).In order to probe the dynamical response of the VM,we need to specify the microscopic equations of motionfor vertices. Assuming the low Reynolds number limit,which is applicable to most cellular systems due to their FIG. 1. Loss and storage shear moduli in the solid (top row) and fluid phase (bottom row). An overlay of the representativereference (grey) and sheared (yellow) configurations in (a) the solid and (b) the fluid phase. The magnitude of the shear is highlyexaggerated for demonstration purposes. (c-d) show the representative storage ( G (cid:48) ) and loss ( G (cid:48)(cid:48) ) shear moduli as functions ofthe shearing frequency, ω , for different values of the cell-shape parameter, p , where we removed the dissipative contributionto G (cid:48)(cid:48) due to the velocity field resulting from the affine part of the deformation that dominates at large frequencies (see text).Dashed curves are the fits based on (c) the Standard Linear Solid (SLS) model in the solid phase and (d) the Burgers model inthe fluid phase. (e-f) The collapse of the moduli curves for different values of p 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. Themajority of the data corresponds to the system of nearly square shape with N x = 15 cells in the horizontal direction, and wealso show examples of larger systems with N x = 37 and N x = 51 cells in the horizontal direction. slow speed, inertial effects can be neglected [51]. Theequations of motion are then a force balance between fric-tion and elastic forces due to deformations of cell shapes,i.e., γ ˙ r i = F i . (2)Here, r i is the position vector of vertex i in a laboratoryframe of reference, F i = −∇ r i E is the mechanical forceon vertex i due to deformation of cells surrounding it, γ is the friction coefficient, and dot denotes the timederivative. In simulations we fixed the value of γ , whichsets the unit of time as γ/ ( KA ).We note that a precise model for dissipation in ep-ithelial tissues is at present not known. While havinga mechanism of dissipating the external energy input iscentral in this study, its precise microscopic details are,however, not important. We therefore assume that eachvertex experiences dissipative drag proportional to itsinstantaneous velocity, which is a common assumptionin discrete models of soft materials [52]. Furthermore,we neglected thermal fluctuations of vertices and henceomit the stochastic term in Eq. (2). This is a reasonable assumption since typical energy scales in tissues signif-icantly exceed the thermal energy, k B T , at room tem-perature T , where k B is the Boltzmann constant. It isimportant, however, to note that in epithelia there areother sources of stochasticity (e.g., fluctuations of thenumber of force-generating molecular motors) that areimportant for tissue scale behaviors [53]. Here, we donot consider such effects but point out that they couldbe directly included in the model as additional forces inEq. (2). Response to a shear deformation.
The hexagonalground state in the solid phase and states correspond-ing to local energy minima in the fluid phase were thenused to investigate the rheological behavior by applyingan oscillatory affine shear deformation (Fig. 1a,b). Ateach time step, we first applied the affine shear defor-mation to the simulation box and all vertices that wasfollowed by internal relaxation of vertices according toEq. (2) (see Methods). The affine shear deformation canbe described by a deformation gradient tensor defined asˆ F = ∂ x /∂ X , where the mapping x = x ( X , t ) mapsthe reference configuration X to a spatial configuration x at time t . For simple shear, the deformation gradienttensor is ˆ F = (cid:0) (cid:15) ( t )0 1 (cid:1) , where (cid:15) ( t ) = (cid:15) sin ( ω t ). Suffi-ciently small amplitude (cid:15) = 10 − (cid:28) σ C ( t ), for each cell C wascomputed using the formalism introduced in [54] (seeSupplementary Information, Sec. I). The average stresstensor ˆ σ ( t ) = (cid:80) C w C ˆ σ C ( t ), with w C = A C / (cid:80) C A C ,was used as a measure for the response of the system.The dynamic shear modulus G ∗ ( ω ) = ˜ τ ( ω ) / ˜ (cid:15) ( ω ) wasthen calculated at a given frequency ω of applied shearstrain, where ˜ τ ( ω ) and ˜ (cid:15) ( ω ) are the Fourier transformsof the response shear stress τ ( t ) = ˆ σ xy ( t ) and the ap-plied strain (cid:15) ( t ), respectively (see Methods). We ensuredthat the simulations were sufficiently long for a steadystate to be reached (see Methods and Supplementary In-formation, Sec. III). The real part of the dynamic shearmodulus, G (cid:48) = Re ( G ∗ ), is the storage shear modulus andthe imaginary part, G (cid:48)(cid:48) = Im ( G ∗ ), is the loss shear mod-ulus. 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 cor-responds to the out-of-phase response and measures theirreversible dissipation [55] (see also Supplementary In-formation, Sec. II). Note that the stress contribution dueto the friction with the substrate resulting from the partof the velocity field due to the externally imposed affinedeformation was removed in the analysis. This is be-cause at high frequencies, the stresses are dominated bythe friction due to external driving that masks the inter-nal material response that we are interested in exploring(see Fig. S3 and Supplementary Information, Sec. I). Forsystems under an oscillatory simple shear, storage andloss shear moduli were obtained for different values of p and different system sizes in the solid and the fluid phaseswith frequencies ω of the applied shear strain spanningover seven orders of magnitude, as shown in Fig. 1c,d.Most simulations were performed for systems with nearlysquare shapes with N x = 15 cells in the horizontal direc-tion. We repeated several simulations for systems with N x = 37 and N x = 51, which showed that the finite sizeeffects are negligible (Fig. 1c-f).In the solid phase there are two different regimes (seeFig. 1c). At low frequencies, ω , the storage shear modu-lus G (cid:48) has a constant value, while the loss shear modulusscales as G (cid:48)(cid:48) ∝ ω . At high frequencies the storage shearmodulus G (cid:48) has a higher constant value, while the lossshear modulus scales as G (cid:48)(cid:48) ∝ ω − . Such rheologicalbehavior is characteristic for the Standard Linear Solid(SLS) model [55]. Storage and loss shear moduli for theSLS model are [55], respectively, G (cid:48) SLS ( ω ) = E + η E ω ( E + E )1 + η E ω , (3a) G (cid:48)(cid:48) SLS ( ω ) = ω η η E ω , (3b) where we used the representation of the SLS model(Fig. 1e, inset) that consists of a spring with elastic con-stant E connected in parallel with a Maxwell element,which comprises a spring with elastic constant E and adashpot with viscosity η connected in series. The aboveexpressions in Eqs. (3) were used to fit the storage andloss shear moduli obtained from simulations. The fittedcurves, represented with dashed lines in Fig. 1c, show anexcellent match with the simulation data, indicating thatthe SLS model is indeed appropriate to describe the shearrheology in the solid phase. This was also confirmed inFig. 1e, where we collapsed the storage and loss shearmoduli for different values of the shape parameter, p ,by rescaling the moduli and frequencies with the fittedvalues of spring and dashpot constants.As the value of the p increases, we observe that thestorage shear modulus reduces at all frequencies and thatthe loss shear modulus reduces at high frequencies. Fur-thermore the crossover between the two regimes shiftstowards lower frequencies (Fig. 1c). This is because theelastic constants E and E decrease linearly with in-creasing p and they become zero exactly at the solid-fluid transition with p = p c ≈ .
722 (Fig. 2a). The dash-pot constant η is nearly independent of p and scaleswith the friction parameter γ , which is the only sourceof dissipation in the VM. The crossover between the tworegimes for both the storage and loss shear moduli cor-responds to a characteristic time scale, η /E , which di-verges as p approaches the solid-fluid transition due tothe vanishing elastic constant (Fig. 2c). Note that thevalues of the elastic constants E and E can be esti-mated analytically. In the quasistatic limit ( ω → et al. showedthat the storage shear modulus is G (cid:48) ( ω →
0) = E = 12 KA (cid:0) − [ α ( p , Γ /KA )] (cid:1) , (4)where α ( p , Γ /KA ) is a scaling factor chosen such thatthe hydrostatic stress vanishes once the system box sizeis rescaled from L to αL [56]. In the high frequency limit( ω → ∞ ), on the other hand, the system follows theexternally imposed affine deformation and has no timefor internal relaxation. Thus, by considering the energycost for a hexagonal tiling under affine deformation, weobtained the storage shear modulus (see SupplementalInformation, Sec. VII) G (cid:48) ( ω → ∞ ) = E + E = 3 √ (cid:18) − p p c (cid:19) . (5)The above Eqs. (4) and (5) were used to extract the val-ues of elastic constants E and E , which showed excel-lent agreement with the fitted values from simulations(Fig. 2a).In the fluid phase, the storage and loss shear modulishow a markedly different behavior (Fig. 1d). There arethree different regimes with two crossover frequencies, FIG. 2. (a-b) Fitted values of spring-dashpot models for thesystem under simple shear. (a) Elastic constants as a func-tion of target cell-shape parameter, p . In the solid phase(i.e., for p < p c ≈ . p . (c-d) Characteristic time scales in (c) the solid and(d) fluid phase obtained from the fitted values of the elasticconstant and the dashpot viscosity. The normalization factor t ∗ = γ/ ( KA ) sets the unit of time. For the fluid phase (i.e.,for p > p c ≈ . which correspond to two characteristic time scales. Atlow frequencies, ω , the storage shear modulus G (cid:48) ∝ ω and the loss shear modulus G (cid:48)(cid:48) ( ω ) ∝ ω . The storagemodulus approaches 0 for ω →
0, which indicates thatthe system is indeed a fluid. At high frequencies the stor-age shear modulus has a constant value, while the lossshear modulus scales as G (cid:48)(cid:48) ( ω ) ∝ ω − . To capture thisbehavior we used the Burgers model, which consists oftwo Maxwell models connected in parallel (Fig. 1f, in-set), to fit the shear moduli measured in the simulations.The storage and loss shear moduli for a Burgers modelare [55], respectively, G (cid:48) Burg ( ω ) = p q ω − q ω (cid:0) − p ω (cid:1) p ω + (1 − p ω ) , (6a) G (cid:48)(cid:48) Burg ( ω ) = p q ω + q ω (cid:0) − p ω (cid:1) p ω + (1 − p ω ) , (6b)where p = η /E + η /E , p = η η / ( E E ), q = η + η , q = η η ( E + E ) / ( E E ). The dashed curvesin Fig. 1d show fits of the storage and loss shear modulifor a range of values of p , which show good agreementwith simulations. Unlike for the solid phase, it is not pos-sible to collapse the data for storage and loss shear mod-uli onto single universal curves because the fluid phase is characterized by two independent timescales η /E and η /E . Thus we show two different collapses for the stor-age and loss shear moduli in the low frequency range(Fig. 1f) and in the high frequency range (Fig. S7 in theSupplementary Information, Sec. V).As the value of the p decreases, we observe that boththe storage and loss shear moduli reduce at intermediateand high frequencies, but they increase at low frequencies(Fig. 1d). We also observe that the first crossover shiftstowards lower frequencies, while the second crossover re-mains at approximately the same frequency. This is be-cause the elastic constants E and E decrease linearlytoward zero as p approaches the solid-fluid transition at p c ≈ .
722 (Fig. 2a). The dashpot constant η also de-creases linearly toward zero, while the dashpot constant η increases but remains finite as p approaches the solid-fluid transition (Fig. 2b). As a consequence, one of thecharacteristic time scales η /E diverges, while the sec-ond time scale η /E remains finite as as p approachesthe solid-fluid transition (Fig. 2d). The diverging char-acteristic time scale captures the macroscopic behaviorof the system, while the second time scale ( ∼ γ/KA )captures the microscopic details of the VM. Finally, wenote that the values of the spring and dashpot constantsare somewhat sensitive to the local energy minimum con-figuration used to probe the response in the fluid phase.The errorbars in Fig. 2 show standard deviation for differ-ent configurations that were obtained by using the samemagnitude of the initial perturbation (see Methods). InFig. S8 in the Supplementary Information, Sec. VI, weshow how the values of the spring and dashpot constantsare affected when configurations were obtained by usingdifferent magnitudes of the initial perturbation. Response to bulk deformations.
We further studiedthe bulk rheological properties of the system by applyingan oscillatory biaxial deformation (Fig. 3a,b) describedby the deformation gradient ˆ F = (cid:0) (cid:15) ( t ) 00 1+ (cid:15) ( t ) (cid:1) , where (cid:15) ( t ) = (cid:15) sin ( ω t ). We applied a sufficiently small am-plitude (cid:15) = 10 − (cid:28) σ ( t ) = [ˆ σ xx ( t ) + ˆ σ yy ( t )], where we again removed thestress contribution due to the friction with the substrateresulting from the part of the velocity field due to theexternally imposed affine deformation. As in the simpleshear test, we then computed the dynamic bulk modulusas B ∗ ( ω ) = ˜ σ ( ω ) / ˜ (cid:15) ( ω ) from which we obtained thestorage bulk modulus B (cid:48) = Re ( B ∗ ) and the loss bulkmodulus B (cid:48)(cid:48) = Im ( B ∗ ) (see Fig. 3c,d).In the solid phase, the storage bulk modulus is inde-pendent of the driving frequency and the loss bulk mod-ulus is zero. Thus the response of the system can becaptured by a single spring E solid (Fig. 3e, inset). Thisis because the hexagonal tiling is stable to biaxial defor-mation in the solid phase and thus there is no internalrelaxation. The measured value of the storage bulk mod-ulus matches the analytical prediction B theory = 2 KA + √ p (7) FIG. 3. Loss and storage bulk moduli in the solid (top row) and fluid phase (bottom row). An overlay of the representativereference (grey) and biaxially deformed (yellow) configurations in (a) the solid and (b) the fluid phase. The magnitude ofthe bulk deformation is highly exaggerated for demonstration purposes. (c-d) show the representative storage ( B (cid:48) ) and loss( B (cid:48)(cid:48) ) bulk moduli as functions of the deformation frequency, ω , for different values of the cell-shape parameter, p , wherewe removed the dissipative contribution to B (cid:48)(cid:48) due to the velocity field resulting from the affine part of the deformation thatdominates at large frequencies (see text). For the solid phase in (c), the loss bulk modulus B (cid:48)(cid:48) ≡
0. For the fluid phase in(d), dashed curves are the fits based on the Standard Linear Solid (SLS) model. (e-f) The collapse of the moduli curves fordifferent values of p for (e) the solid phase and (f) the fluid phase. The insets show the representation of (e) the spring modeland (f) the SLS model in terms of the springs and dashpots. by Staple, et al. [47], where the hexagonal tiling is as-sumed to undergo affine deformation under biaxial de-formation. Storage bulk moduli, normalized by B theory ,for different values of p all collapse to 1 (Fig. 3e).In the fluid phase, the bulk response behavior of thesystem can be described by the SLS model (Fig. 3f, in-set). While it might appear counter-intuitive to modela fluid with the SLS model, this is a direct consequenceof the fact that in the fluid state, the bulk modulus isfinite but the shear modulus vanishes, i.e., the fluid flowsin response to shear but resists bulk deformation. Thefitted storage and loss bulk moduli for the SLS model[see Eq. (3)] show an excellent match with the simula-tion data (Fig. 3d). This was also confirmed in Fig. 3f,where we collapsed the storage and loss bulk moduli fordifferent values of p .The fitted values of spring elastic and dashpot vis-cosity constants for different values of p are plotted inFig. 4. In the fluid phase, the storage bulk modulus in thehigh frequency limit ( B (cid:48) ( ω → ∞ ) = E + E ) continu-ously increases from the value for the solid phase B theory in Eq. (7) as the system transitions from solid to fluid(Fig. 4a). The storage bulk modulus in the quasistatic FIG. 4. Fitted values of spring-dashpot models for the systemunder bulk deformation as a function of the target cell-shapeparameter, p . (a) In the solid phase ( p < p c ≈ . E solid agrees with the analyticalprediction B theory in Eq. (7). At the solid-fluid transitionpoint ( p = p c ≈ . B (cid:48) ( ω →∞ ) = E + E , of the fluid phase. The low frequency limitof the bulk storage modulus is B (cid:48) ( ω →
0) = E in the fluidphase. (b) Dashpot viscosity constant as a function of p .For the fluid phase ( p > p c ≈ . solid Hessiantheory fluid
FIG. 5. Tuning the solid to fluid transition by applying uniaxial deformation. (a) The solid-fluid transition boundary in the a − p plane, where a measures the amount of uniaxial deformation described by the deformation gradient ˆ F = (cid:0) a
00 1 (cid:1) . Blueline shows the analytical prediction from Eq. (9), 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 and the Burgers modelin the fluid phase when the system is under uniaxial compression ( a = 0 . a = 1), and under uniaxialtension ( a = 1 . limit ( B (cid:48) ( ω →
0) = E ) emerges at the transition pointwith a finite value and increases as p increases from p c (Fig. 4a). Fig. 4b shows that the dashpot constant di-verges as the p decreases toward p c . Thus, the charac-teristic time scale η /E also diverges, but for a differentreason than for the shear deformation, where the springconstant is vanishing (see Fig. 2). Finally, we note that,unlike for the response to shear, the values of the springand dashpot constants for bulk deformation are not sen-sitive to the local energy minimum configuration used toprobe the response in the fluid phase, which is reflectedin very small errorbars in Fig. 4. This is because the bulkmoduli are dominated by the changes in cell areas. Response to a shear deformation of a uniaxi-ally pre-deformed system.
The solid-fluid transitionfor the regular hexagonal tiling occurs when p ≈ . a , which is described by the deformation gradientˆ F = (cid:0) a
00 1 (cid:1) , then the high frequency limit of the storageshear modulus that is dominated by affine deformationbecomes (see Supplemental Information, Sec. VII), G (cid:48) affine ( a ) = 2 √ / a (cid:32) a ) (cid:33) × (cid:16) − p + √ (cid:16) (cid:112) a (cid:17)(cid:17) . (8)By setting the affine shear modulus to 0, we obtained thesolid-fluid transition boundary in the a − p plane as a ( p ) = (cid:115) √ p − p √ . (9)The above analytical prediction for the phase boundary(Fig. 5a, blue line) shows an excellent agreement withthe stability analysis in terms of the eigenvalues of the Hessian matrix of the energy function (Fig. 5a, red dots).A given configuration is stable if all eigenvalues of theHessian matrix are positive and the loss of mechanicalstability occurs when the lowest eigenvalue becomes 0.For a given p , the value of the lowest eigenvalue reduceswith decreasing a , i.e., as the magnitude of compression isincreased. Thus, the compression (stretching) shifts thesolid-fluid transition towards the lower (higher) values of p (see Fig. 5a).We also probed the response to oscillatory shear ap-plied to uniaxially compressed and stretched systems.This analysis was done on both the uniaxially deformedhexagonal tiling in the solid phase as well as a systemin the fluid phase obtained by relaxing the unstable uni-axially deformed hexagonal tiling after an initial randomperturbation (see Methods). The response to the sheardeformation is qualitatively similar and can still be de-scribed by the SLS model in the solid phase and theBurgers model in the fluid phase. Fig. 5b,c shows fittedvalues of the parameters for spring-dashpot models whenthe system is under uniaxial compression ( a = 0 . a = 1, i.e., same as Fig. 2a,b), and uni-axial tension ( a = 1 . p approaches the critical value predicted by Eq. (9). Thedashpot viscosity remains constant in the solid phase.Once the system enters the fluid phase as p increases, anew second dashpot constant emerges and increases from0, while the value of the first dashpot constant decreases.As in the simple shear case, we note that the values ofthe spring and dashpot constants are somewhat sensitiveto the local energy minimum configuration used to probethe response in the fluid phase. The errorbars in Fig. 2show standard deviation for configurations that were ob-tained by using different random initial perturbation (seeMethods). Finally, we note that besides uniaxial defor-mation, the solid-fluid transition point can be tuned byother modes of deformation (see Fig. S9 and Supplemen-tary Information, Sec. VII). DISCUSSION
We have performed a detailed analysis of the rheo-logical properties of the Vertex Model subject to small-amplitude oscillatory deformations over seven orders ofmagnitude in the driving frequency. Our analysis showsthat the VM exhibits non-trivial viscoelastic behaviorthat can be tuned by a single dimensionless geometricparameter - the shape parameter, p . In order to char-acterize the response, we constructed constitutive rhe-ological models that use combinations of linear springsand dashpots connected in series and in parallel. Thesemodels allowed us to match the shear response of theVM to that of the Standard Linear Solid model in thesolid phase and the Burgers model in the fluid phase.In the low-frequency, i.e., quasistatic regime, our re-sults are fully consistent with many previous studies[20, 21, 47, 56]. Our work, however, provides insightsinto the time-dependent response of the VM over a broadrange of driving frequencies, which is important if one isto develop full understanding of the rheological proper-ties of the VM and how those inform our understandingof the rheology of epithelial tissues.We also showed that the critical value for the solid-fluidtransition can be tuned by applying pre-deformation. In-terestingly, under uniaxial and biaxial (i.e., isotropic)compression the solid to fluid transition shifts to lowervalues of p , leading to the non-intuitive prediction thatone can fluidize the system by compressing it. This is,however, not surprising since the transition is driven bya geometric parameter that is inversely proportional tothe square root of the cell’s native area. Compressing thesystem reduces its area and, hence, effectively increases p . It is, however, important to note that this is justa property of the VM and it does not necessarily implythat actual epithelial tissue would behave in the sameway. Cells are able to adjust their mechanical propertiesin response to applied stresses and it would be too sim-plistic to assume that compression would directly lead tochanges in the native area. In fact, experiments on hu-man bronchial epithelial cells show that applying apical-to-basal compression, which effectively expands the tis-sue laterally (i.e., corresponds to stretching in our model)fluidizes the tissue [19].Furthermore, the transition from solid phase to fluidphase is accompanied by the emergence of a large numberof soft modes. As we have already alluded, it has recentlybeen argued that these soft modes lead to a nonlinear re-sponse distinct from that obtained in classical models ofelasticity [44]. Approximately half of the eigenmodes arezero modes (see Fig. S10 in the Supplementary Informa-tion, Sec. VIII). While the analysis of soft modes in theVM is an interesting problem [49], it is beyond the scopeof this work. Other models in this class have intriguingnon-trivial mechanical properties, such as the existence of topologically protected modes [57–62].We note that while we studied the dynamical responseover a wide range of frequencies, our work focused onthe behavior in the linear response regime, where thereare effectively no plastic events, i.e., while being allowed,T1 transformations typically did not occur during theprocess of probing the rheology. A full understanding ofthe VM rheology would also need to allow for cell rear-rangements. This is, however, a very challenging problemand first steps in addressing it have only recently beenmade [45].Regardless of whether cells in an epithelial tissue arearrested or able to move, the rheological response of thetissue is viscoelastic with multiple time scales [38]. Thisresponse arises as a result of the complex material prop-erties of individual cells combined with four basic cel-lular behaviors: movement, shape change, division, anddifferentiation. The tissue not only has a non-trivial rhe-ological response but is also able to tune it. There isgrowing evidence that this ability of biological systemsto tune their rheology, and in particular, transition be-tween solid-like and fluid-like behaviors, plays a key roleduring morphogenesis [4]. How such cellular processesare regulated and coordinated to form complex morpho-logical structures is only partly understood. It is, how-ever, clear that the process involves mechano-chemicalfeedback between mechanical stresses and the expressionof genes that control the force-generating molecular ma-chinery in the cell. Any models that aim to describemorphological processes, therefore, need to include cou-pling between biochemical processes and mechanical re-sponse. The base mechanical model, however, must beable to capture the underlying viscoelastic nature of tis-sues. Our work provides evidence that the vertex model,a commonly used model to study mechanics of epithelialtissues, has interesting non-trivial rheological behavior.This, combined with its ability to capture both fluid-and solid-like behavior by tuning a single geometric pa-rameter shows it to be an excellent base model to buildmore complex descriptions of real tissues. METHODS
Most simulations were performed for systems withnearly square shapes with N x = 15 cells in horizontaldirection (i.e., N = 240 cells in total) subject to periodicboundary conditions. The shape of the simulation boxwas chosen to be as near to a square as allowed by the ge-ometry of a hexagon and the size of the box was such thatit accommodated N cells of area A C that matched thepreferred areas A . Simulations of larger system sizes( N x = 37, 51, i.e., N = 1406, 2652 cells, respectively)were performed for a subset of values of p to explore thefinite size effects. No quantitative differences between thesystem with N = 240 cells and larger systems were ob-served. Since the ratio between the perimeter and elasticmoduli does not qualitatively change the behavior of theVM [20, 41], we used Γ / ( KA ) ≈ .
289 for all simula-tions.We applied an oscillatory affine deformation of fre-quency ω either as simple shear described by the defor-mation gradient ˆ F = (cid:0) (cid:15) ( t )0 1 (cid:1) or biaxial deformation de-scribed by ˆ F = (cid:0) (cid:15) ( t ) 00 1+ (cid:15) ( t ) (cid:1) , where (cid:15) ( t ) = (cid:15) sin( ω t ).In all simulations, we used a small magnitude of defor-mation, i.e., (cid:15) = 10 − , so that we probed the linearresponse and the measured moduli were independent ofthe magnitude of the deformation. In every time step,the system evolved according to the overdamped dynam-ics in Eq. (2) after the affine deformation was applied.Equations of motion were integrated using the first orderEuler method with the time step, ∆ t ≈ . γ/ ( KA ).Measurements of the response stresses for each cell andthe entire system were taken 25 times within each cy-cle of oscillatory deformation by using Eq. (S17) in theSupporting Information, Sec. I.To ensure that we were probing the steady state, weperformed the following analysis. For example, in thecase of shear deformation, the shear stress signal τ ( t ) wasdivided into blocks of length T = 3 T , each containing3 cycles of the time period T = 2 π/ω of the drivingshear deformation. Within each block n , we performedthe Fourier transform of τ ( t ) and obtained ˜ τ n ( ω ) as˜ τ n ( ω ) = 1 T (cid:90) nT ( n − T τ ( t ) e iωt dt, (10)where n is a positive integer. Similar Fourier transformanalysis was performed for the strain, (cid:15) ( t ). The length ofthe simulation was chosen such that it contained a suffi-cient number of blocks in order for the ˜ τ n ( ω ) to reach asteady state value ˜ τ ( ω ). The obtained steady state valueof ˜ τ ( ω ) was used to calculate the storage and loss shearmoduli as described in the main text. Analogous proce-dure was applied to the normal stress, σ ( t ), in the case ofthe bulk deformation. Please refer to the SupplementaryInformation, Sec. III for a representative example of thesteady state analysis.In the solid phase, we performed rheological tests on ahexagonal tiling. In the fluid phase, however, the hexag-onal tiling is unstable. Instead, the hexagonal tiling wasrandomly perturbed and then relaxed to a nearby localstable state using the FIRE optimizer [50]. The localenergy minimum was determined with the relative ac-curacy of 10 − . A random perturbation was appliedto each vertex i , i.e., each vertex was displaced fromits original position in the hexagonal tiling by a vector δ r i = δx i e x + δy i e y , where x i and y i were Gaussian ran- dom variables with zero mean and standard deviation1 . × − √ A . The rheology of a local stable state wasthen probed following the same procedure as in the solidcase.During the energy minimization and oscillatory defor-mations, T1 transitions were allowed but were not com-mon. T1 transitions were implemented following the pro-cedure introduced by Spencer, et al. [63]. ACKNOWLEDGEMENTS
This research was primarily supported by NSF throughthe Princeton University’s Materials Research Scienceand Engineering Center DMR-2011750 and by theProject X Innovation Research Grant from the Prince-ton School of Engineering and Applied Science. RSacknowledges support by the UK BBSRC (AwardBB/N009789/1). This project was initiated during theKITP program “Symmetry, Thermodynamics and Topol-ogy in Active Matter” (ACTIVE20), and it is supportedin part by the National Science Foundation under GrantNo. NSF PHY-1748958. We would like to acknowledgeuseful discussions with Mikko Haataja.
DATA AVAILABILITY
The data supporting the results and findings of thisstudy is available from the corresponding authors uponreasonable request.
CODE AVAILABILITY
Simulation and analysis codes used in this study areavailable from the corresponding authors upon reason-able request.
CONTRIBUTIONS
ST, RS, and AK conceived the study and designed theproject. ST and NS performed numerical simulationsand analysed the data. RS developed the Vertex Modelcode used in numerical simulations. ST, NS, RS, and AKwrote the paper.
COMPETING INTERESTS
The authors declare no competing interests. [1] T. Lecuit, P.-F. Lenne, and E. Munro, Force generation,transmission, and integration during cell and tissue mor-phogenesis, Annual Review of Cell and DevelopmentalBiology , 157 (2011). [2] C.-P. Heisenberg and Y. Bella¨ıche, Forces in tissue mor-phogenesis and patterning, Cell , 948 (2013).[3] E. Hannezo and C.-P. Heisenberg, Mechanochemicalfeedback loops in development and disease, Cell , 12 (2019).[4] N. I. Petridou and C.-P. Heisenberg, Tissue rheology inembryonic organization, The EMBO Journal , e102497(2019).[5] M. H. Ross and W. Pawlina, Histology (LippincottWilliams & Wilkins, 2006).[6] D. Krndija, F. El Marjou, B. Guirao, S. Richon, O. Leroy,Y. Bellaiche, E. Hannezo, and D. M. Vignjevic, Ac-tive cell migration is critical for steady-state epithelialturnover in the gut, Science , 705 (2019).[7] L. Wolpert, C. Tickle, and A. M. Arias,
Principles ofdevelopment (Oxford University Press, USA, 2015).[8] R. A. Weinberg,
The biology of cancer (Garland Science,2013).[9] G. Forgacs, R. A. Foty, Y. Shafrir, and M. S. Steinberg,Viscoelastic properties of living embryonic tissues: aquantitative study, Biophysical Journal , 2227 (1998).[10] M. Marchetti, J. Joanny, S. Ramaswamy, T. Liverpool,J. Prost, M. Rao, and R. A. Simha, Hydrodynamics ofsoft active matter, Reviews of Modern Physics , 1143(2013).[11] P. Friedl and D. Gilmour, Collective cell migration inmorphogenesis, regeneration and cancer, Nature ReviewsMolecular Cell Biology , 445 (2009).[12] R. Alert and X. Trepat, Physical models of collective cellmigration, Annual Review of Condensed Matter Physics , 77 (2020).[13] M. Poujade, E. Grasland-Mongrain, A. Hertzog,J. Jouanneau, P. Chavrier, B. Ladoux, A. Buguin, andP. Silberzan, Collective migration of an epithelial mono-layer in response to a model wound, Proceedings of theNational Academy of Sciences , 15988 (2007).[14] X. Trepat, M. R. Wasserman, T. E. Angelini, E. Millet,D. A. Weitz, J. P. Butler, and J. J. Fredberg, Physicalforces during collective cell migration, Nature Physics ,426 (2009).[15] D. T. Tambe, C. C. Hardin, T. E. Angelini, K. Rajen-dran, C. Y. Park, X. Serra-Picamal, E. H. Zhou, M. H.Zaman, J. P. Butler, D. A. Weitz, et al. , Collective cellguidance by cooperative intercellular forces, Nature Ma-terials , 469 (2011).[16] A. Brugu´es, E. Anon, V. Conte, J. H. Veldhuis,M. Gupta, J. Colombelli, J. J. Mu˜noz, G. W. Brod-land, B. Ladoux, and X. Trepat, Forces driving epithelialwound healing, Nature Physics , 683 (2014).[17] R. Etournay, M. Popovi´c, M. Merkel, A. Nandi,C. Blasse, B. Aigouy, H. Brandl, G. Myers, G. Salbreux,F. J¨ulicher, et al. , Interplay of cell dynamics and epithe-lial tension during morphogenesis of the drosophila pupalwing, eLife , e07090 (2015).[18] T. E. Angelini, E. Hannezo, X. Trepat, M. Marquez, J. J.Fredberg, and D. A. Weitz, Glass-like dynamics of collec-tive cell migration, Proceedings of the National Academyof Sciences , 4714 (2011).[19] J.-A. Park, J. H. Kim, D. Bi, J. A. Mitchel, N. T.Qazvini, K. Tantisira, C. Y. Park, M. McGill, S.-H.Kim, B. Gweon, et al. , Unjamming and cell shape in theasthmatic airway epithelium, Nature Materials , 1040(2015).[20] D. Bi, J. Lopez, J. M. Schwarz, and M. L. Manning, Adensity-independent rigidity transition in biological tis-sues, Nature Physics , 1074 (2015).[21] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning,Motility-driven glass and jamming transitions in biolog- ical tissues, Physical Review X , 021011 (2016).[22] L. Atia, D. Bi, Y. Sharma, J. A. Mitchel, B. Gweon, S. A.Koehler, S. J. DeCamp, B. Lan, J. H. Kim, R. Hirsch, et al. , Geometric constraints during epithelial jamming,Nature Physics , 613 (2018).[23] D. M. Sussman, M. Paoluzzi, M. C. Marchetti, and M. L.Manning, Anomalous glassy dynamics in simple mod-els of dense biological tissue, EPL (Europhysics Letters) , 36001 (2018).[24] M. Czajkowski, D. M. Sussman, M. C. Marchetti, andM. L. Manning, Glassy dynamics in models of confluenttissue with mitosis and apoptosis, Soft Matter , 9133(2019).[25] S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknep-nek, and E. Bertin, Dense active matter model of motionpatterns in confluent cell monolayers, Nature Communi-cations , 1 (2020).[26] B. Szabo, G. Sz¨oll¨osi, B. G¨onci, Z. Jur´anyi, D. Selmeczi,and T. Vicsek, Phase transition in the collective migra-tion of tissue cells: experiment and model, Physical Re-view E , 061908 (2006).[27] M. Sadati, N. T. Qazvini, R. Krishnan, C. Y. Park, andJ. J. Fredberg, Collective migration and cell jamming,Differentiation , 121 (2013).[28] S. Garcia, E. Hannezo, J. Elgeti, J.-F. Joanny, P. Sil-berzan, and N. S. Gov, Physics of active jamming duringcollective cellular motion in a monolayer, Proceedings ofthe National Academy of Sciences , 15314 (2015).[29] M. Merkel and M. L. Manning, A geometrically con-trolled rigidity transition in a model for confluent 3d tis-sues, New Journal of Physics , 022002 (2018).[30] J. A. Mitchel, A. Das, M. J. O’Sullivan, I. T. Stancil,S. J. DeCamp, S. Koehler, O. H. Oca˜na, J. P. Butler,J. J. Fredberg, M. A. Nieto, et al. , In primary airwayepithelial cells, the unjamming transition is distinct fromthe epithelial-to-mesenchymal transition, Nature Com-munications , 1 (2020).[31] B. B´enaz´eraf, P. Francois, R. E. Baker, N. Denans, C. D.Little, and O. Pourqui´e, A random cell motility gradi-ent downstream of fgf controls elongation of an amnioteembryo, Nature , 248 (2010).[32] A. K. Lawton, A. Nandi, M. J. Stulberg, N. Dray, M. W.Sneddon, W. Pontius, T. Emonet, and S. A. Holley, Reg-ulated tissue fluidity steers zebrafish body elongation,Development , 573 (2013).[33] A. Mongera, P. Rowghanian, H. J. Gustafson, E. Shelton,D. A. Kealhofer, E. K. Carn, F. Serwane, A. A. Lucio,J. Giammona, and O. Camp`as, A fluid-to-solid jammingtransition underlies vertebrate body axis elongation, Na-ture , 401 (2018).[34] N. Desprat, A. Guiroy, and A. Asnacios, Microplates-based rheometer for a single living cell, Review of Scien-tific Instruments , 055111 (2006).[35] G. Salbreux, G. Charras, and E. Paluch, Actin cortexmechanics and cellular morphogenesis, Trends in Cell Bi-ology , 536 (2012).[36] H. Berthoumieux, J.-L. Maˆıtre, C.-P. Heisenberg, E. K.Paluch, F. J¨ulicher, and G. Salbreux, Active elastic thinshell theory for cellular deformations, New Journal ofPhysics , 065005 (2014).[37] A. R. Harris, L. Peter, J. Bellis, B. Baum, A. J. Kabla,and G. T. Charras, Characterizing the mechanics ofcultured cell monolayers, Proceedings of the NationalAcademy of Sciences , 16449 (2012). [38] A. Bonfanti, J. Fouchard, N. Khalilgharibi, G. Char-ras, and A. Kabla, A unified rheological model for cellsand cellularised materials, Royal Society Open Science ,190920 (2020).[39] D. Matoz-Fernandez, E. Agoritsas, J.-L. Barrat,E. Bertin, and K. Martens, Nonlinear rheology in a modelbiological tissue, Physical Review Letters , 158105(2017).[40] T. Nagai and H. Honda, A dynamic cell model for theformation of epithelial tissues, Philosophical Magazine B , 699 (2001).[41] R. Farhadifar, J.-C. R¨oper, B. Aigouy, S. Eaton, andF. J¨ulicher, The influence of cell mechanics, cell-cell inter-actions, and proliferation on epithelial packing, CurrentBiology , 2095 (2007).[42] A. G. Fletcher, M. Osterfield, R. E. Baker, and S. Y.Shvartsman, Vertex models of epithelial morphogenesis,Biophysical Journal , 2291 (2014).[43] D. L. Barton, S. Henkes, C. J. Weijer, and R. Sknepnek,Active vertex model for cell-resolution description of ep-ithelial tissue mechanics, PLoS Computational Biology , e1005569 (2017).[44] M. Moshe, M. J. Bowick, and M. C. Marchetti, Geometricfrustration and solid-solid transitions in model 2d tissue,Physical Review Letters , 268105 (2018).[45] M. Popovi´c, V. Druelle, N. Dye, F. J¨ulicher, andM. Wyart, Inferring the flow properties of epithelial tis-sues from their geometry, New Journal of Physics xx ,accepted manuscript (2020).[46] I. Tah, T. Sharp, A. Liu, and D. M. Sussman, Quantifyingthe link between local structure and cellular rearrange-ments using information in models of biological tissues,Soft Matter (2021).[47] D. Staple, R. Farhadifar, J.-C. R¨oper, B. Aigouy,S. Eaton, and F. J¨ulicher, Mechanics and remodelling ofcell packings in epithelia, The European Physical JournalE , 117 (2010).[48] D. Bi, J. H. Lopez, J. Schwarz, and M. L. Manning, En-ergy barriers and cell migration in densely packed tissues,Soft Matter , 1885 (2014).[49] L. Yan and D. Bi, Multicellular rosettes drive fluid-solidtransition in epithelial tissues, Physical Review X ,011029 (2019).[50] E. Bitzek, P. Koskinen, F. G¨ahler, M. Moseler, andP. Gumbsch, Structural relaxation made simple, Phys-ical Review Letters , 170201 (2006).[51] E. M. Purcell, Life at low reynolds number, Americanjournal of physics , 3 (1977).[52] D. Frenkel and B. Smit, Understanding molecular simu-lation: from algorithms to applications , Vol. 1 (Elsevier,2001).[53] S. Curran, C. Strandkvist, J. Bathmann, M. de Gennes,A. Kabla, G. Salbreux, and B. Baum, Myosin ii controlsjunction fluctuations to guide epithelial tissue ordering,Developmental Cell , 480 (2017).[54] A. Nestor-Bergmann, G. Goddard, S. Woolner, and O. E.Jensen, Relating cell shape and mechanical stress ina spatially disordered epithelium using a vertex-basedmodel, Mathematical Medicine and Biology: A Journalof the IMA , i1 (2018).[55] R. G. Larson, The structure and rheology of complex flu-ids , Vol. 150 (Oxford University Press New York, 1999).[56] N. Murisic, V. Hakim, I. G. Kevrekidis, S. Y. Shvarts-man, and B. Audoly, From discrete to continuum mod- els of three-dimensional deformations in epithelial sheets,Biophysical Journal , 154 (2015).[57] C. Kane and T. Lubensky, Topological boundary modesin isostatic lattices, Nature Physics , 39 (2014).[58] T. Lubensky, C. Kane, X. Mao, A. Souslov, and K. Sun,Phonons and elasticity in critically coordinated lattices,Reports on Progress in Physics , 073901 (2015).[59] J. Paulose, B. G.-g. Chen, and V. Vitelli, Topologicalmodes bound to dislocations in mechanical metamateri-als, Nature Physics , 153 (2015).[60] S. D. Huber, Topological mechanics, Nature Physics ,621 (2016).[61] D. Z. Rocklin, S. Zhou, K. Sun, and X. Mao, Trans-formable topological mechanical metamaterials, NatureCommunications , 1 (2017).[62] X. Mao and T. C. Lubensky, Maxwell lattices and topo-logical mechanics, Annual Review of Condensed MatterPhysics , 413 (2018).[63] M. A. Spencer, Z. Jabeen, and D. K. Lubensky, Vertexstability and topological transitions in vertex models offoams and epithelia, The European Physical Journal E , 2 (2017). Linear Viscoelastic Properties of the Vertex Model for Epithelial TissuesSupplementary Information
I. STRESS TENSOR FOR EACH CELL IN THE VERTEX MODEL
The expression for the stress tensor of the vertex model (VM) follows for the most part the derivation in the Ref. [54]in the main text, but there is one important difference as discussed below. For completeness, we here reproduce thekey steps of the calculation. The energy functional of the VM can be rewritten as E = (cid:88) C (cid:20) K A C − A ) + Γ2 P C − Λ P C (cid:21) , (S1)where we have defined Λ = Γ P and we have omitted the constant term Γ ( P ) compared to the Eq. (1) in the maintext. Other variables have the same meaning as in the main text. Force on the vertex i is then given as F i = −∇ r i E = F area i + F perim i , (S2)where F area i and F perim i are the contributions to the force due to area and perimeter terms, respectively. If we introducepressure p C = − K ( A C − A ) , (S3)for a cell C then it is straightforward to show that F area i = (cid:88) e (cid:0) p C e,r − p C e,l (cid:1) e z × l e ≡ (cid:88) C f area i,C , (S4)and F perim i = (cid:88) e (cid:0) Γ (cid:0) P C e,l + P C e,r (cid:1) − (cid:1) ˆ l e ≡ (cid:88) C f perim i,C . (S5)In these expressions, the e − sum loops over the edges originating at vertex i in the clockwise direction, index C e,l ( C e,r ) corresponds to the cell being to the left (right) of edge e when looking away from the vertex i , and l e , is a vectorpointing along the edge e , whose magnitude is equal to the length of the edge (Fig. S1a). e z is the unit-length vectorperpendicular to the plane of the cell sheet that is assumed to be in the xy plane. Finally, hat ˆ · denotes unit-lengthvectors. The summations over edges in Eqs. (S4) and (S5) couple neighboring cells that share a junction, and theycan be rewritten to equations in terms of contributions of individual cells to a force on a given vertex. Using the FIG. S1. (a) Blue arrow denotes an edge vector l e shared by cells C e,l and C e,r . The arc arrow marks the direction of in whichterms are summed in Eqs. (S4) and (S5). (b) Vertices in a cell are ordered in the counterclockwise direction and red arrowsshow the definitions for edge vectors between the neighboring vertices. ∇ r i A C = e z × ( l i − ,i − l i +1 ,i ), where l i ± ,i denotes the edge vector between the vertex i and the neighboringvertex i ± C counted in the counterclockwise directions (Fig. S1b), gives f area i,C = 12 δ C,i p C e z × ( l i − ,i − l i +1 ,i ) , (S6)where δ C,i = 1 if the vertex i belongs to cell C and δ C,i = 0 otherwise. Similarly, we find the perimeter component ofthe force of cell C on the vertex i as f perim i,C = δ C,i (Γ P C − Λ) (cid:16) ˆl i − ,i + ˆ l i +1 ,i (cid:17) . (S7)Further, using the assumption that dissipative force on vertex i is equally distributed on neighboring cells, we havethat the dissipative force on vertex i with coordination number z i due to cell C is f dissip i,C = − γ v i /z i . In the overdampedlimit, the total mechanical force on each cell is balanced by dissipative forces, so for vertex i , the total mechanicalforce f i = F area i + F perim i is balanced by the dissipative forces (see Fig. S2), i.e. f i − z i γ v i (cid:88) C δ C,i = f i − γ v i = 0 , (S8)where v i is the velocity of vertex i . However, this is not the case for components of vertex forces due to individualcells. In order to construct the stress tensor defined for a cell, we need to change perspective and account for all theforces acting on the cell from the surroundings. The resultant force F C,i on the cell C due to the vertex i is F C,i = f C,i − z i γ v i δ C,i = − f i,C − z i γ v i δ C,i , (S9)where we used the 3rd Newton’s law to relate the elastic part of the force f C,i from the vertex i acting on cell C tothe equal and opposite force f i,C = f area i,C + f perim i,C from the cell C acting on the vertex i . This change of sign is thedifference compared to the Ref. [54]. Note that the viscous drag force − z i γ v i δ C,i acts at the external force on thecell and we thus keep the same sign. f i,C'' C C'C'' f i,C f i,C' i f dissip FIG. S2. At vertex i , the total mechanical force, f i,C + f i,C (cid:48) + f i,C (cid:48)(cid:48) , in terms of the sum of the contributions from individualcells, is balanced by the dissipative force, f dissip , from the external environment. Z C vertices as r C = Z C (cid:80) i ∈ C r i , its velocityas v C = ˙ r C = Z C (cid:80) i ∈ C v i , and the position of vertex i relative to the cell center C as ˜ r i = r i − r C . Thus we cancompute the total force on the cell C as F C = Z C (cid:88) i =1 F C,i = Z C (cid:88) i =1 (cid:18) − f i,C − γ z i v i (cid:19) = f C − γ Z c z i v C = 0 , (S10)which is equal to zero due to the overdamped limit.We can now construct the stress tensor. For a tensor ˆ σ that is symmetric and divergence-free, defined over an area A with perimeter ∂A , we have ˆ σ = ∇ · ( R ⊗ ˆ σ ), where R is an arbitrary position vector. Thus, we write (cid:90) A ˆ σ dA = (cid:90) A ∇ · ( R ⊗ ˆ σ ) dA = (cid:90) ∂A R ⊗ ˆ σ · ˆ n dl, (S11)where we used the Stokes theorem to convert the integral over area A into the integral over the boundary ∂A with theoutwards pointing unit normal vector ˆ n and the length element dl . If we note that ˆ σ · n is the force density acting onthe boundary and assume that stress is constant over the entire area of the cell, then the Eq. (S11) can be rewrittenas A C ˆ σ C = (cid:88) i ∈ C r i ⊗ F C,i = (cid:88) i ∈ C ( r C + ˜ r i ) ⊗ F C,i = r C ⊗ (cid:88) i ∈ C F C,i + (cid:88) i ∈ C ˜ r i ⊗ F C,i A C ˆ σ C = r C ⊗ F C + (cid:88) i ∈ C ˜ r i ⊗ (cid:20) − f i,C − γ z i v i (cid:21) = (cid:88) i ∈ C ˜ r i ⊗ (cid:20) − f i,C − γ z i v i (cid:21) . (S12)Therefore the stress tensor ˆ σ C for cell C isˆ σ C = 1 A C (cid:88) i ∈ C ˜ r i ⊗ (cid:20) − f i,C − γ z i v i (cid:21) , ˆ σ C = − A C (cid:88) i ∈ C ˜ r i ⊗ (cid:18) p C e z × ( l i − ,i − l i +1 ,i ) + (Γ P C − Λ) (cid:16) ˆl i − ,i + ˆ l i +1 ,i (cid:17) − γ z i v i (cid:19) , ˆ σ C = − p C A C (cid:88) i ∈ C ˜ r i ⊗ [ e z × ( l i − ,i − l i +1 ,i )] − (Γ P C − Λ) A C (cid:88) i ∈ C ˜ r i ⊗ (cid:16) ˆl i − ,i + ˆ l i +1 ,i (cid:17) − γA C (cid:88) i ∈ C z i ˜ r i ⊗ v i . (S13)Now, let us compute the tensor ˜ r i ⊗ [ e z × ( l i − ,i − l i +1 ,i )]. We first note that l i ± ,i = ˜ r i ± − ˜ r i . Further,[˜ r i ⊗ ( e z × ( l i − ,i − l i +1 ,i ))] αβ = ˜ r i,α [ e z × ( ˜r i − − ˜r i +1 )] β = ˜ r i,α ( δ β,y (˜ r i − ,x − ˜ r i +1 ,x ) − δ β,x (˜ r i − ,y − ˜ r i +1 ,y )) . (S14)If we now multiply the last expression by and sum over all vertices, we obtain (cid:34) (cid:88) i ∈ C ˜ r i ⊗ ( e z × ( l i − ,i − l i +1 ,i )) (cid:35) αβ = A C δ αβ . (S15)We proceed to find ˜ r i ⊗ (cid:16) ˆl i − ,i + ˆ l i +1 ,i (cid:17) . Since all these terms appear under a sum that is cyclic (i.e., over all vertices)we can shift indices as convenient to write, (cid:88) i ∈ C ˜ r i ⊗ (cid:16) ˆl i − ,i + ˆ l i +1 ,i (cid:17) = (cid:88) i ∈ C (cid:20) ˜ r i +1 ⊗ (˜ r i − ˜ r i +1 ) | ˜ r i − ˜ r i +1 | + ˜ r i ⊗ (˜ r i +1 − ˜ r i ) | ˜ r i +1 − ˜ r i | (cid:21) = (cid:88) i ∈ C ( − ˜ r i +1 + ˜ r i ) ⊗ (˜ r i +1 − ˜ r i ) | ˜ r i +1 − ˜ r i | = − (cid:88) i ∈ C | l i +1 ,i | (cid:104) ˆ l i +1 ,i ⊗ ˆ l i +1 ,i (cid:105) . (S16)This enables us to write the stress tensor for cell C asˆ σ C = − p C ˆ I + (Γ P C − Λ) A C (cid:88) i ∈ C | l i +1 ,i | ˆ l i +1 ,i ⊗ ˆ l i +1 ,i − γ A c (cid:88) i ∈ C z i (˜ r i ⊗ v i + v i ⊗ ˜ r i ) , (S17)5where we symmetrized the last viscous part. When calculating stresses according the the above Eq. (S17) we removedthe affine part of the velocity field by replacing the vertex velocities v i with velocities v i − v affine i , where the affinepart of the velocity field is v affine i = d ˆ F dt r i . (S18)Here, ˆ F is the deformation gradient due to the applied external deformation. For the oscillatory external deformationwith the frequency ω , the stress contribution due to the friction resulting from the affine part of the velocity field v affine i is linearly proportional to the driving frequency ω and it dominates at large frequencies (see Fig. S3). FIG. S3. The effect of the friction resulting from the affine velocity v affine i in Eq. (S18) on the loss shear modulus for p = 3 . II. CONNECTION BETWEEN STRESS RESPONSE AND RHEOLOGY
Fig. S4 shows a typical average shear stress τ ( t ) in response to an applied oscillatory simple shear with thestrain (cid:15) = (cid:15) sin( ω t ). The shear stress response can be represented as τ ( t ) = τ sin( ω t + δ ) = τ cos( δ ) sin( ω t ) + τ sin( δ ) cos( ω t ). The storage shear modulus is related to the in-phase response and is defined as G (cid:48) = ( τ /(cid:15) ) cos δ .The loss shear modulus is related to the out-of-phase response and is defined as G (cid:48)(cid:48) = ( τ /(cid:15) ) sin δ [55]. FIG. S4. Typical shear stress (red curve) as a function of time in response to a periodic shear strain (blue curve) in (a) thesolid phase and (b) the fluid phase. The shear stress is averaged over all cells.
III. APPROACH OF THE RESPONSE STRESS TOWARDS THE STEADY STATE
Here, we show an example of how the steady state shear stress ˜ τ ( ω ) is measured in response to the appliedoscillatory simple shear with a time period T = 27 . γ/ ( KA ) = 2 π/ω for the shape parameter p = 3 . p c ≈ .
722 for the solid-fluid transition. The shear stress signal τ ( t ) was divided intoblocks of length T = 3 T , each containing 3 cycles of the time period of the driving shear deformation (see Fig. S5a).Within each block n , we performed the Fourier transform of τ ( t ) and obtained ˜ τ n ( ω ) as˜ τ n ( ω ) = 1 T (cid:90) nT ( n − T τ ( t ) e iωt dt, (S19)where n is a positive integer. The value of ˜ τ n ( ω ) converges exponentially to the steady state value (see Fig. S5b),where the relaxation time is related to the characteristic time scales of the viscoelastic models (see Fig. 2c,d in themain text). For values of p far away from p c , the system quickly reaches a steady state (within 3–6 cycles). As p approaches p c the relaxation times become much longer, which is reflecting the diverging characteristic time scales ofthe viscoelastic models (see Fig. 2c,d in the main text). FIG. S5. Approach of the response shear stress towards the steady state. (a) The shear stress signal τ ( t ) was divided intoblocks indicated by the vertical dashed lines. (b) Fourier transform of the response shear stress ˜ τ n ( ω ) at the driving frequency, ω , as a function of the block number, n . IV. RESIDUAL HYDROSTATIC STRESS IN THE SOLID PHASE
In the solid phase, we studied the rheology of the hexagonal lattice with each cell of area A C = A but with theperimeter P C not equal to the preferred perimeter P , which induces residual hydrostatic stress in equilibrium. Thisresidual stress can be eliminated if the lattice is uniformly rescaled by a factor α , which minimizes the followingdimensionless energy per cell, e C ( α ) = 12 (cid:0) α − (cid:1) + ˜Γ2 ( αp C − p ) , (S20)where e C = E C KA , ˜Γ = Γ KA , p C = P C √ A = √ ≈ . α is a root of e (cid:48) C ( α ) = 0. In the solid phase, α <
1, andthe system shrinks to relax the residual stress if it is not held by the periodic boundary conditions. At the solid-fluidtransition point, α = 1 since the area and perimeter of each cell match their preferred values simultaneously. If theresidual stress is eliminated by rescaling the box, the rheology of the system subject to a simple shear can still bedescribed by the SLS model, although the fitted values of spring constants are different, as shown in Fig. S6. FIG. S6. The fitted spring constants in the solid phase when the simulation box is not rescaled (closed symbols) and rescaled(open symbols) to eliminate residual stresses. V. COLLAPSE OF STORAGE AND LOSS SHEAR MODULI IN THE FLUID PHASE
In Fig. 1f in the main text, we showed the collapse of storage and loss shear moduli for the fluid phase in the lowfrequency regime. Here we show the collapse in the high frequency range (see Fig. S7), where we took into accountthat the relevant characteristic timescale scales as η /E ∼ γ/ ( KA ). FIG. S7. The collapse of the storage ( G (cid:48) ) and loss ( G (cid:48)(cid:48) ) shear moduli curves in the high frequency regime for different valuesof p for the fluid phase. VI. EFFECTS OF THE INITIAL PERTURBATION IN THE FLUID PHASE
We note that the rheological behavior in the fluid phase is sensitive to the magnitude σ D of the initial perturbationthat was used to obtain different local energy minima configurations. In the main text, we showed the fitted valuesof spring and dashpot constants (Fig. 2) for the local energy minima configurations that were obtained by displacingeach vertex coordinate of the hexagonal tiling by a Gaussian random variable with zero mean and standard deviation σ D = 1 . × − √ A . Here, we show that the fitted values of the spring and dashpot constants are somewhat sensitiveto the magnitude σ D of the random perturbation (see Fig. S8). FIG. S8. Fitted values of (a) spring and (b) dashpot constants for the system under simple shear deformation as a function ofthe target cell-shape parameter, p , and the magnitude σ D of the random perturbation that was used to obtain different localenergy minima configurations in the fluid phase. Errorbars correspond to the standard deviation for simulations that wererepeated for configurations that correspond to different local energy minima. VII. TUNING PHASE TRANSITION WITH DIFFERENT MODES OF DEFORMATION.
In the main text, we showed that the solid-fluid transition point can be tuned by uniaxially pre-compressing/stretching the system. Here we discuss other deformation modes that can also tune the transition.If the hexagonal lattice is deformed biaxially by ˆ F = (cid:0) a a (cid:1) , the shear modulus due to the affine deformation becomes G affine = 3 √ (cid:18) − p ap c (cid:19) . (S21)By setting G affine to 0, the phase boundary in the a − p plane is a ( p ) = p p c . (S22)Similarly, consider a pure shear deformation ˆ F = (cid:0) a
00 1 /a (cid:1) . The shear modulus due to the affine deformation thenbecomes G affine = 2 √ (cid:0) √ a + 3 a √ a (cid:1) (cid:0) √ / (1 + √ a ) − ap (cid:1) Γ3 / a (1 + 3 a ) / , (S23)and the phase boundary is p ( a ) = 2 √ / (1 + √ a )3 a . (S24)The phase diagrams for a system that is under biaxial deformation or pure shear are shown in Fig. S9. The phaseboundary in a − p plane follows Eq. (S22) for biaxial deformation and Eq. (S24) for pure shear. The system can berigidified by stretching or shearing. FIG. S9. Phase diagrams when the system is under (a) biaxial deformation and (b) pure shear.
Note that there are two equivalent ways of derive the shear modulus due to affine deformation. The first one isto calculate the energy density of the system perturbed by an additional simple shear ˆ F = (cid:0) (cid:15) (cid:1) where (cid:15) (cid:28)
1. Forexample, for the hexagonal tiling without any pre-deformation (i.e., regular hexagons), the energy density can beexpanded in a power series in (cid:15) as EA = 12 3 √ (cid:18) − p p c (cid:19) Γ (cid:15) + o (cid:0) (cid:15) (cid:1) ≡ G affine (cid:15) + o (cid:0) (cid:15) (cid:1) , (S25)where p c = 2 √ / ≈ .
722 and we neglected the constant term. The quadratic term characterizes the linearresponse of the system, which gives the shear modulus as in Eq. (5) in the main text. The second approach is todirectly use the elastic part of the stress formula Eq. (S17) by setting the velocity v i to be 0. Assume the system isperturbed by a simple shear ˆ F = (cid:0) (cid:15) (cid:1) where (cid:15) (cid:28)
1, and calculate the shear stress. The coefficient of the leadingorder term in (cid:15) is the shear modulus, which coincides with the modulus from the energy calculation. Similar derivationof the shear modulus can be carried out for the pre-deformed hexagonal tilings.0
VIII. SPECTRUM OF THE NORMAL MODES
We calculated the eigenvalues λ of the Hessian matrix associated with the energy functional of the VM. We associateeach positive eigenvalue λ with a corresponding eigenfrequency ω = √ λ , which describes the oscillations of that modeas the system is perturbed about its stable point. Fig. S10 shows the cumulative density of states, which is defined as N ( ω ) = (cid:90) ∞ + D ( ω (cid:48) ) dω (cid:48) + N ( λ = 0) θ ( ω ) , (S26)where D ( ω ) is density of states, N ( λ = 0) is the fraction of zero eigenvalues and θ ( ω ) is the Heaviside step function.In the solid phase, there are no zero modes other than the two translational rigid body motions. In the fluid phase,however, approximately half of the eigenmodes are zero modes. As p approaches the critical value p c in both solidand fluid phase, N ( ω ) curves move to the left so the system becomes softer, which is consistent with the dependenceof the spring constants on p shown in Fig. 2a in the main text.shown in Fig. 2a in the main text.