Nonlinear and Nonlocal Elasticity in Coarse-Grained Differential-Tension Models of Epithelia
NNonlinear and Nonlocal Elasticity in Coarse-Grained Differential-Tension Models of Epithelia
Pierre A. Haas ∗ and Raymond E. Goldstein † Department of Applied Mathematics and Theoretical Physics, Centre for Mathematical Sciences,University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom (Dated: October 23, 2018)The shapes of epithelial tissues result from a complex interplay of contractile forces in the cytoskeleta of thecells in the tissue, and adhesion forces between them. A host of discrete, cell-based models describe these forcesby assigning different surface tensions to the apical, basal, and lateral sides of the cells. These differential-tensionmodels have been used to describe the deformations of epithelia in different living systems, but the underlyingcontinuum mechanics at the scale of the epithelium are still unclear. Here, we derive a continuum theory fora simple differential-tension model of a two-dimensional epithelium and study the buckling of this epitheliumunder imposed compression. The analysis reveals how the cell-level properties encoded in the differential-tensionmodel lead to linear, nonlinear as well as nonlocal elastic behavior at the continuum level.
I. INTRODUCTION
Intercellular adhesion proteins and cortical actin networksare well established as regulators of cell surface mechanics,and hence of the deformations of epithelia during morpho-genesis [1]. Ever since the seminal work of Odell et al. [2],these cellular components have therefore underlain mathemat-ical models of epithelia [3, 4]. One large class of such modelsare differential-tension models [5–12], in which cell polarity,cell-cell adhesion properties and the actomyosin network in-duce different surface tensions in different sides of the discretecells. Coupling the mechanical models describing epithelialdeformations to models of the intracellular biochemistry isa key challenge in the field [13], but some progress has re-cently been made by coupling the differential-tension modelof Ref. [8] to the diffusion of a ‘mechanogen’ that inducescontractility [14].While the differential-tension models can quantitatively re-produce the morphology of epithelial folds in many differentliving systems [11], it is likely that, in general multilayeredepithelia, the formation of epithelial folds must be ascribed toa combination of these intra-epithelial stresses and differentialgrowth of different parts of the tissue. Models based on the lat-ter only have for example been invoked to describe, at the scaleof the epithelium, the formation of cortical convolutions in thebrain [15–20] and of the intestinal villi [21–23], the ‘fication’of which lends itself to the pun that gave Ref. [23] its title. Thecoarse-grained limit of the differential-tension models at thisscale is less well-studied, however, and this continuum limit isthe topic of this paper.We shall focus on the most basic setup of these differential-tension models [5, 7, 10] in an epithelial monolayer: the apical,basal, and lateral sides of the cells have respective areas A a , A b , and A (cid:96) . The internal state of the cells induces differentsurface tensions Γ a , Γ b , and Γ (cid:96) in the apical, basal, and lateralsides of the cells, respectively. The energy of a single celltherefore reads E = Γ a A a + Γ b A b + Γ (cid:96) A (cid:96) , (1)where the factor of 1 / ∗ [email protected] † [email protected] additional physics such as a basement membrane [11] or aconfining vitelline membrane [6]. In a full three-dimensionalsetup, this leads to the study of the shapes of prism-shapedcells [8, 12]. Here, we shall restrict to the two-dimensionalsetup [7, 10] of an epithelium consisting of isosceles trape-zoidal cells of parallel apical and basal sides of lengths L a and L b . These trapezoids are joined up across their lateral sides,which have equal length L (cid:96) (Fig. 1). Since there are two suchlateral sides, the cell energy (1) reduces to E = Γ a L a + Γ b L b + Γ (cid:96) L (cid:96) , (2)per unit extent in the third dimension [7, 10]. A different two-dimensional limit is obtained by averaging over the thicknessof the cell sheet and describing in-plane deformations only.Such models, termed area- and perimeter-elasticity models,have been studied extensively [24–28].The simplest problem in the mechanics of elastic rods istheir Euler buckling under applied forces [29]; it is thereforemeet to ask how the buckling behavior of an active materialsuch as this model epithelium differs from that of an elasticrod. This problem was considered in Ref. [10], where thecontinuum limit of Eq. (2) was mapped to Euler’s Elasticaequation [30]. These calculations, complementing simula-tions of the discrete model in Ref. [7], were phrased in termsof spontaneous buckling of the epithelium, but an additionalcompressive or extensile force is required to produce these de-formations, making the analysis of Ref. [10] more appropriateto the present setup of buckling under imposed forces. More-over, the analysis of Ref. [10] is not completely consistent withthe discrete model, since it does not impose the condition thatthe trapezoidal cells match up exactly along their lateral sides. basal sideapical side L b L a L ‘ FIG. 1. Model epithelium. Isosceles trapezoidal cells of apical andbasal bases L a and L b are connected along their lateral sides, whichhave length L (cid:96) . The continuous line and shaded area provide acartoon of the continuum limit. a r X i v : . [ c ond - m a t . s o f t ] O c t L a L ‘ L φ L b K ψ ψ ψ = ◦ − φ ◦ − φ ◦ − φ ◦ − φ φ + φ + φ φ + φ φ + φ xy (a) (b)FIG. 2. Cell Geometry. (a) Geometry of a single isosceles trapezoidal cell of mean base K and height L , and sidelengths L a , L b , L (cid:96) , the lateralsides being at an angle 2 φ to each other. (b) Definition of the tangent angle ψ of the midline below the horizontal. The geometry of contiguouscells defines the relation between φ and ψ , as expressed in Eq. (11). In this paper, we perform a consistent asymptotic expansionof the discrete geometry of this model, revealing nonlinear andnonlocal elastic terms in the continuum limit. We then analyzethe buckling behavior of the continuum model under imposedcompression analytically and numerically.
II. CONTINUUM MODELA. Single-cell Energy
To obtain the energy of a single cell, we express the side-lengths of the trapezoidal cells in terms of their mean base K ,their height L , and the angle 2 φ that their lateral sides makewith each other (Fig. 2a): L a = K + L (cid:96) sin φ, (3a) L b = K − L (cid:96) sin φ, (3b) L = L (cid:96) cos φ. (3c)Incompressibility implies the cell area conservation constraint A c = K L . (3d)Upon eliminating K and L using these relations, the energy ofa single cell is expressed, from Eq. (2) and as a function of L (cid:96) and φ , as E = ( Γ a + Γ b ) A c L (cid:96) sec φ + ( Γ a − Γ b ) L (cid:96) sin φ + Γ (cid:96) L (cid:96) . (4)
1. Non-dimensionalization
We non-dimensionalize this expression by scaling lengthswith the square root of the cell area and thus define the non-dimensional length of the lateral sides of the trapezoidal cells, λ = L (cid:96) / A / . We further set (cid:96) = L / A / and κ = K / A / . Fi-nally, following Ref. [10], we introduce the parameters α = Γ a Γ (cid:96) , β = Γ b Γ (cid:96) , (cid:96) = (cid:112) α + β, δ = α − β. (5) We note that (cid:96) is the (uniform and non-dimensionalized)thickness of the epithelium in the flat configuration. Hence s = / (cid:96) is the (non-dimensionalized) width of a singlecell (i.e. its arclength in the flat configuration). The non-dimensional energy e = E / Γ (cid:96) A / of a single cell is therefore e = (cid:96) λ sec φ + λ (cid:0) δ sin φ + (cid:1) . (6)Without loss of generality, we assume that δ > φ <
2. Transition to Constricted Cells
The transition to constricted triangular cells is a geometricsingularity in the discrete model. These triangular cells arizeas limiting cases of the trapezoidal cells when L a = L b =
0. Using Eq. (3), the conditions L a , L b (cid:62) λ (cid:54) φ | sin φ | = | sin 2 φ | . (7) B. Energy of an Epithelium
In the continuum limit, we take φ to be a function of thearclength s of the midline of the undeformed, flat epithelium.Summing over all the cells, we obtain the non-dimensionalenergy E of the epithelium, E = ∫ e ( φ ) (cid:96) d s = (cid:96) ∫ (cid:18) (cid:96) λ sec φ + λ(cid:96) (cid:0) + δ sin φ (cid:1)(cid:19) d s , (8)By imposing an energy density equal to 1 / s = (cid:96) times the(non-dimensional) energy of a single cell and integrating withrespect to the reference arclength in this manner, we haveimposed local cell area conservation [31].The boundary conditions are most naturally expressed interms of the angle ψ of the deformed midline of the epithelium below the horizontal (Fig. 2b). We therefore express the energyin terms of ψ . This is usefully done in the scaling limit (cid:96) (cid:29)
1. Asymptotic Expansion
We make two further scaling assumptions: λ = O ( (cid:96) ) , φ = O (cid:0) (cid:96) − (cid:1) . (9)These scalings correspond to the regime (cid:96) ∼ (cid:96) and φ (cid:28) λ ∼ (cid:96) ∼ (cid:96) .Further, area conservation (3d) requires that κ ∼ / (cid:96) , andthus, from Eqs. (3a,3b), we must have φ (cid:46) κ / λ ∼ / (cid:96) . Thesecond scaling thus corresponds to the largest deformationsallowed. We therefore introduce the parameter Λ = λ / (cid:96) = O ( ) . (10)We are now set up to relate φ and ψ , for which purpose we usethe geometric relation ψ ( s + ks ) − ψ ( s ) = φ ( s ) + φ ( s + s ) + · · · + φ (cid:0) s + ( k − ) s (cid:1) + φ ( s + ks ) , (11)valid for any positive integer k , as sketched in Fig. 2b. InAppendix A, we show that, with our scaling assumptions, thecontinuum limit of this relation is ψ (cid:48) ( s ) = ∞ (cid:213) m = B m ( m ) ! φ ( m ) ( s ) (cid:96) m − , (12)wherein B = , B = − , . . . are the Bernoulli numbers (ofthe first kind) [32]. The next step is to invert this series, toexpress φ in terms of the derivatives of ψ . While we arenot aware of any explicit expression for the coefficients of theinverted series, it is straightforward to invert the series order-by-order by substituting back and forth, and thus obtain φ ( s ) = ψ (cid:48) ( s ) (cid:96) − ψ (cid:48)(cid:48)(cid:48) ( s ) (cid:96) + ψ ( v ) ( s ) (cid:96) + · · · , (13)where dashes denote differentiation with respect to s . InRef. [10], only the first term of this expansion was obtained.Inclusion of the second term will enable us to analyze thebuckling behavior of the epithelium in what follows.
2. Shape Equations for the Buckled Epithelium
We describe the shape of the buckled epithelium by the coor-dinates (cid:0) x ( s ) , y ( s ) (cid:1) of the centreline of the epithelium, definedby the axes in Fig. 2b. To derive the continuum equationsdescribing the centreline, we begin by projecting the discretegeometry onto the axes, x ( s + s ) − x ( s ) = (cid:16) κ ( s ) cos ψ ( s ) + κ ( s + s ) cos ψ ( s + s ) (cid:17) , (14a) y ( s + s ) − y ( s ) = (cid:16) κ ( s ) sin ψ ( s ) + κ ( s + s ) sin ψ ( s + s ) (cid:17) . (14b) Using κ ( s ) = ( (cid:96) Λ ) − sec φ ( s ) and expanding these equationsorder-by-order in inverse powers of (cid:96) using Eq. (13), we ob-tain, after a considerable amount of algebra [33], Λ d x d s = f cos ψ − g sin ψ, Λ d y d s = f sin ψ + g cos ψ, (15)where f = + ψ (cid:48) (cid:96) + ψ (cid:48) + ψ (cid:48)(cid:48) + ψ (cid:48) ψ (cid:48)(cid:48)(cid:48) (cid:96) + O (cid:0) (cid:96) − (cid:1) , (16a) g = ψ (cid:48)(cid:48) (cid:96) + ψ (cid:48) ψ (cid:48)(cid:48) − ψ (iv) (cid:96) + O (cid:0) (cid:96) − (cid:1) . (16b)Integrating these differential equations yields the shape of thebuckled epithelium. Deviations from the ‘standard’ values f = g = O (cid:0) (cid:96) − (cid:1) .
3. Derivation of the Governing Equation
We shall seek to describe buckled configurations of an ep-ithelium of undeformed length 2 Σ (cid:29) s . We change variablesby introducing σ = s / Σ , use dots to denote differentiationwith respect to σ , and define Ξ = (cid:96) Σ (cid:38) O ( ) . (17)Since Ξ = Σ / s , the number of cells in the epithelium issimply N = Ξ .We shall seek buckled solutions with clamped boundaryconditions and a prescribed relative compression D , so that x ( ) − x ( ) = ( − D ) , y ( ) = y ( ) , (18)where the coordinates are now expressed relative to the scaledarclength σ . We shall restrict to symmetrically buckled con-figurations for which ψ ( σ ) = − ψ ( − σ ) . The second conditionabove is then satisfied. We may further reduce the solution tothe range 0 (cid:54) σ (cid:54)
1, with the condition of prescribed com-pression reading x ( ) − x ( ) = − D . To minimize the energyof the epithelium at this imposed displacement, we thereforeconsider the Lagrangian L = ∫ (cid:18) sec φ ( σ ) Λ + Λ (cid:0) + δ sin φ ( σ ) (cid:1)(cid:19) d σ + µΣ ∫ (cid:219) x ( σ ) d σ, (19)where the Lagrange multiplier µ imposes the displacementcondition and has the interpretation of a horizontal, compres-sive force.Upon substituting for φ using Eq. (13), expanding in inversepowers of Ξ , discarding terms that vanish upon integration,and integrating by parts, we find L = ∫ (cid:20) Λ + Λ + (cid:219) ψ ΛΞ − δΛ (cid:219) ψ Ξ + (cid:219) ψ ΛΞ + (cid:220) ψ ΛΞ + µ cos ψΛ (cid:18) + (cid:219) ψ Ξ + (cid:219) ψ Ξ + (cid:220) ψ Ξ + (cid:219) ψ ...ψ Ξ (cid:19) + O (cid:0) Ξ − (cid:1)(cid:21) d σ. (20)To analyze the dependence of the energy on the differential tension δ , we must go beyond lowest order [34]. We therefore truncatethe expansion at fourth order to obtain a description of the epithelium in the spirit of a Landau theory. Not only do nonlinearelastic terms arise at this order of truncation, but nonlocal terms appear, too: the theory is not elastic [35], since the energydepends not only on strain (i.e. curvature), but also on its (spatial) derivatives, introducing a nonlocal dependence on strain.To obtain the governing equation, we vary the truncated expansion (20) with respect to ψ , noting that Λ is a constant since thetrapezoidal cells are required to match up exactly [36]. After a considerable amount of algebra, we find ....ψ = Ξ (cid:220) ψ − ∆Ξ (cid:219) ψ (cid:220) ψ + µ cos ψ + + µ cos ψ ( + µ cos ψ ) (cid:219) ψ (cid:220) ψ + µΞ sin ψ + µ cos ψ (cid:20) − Ξ (cid:219) ψ − Ξ (cid:18) (cid:219) ψ − (cid:220) ψ − (cid:219) ψ ...ψ (cid:19)(cid:21) , (21)wherein ∆ = δΛ , subject to the boundary conditions ψ ( ) = ψ ( ) = , (cid:220) ψ ( ) = (cid:220) ψ ( ) = , (22a)and the integral condition ∫ cos ψ (cid:18) + (cid:219) ψ Ξ + (cid:219) ψ Ξ + (cid:220) ψ Ξ + (cid:219) ψ ...ψ Ξ (cid:19) d σ = Λ ( − D ) . (22b)The last condition imposes the fixed end-to-end shortening ofthe epithelium. These equations have a trivial solution ψ = Λ = ( − D ) − , corresponding to the compressed but unbuckledstate of the epithelium.We note that, although Eq. (21) only depends on δ and Λ through their agglomerate ∆ , a separate dependence on Λ arisesin condition (22b). Minimising the energy of buckled solutionsof Eqs. (21,22) with respect to Λ finally determines Λ . III. BUCKLING ANALYSIS
In this section, we analyze the buckling behavior of theepithelium, first determining the threshold for buckling ana-lytically and then discussing the post-buckling behavior usinga weakly nonlinear analysis of the governing equations.The buckling analysis naturally divides into two parts: wefirst seek buckled configurations of small amplitude for eachvalue of Λ , and then minimize the energy of these configura-tions with respect to Λ . A. Solution of the buckling problem
The form of the trivial solution and of condition (22b) sug-gest that the appropriate small parameter for the first part ofthe analysis is ε = − Λ ( − D ) . (23)We therefore expand ψ ( σ ) = ε (cid:16) ψ ( σ ) + εψ ( σ ) + ε ψ ( σ ) + O (cid:0) ε (cid:1)(cid:17) , (24) and write µ = µ + ε µ + ε µ + O (cid:0) ε (cid:1) . It is important tonote that, while the governing equations derived above are onlyvalid in the limit Ξ (cid:29)
1, this parameter is not an asymptoticparameter for the buckling analysis. It will however be usefulto introduce ξ = π Ξ (cid:28) . (25)Next, we solve Eq. (21), subject to the boundary and integralconditions (22), order-by-order.
1. Solution at order O ( ε ) At lowest order, the problem becomes ....ψ − Ξ (cid:220) ψ − Ξ µ + µ ψ = , (26)subject to ψ ( ) = ψ ( ) = (cid:220) ψ ( ) = (cid:220) ψ ( ) =
0. The lowesteigenvalue of this problem is µ = z − z , where z = ξ (cid:18) + ξ (cid:19) , (27)the corresponding solution for ψ being ψ ( σ ) = Ψ sin π σ. (28)
2. Solution at order O (cid:0) ε (cid:1) At next order, upon substituting for µ , ...ψ − Ξ (cid:220) ψ − z Ξ ψ = µ Ξ ( − z ) ψ − ∆Ξ ( − z ) (cid:219) ψ (cid:220) ψ . (29)subject to ψ ( ) = ψ ( ) = (cid:220) ψ ( ) = (cid:220) ψ ( ) =
0. These conditions imply that µ =
0, which is the usual result for thesupercritical pitchfork bifurcation expected for this buckling problem. Thence ψ ( σ ) = Ψ sin π σ + ∆Ψ (cid:18) ξ − ξ (cid:19) sin 2 π σ + O (cid:0) ξ (cid:1) . (30)
3. Solution at order O (cid:0) ε (cid:1) Finally, upon substituting for µ and µ = ....ψ − Ξ (cid:220) ψ − z Ξ ψ = Ξ µ ( − z ) ψ − ∆Ξ ( − z ) (cid:16) (cid:219) ψ (cid:220) ψ + (cid:220) ψ (cid:219) ψ (cid:17) + (cid:16) + z (cid:17) (cid:219) ψ (cid:220) ψ − z ( − z ) Ξ ψ − z ψ (cid:16) Ξ (cid:219) ψ − (cid:220) ψ − (cid:219) ψ ...ψ (cid:17) , (31)subject to ψ ( ) = ψ ( ) = (cid:220) ψ ( ) = (cid:220) ψ ( ) =
0. After a considerable amount of algebra, we obtain µ = Ψ (cid:20) ξ + (cid:18) − ∆ (cid:19) ξ (cid:21) + O (cid:0) ξ (cid:1) . (32)and thence ψ ( σ ) = Ψ sin π σ + Ψ Ψ ∆ (cid:18) ξ − ξ (cid:19) sin 2 π σ + Ψ (cid:18) + ∆ − ξ − ∆ − ξ (cid:19) sin 3 π σ + O (cid:0) ξ (cid:1) . (33)
4. Calculation of the amplitudes
The constants Ψ ,Ψ ,Ψ left undetermined by the above cal-culation are obtained by expanding both sides of the integralcondition (22b). Solving order-by-order, we obtain Ψ = √ − z , Ψ = , (34)where we have chosen Ψ > Ψ = ψ ≡
0; here, a non-zero ψ is required because of thesymmetry breaking resulting from the term proportional to (cid:219) ψ in the Lagrangian (20). Further, we obtain Ψ = + (cid:18) − ∆ (cid:19) ξ + (cid:18) + ∆ (cid:19) ξ + O (cid:0) ξ (cid:1) . (35) B. Minimization of the Energy
In the second part of the buckling analysis, we determinethe buckling threshold, and then analyze the post-bucklingbehavior. Substituting for Λ using Eq. (23) in the energyterm in Eq. (20) and expanding, the energy of the buckledconfiguration is (cid:18) − D + − D (cid:19) + (cid:18) − D − z − − D (cid:19) ε + O (cid:0) ε (cid:1) , (36)wherein the first bracketed term is the energy of the trivialsolution ψ = Λ = ( − D ) − . Accordingly, buckled config- urations become energetically favorable if1 − D − z − − D < ⇐⇒ D > D ∗ ≡ − √ − z . (37)In particular, the buckling threshold is independent of the dif-ferential tension δ .We are left to determine the value of Λ that minimizes theenergy of the buckled configuration. This is equivalent withrelating ε to the excess compression d = D − D ∗ >
0. Wetherefore write ε = ε d / + O ( d ) , and obtain an expansion ofthe energy in d (cid:28) E (cid:96) = (cid:32) C ε − ε − z (cid:33) d + O (cid:0) d / (cid:1) , (38)where C = + ξ + (cid:18) − δ (cid:19) ξ + O (cid:0) ξ (cid:1) . (39)The energy is thus minimized when ε = (cid:0) C ( − z ) (cid:1) − / .Substituting this result into the expression for µ , we finallyobtain µ ∼ µ ( ξ ) + (cid:18) ξ + − δ ξ + O (cid:0) ξ (cid:1)(cid:19) ( D − D ∗ ) . (40)This is the main result of our asymptotic analysis of the buck-ling: the force required to compress the epithelium decreaseswith increasing differential tension δ > (cid:219) ψ > δ >
0) and unfavorableregions ( (cid:219) ψ < δ > δ = . . . . . D .
005 0 . . . D ∗ ( − D ) − Λ D Λ D Λ . . . . . . . . . . D µ .
005 0 . D ∗ D µ D µ . . . ∂ ˙ x + ( ) ∂Λ = D ≈ ( r / π ) δ r = Ξ/‘ D (a) (b) (c)FIG. 3. Numerical buckling results. (a) Plot of Λ against relative compression D . Above a critical compression D (cid:52) (buckled shape for D = D (cid:52) shown), solution shapes at the energy minimum begin to self-intersect (dotted line for D > D (cid:52) ). Inset: zoomed plot of Λ (filled marks)against D close to the buckling threshold D ∗ . (b) Plot of compressive force µ against relative compression D , showing numerical results (solidline and dotted line for D > D (cid:52) ) in agreement with asymptotic results (dashed line). Inset: zoomed plot of µ against D close to D ∗ . Parametervalues for numerical calculations: Ξ = δ = (cid:96) = √
10. (c) Critical compression D (cid:52) against r = Ξ / (cid:96) , for different values of δ , at fixed (cid:96) = √
10, and approximation (44) thereof. Insets show buckled shapes at D = D (cid:52) and δ =
1, for different values of Ξ . IV. POST-BUCKLING BEHAVIOR
While asymptotic analysis can describe the deformationsof the epithelium just beyond the buckling threshold, largercompressions must be studied numerically. We solve the gov-erning equation (21), complemented by the boundary and in-tegral conditions (22), numerically using the boundary-value-problem solvers bvp4c , bvp5c of Matlab (The MathWorks,Inc.) and the continuation package auto [37]. A. Transition to Constricted Cells
For the numerical solution, we fix Ξ and δ , and obtain solu-tions for different values of Λ . By interpolation, we determinethe value of Λ that minimizes the energy (Fig. 3a). Thencewe obtain the corresponding value of the compressive force µ (Fig. 3b), in agreement with the asymptotic results of the previ-ous section. This also validates our numerical implementationof the system.There is, however, one extra constraint that has not beenincorporated into the continuum equations: the constraint,related to the transition to constricted cells that we have brieflydiscussed before, that the lateral sides of the trapezoidal cellscannot self-intersect. At the level of the continuum description,this constraint translates to the condition that the apical andbasal surfaces of the epithelium cannot self-intersect. Theapical and basal surfaces of the cell sheet are described by thecurves x ± = x ∓ (cid:96) sin ψ, y ± = y ± (cid:96) cos ψ, (41)where (cid:96) = (cid:96) Λ cos φ is the local thickness of the cell sheet. Itis important to note that Eqs. (41) are exact equations sincethe Kirchhoff ‘hypothesis’ of the analysis of slender elasticstructures, the asymptotic result [30] that the normal to the undeformed midline remains normal in the deformed config-uration, is an exact result in the discrete model that underliesour continuum theory. (For this same reason, an analogousanalysis for an elastic object beyond asymptotically small de-formations would, rather more intricately, require solving forthe stretches in each parallel to the midline.)Numerically, we find that, as D is increased at fixed Λ ,the shapes of minimal energy self-intersect above a criticalcompression D (cid:52) . The numerical solutions also reveal thatself-intersections first arise at σ =
0, when (cid:219) x + = ψ ( ) = (cid:219) φ ( ) =
0, ofwhich the latter follows from Eq. (13) by symmetry, f ( ) Λ − Λ(cid:96) Ξ (cid:219) ψ ( ) cos φ ( ) = , (42)where f is defined as in Eq. (16a). We note that, with thiscondition, an explicit dependence on both (cid:96) and Ξ has arisenfor the first time in our analysis. Estimating the critical compression D (cid:52) The numerical data in Fig. 3a,b suggest that the asymptoticresults of the previous section can approximate the bucklingbehavior of the epithelium well up to compressions as large as D (cid:52) . We therefore use our asymptotic results to estimate thecritical compression D (cid:52) . For this purpose, we treat r = Ξ / (cid:96) as an O ( ) quantity. Then, using Λ = ( − z ) − / + O ( d ) , f ( ) Λ − Λ(cid:96) Ξ (cid:219) ψ ( ) cos φ ( ) = √ − z − π ε Ψ r √ − z d / + O ( d ) , (43)whence, to lowest order in ξ , D (cid:52) ≈ r π . (44) Λ E ˙ x ( ) Ξ > Ξ ∗ ( δ ) Λ E ˙ x ( ) Ξ < Ξ ∗ ( δ ) . . . . . . . . . . . D ( Ξ ) D ( Ξ ) Ξ > Ξ ∗ ( δ ) Ξ < Ξ ∗ ( δ ) ( − D ) − D Λ . . . . . . . . . . . . D ( Ξ ) D ( Ξ ) Ξ < Ξ ∗ ( δ ) Ξ > Ξ ∗ ( δ ) D µ/ξ (a)(b) (c) (d)FIG. 4. Buckling for D > D (cid:52) . Two scenarios are possible: (a) If ∂ (cid:219) x ( )/ ∂Λ > Ξ < Ξ ∗ ( δ ) , buckled shapes with increased Λ do notself-intersect. (b) If ∂ (cid:219) x ( )/ ∂Λ < Ξ > Ξ ∗ ( δ ) , buckled shapes with decreased Λ do not self-intersect. Numerical results: (c) Plot of Λ (for buckled shapes of minimal energy without self-intersections) against relative compression D , for Ξ < Ξ ∗ (thick lines) and Ξ > Ξ ∗ (thinlines). For Ξ < Ξ ∗ , Λ → ( − D ) − for D > D (cid:52) . Insets show buckled shapes at D = D (cid:52) . (d) Corresponding plots of scaled compressiveforce µ / ξ against D . Dotted lines for D > D (cid:52) correspond to self-intersecting shapes at the energy minimum. Parameter values for numericalcalculations: δ = Ξ = Ξ = (cid:96) = √ This approximation is not itself an asymptotic result, yet, forsmall enough values of r , it compares well to numerical es-timates of D (cid:52) obtained by a bisection search (Fig. 3c). Thenumerical results also show that, at fixed Ξ , D (cid:52) decreaseswith increasing δ , with relative variations of about 10% for therange of δ under consideration.We also find numerically that, at large values of r , stericinteractions between different parts of the cell sheet becomeimportant before D reaches D (cid:52) . A detailed analysis of theseinteractions is beyond the scope of this discussion.We have tacitly assumed that, at fixed Ξ and at fixed D > D ∗ ,the energy E has a single local minimum as a function of Λ . This is indeed the case for small enough values of δ ,but fails at compressions D > D (cid:52) as δ is increased, so thispossibility is not of direct relevance to the present discussion.Interestingly, eigenmodes (buckled solutions with zero force)of the epithelium arise at large δ . These eigenmodes are notenergy minimizers, but, for completeness, we discuss thesesolutions in Appendix B. B. Buckled shapes for D > D (cid:52) As D is increased beyond D (cid:52) , we might expect fans of con-stricted cells to expand around the trough and (later) the crestof the buckled shape, but deriving the equations describingthese fans and solving for these shapes is beyond the scope ofthe present discussion.Here, we note simply that, for D > D (cid:52) , buckled shapeswithout self-intersections can be found. Two scenarios canbe envisaged a priori , depending on the sign of ∂ (cid:219) x ( )/ ∂Λ atthe energy minimum (Fig. 4a,b): if ∂ (cid:219) x ( )/ ∂Λ >
0, buckledsolutions without self-intersection arise as Λ is increased; if ∂ (cid:219) x ( )/ ∂Λ >
0, such shapes are found as Λ is decreased. Inter-estingly, both of these possibilities arise in the system: thereexists a critical value Ξ = Ξ ∗ ( δ ) at which ∂ (cid:219) x ( )/ ∂Λ =
0. The first possibility occurs in the case
Ξ < Ξ ∗ , and the second onein the case Ξ > Ξ ∗ (Fig. 4c).If Ξ < Ξ ∗ , the buckling amplitude decreases for these so-lutions as D is increased beyond D (cid:52) : Λ tends to its value ( − D ) − in the unbuckled, compressed configuration, while µ decreases (Fig. 4c,d). By contrast, if Ξ > Ξ ∗ , the bucklingamplitude is increased: as D increases and Λ decreases, µ increases faster than in the self-intersecting configurations atthe energy minimum (Fig. 4c,d).While these qualitative considerations cannot capture theexact mechanics of the fans of constricted cells near the throughand crest of the buckled shape, we expect them to give aqualitative indication of the buckling behavior as D is increasedjust beyond D (cid:52) . For larger values of D , there are more intricatepossibilities: there are in general two values of Λ such that (cid:219) x ( ) =
0, on either side of the energy minimum. One of thesesolutions defines the branch shown in Fig. 4c, but the secondsolution may become energetically favorable over the first oneas D is increased sufficiently. This actually happens on thebranch with Ξ = Ξ < Ξ ∗ in Fig. 4c at D ≈ .
44, but, for thesecond solution, different parts of the cell sheet start to touchbefore this value of D is reached and hence we do not pursuethis further here. V. CONCLUSION
In this paper, we have derived, by taking a rigorous asymp-totic limit, the continuum limit of a simple discrete differential-tension model of a two-dimensional epithelium. If the expan-sion is carried to high enough order for the differential tensionbetween the apical and basal sides of the epithelium to arisein the energy, nonelastic terms that are nonlocal in the strainsappear. This is the key lesson to be drawn from taking the con-tinuum limit. We have gone on to use this continuum model tostudy the buckling of the epithelium under imposed confine-ment, showing how, post-buckling, the compressive force isreduced with increasing differential tension. A second buck-ling transition occurs when constricted cells start to form nearthe troughs and crests of the buckled shape; we have discussedthe behavior close to this transition qualitatively. Taking theanalysis of the buckling behavior of epithelium in this contin-uum framework beyond the transition to constricted cells is thekey challenge for future work on this problem.Possible extensions of the continuum framework includemimicking the setup of studies of the discrete model [6, 9, 11]by coupling the epithelium to an elastic substrate or incorpo-rating fixed-volume constraints for a closed one-dimensionalepithelium. The question of how to extend the continuummodel to describe a two-dimensional epithelium also remainsopen. In particular, how are the deviations from elasticityaffected by the increase of the dimensionality of the system?Cell sheet deformations during development commonly fea-ture large geometric deformations, but the elastic deformationscan remain small provided that the deformed geometry remainsclose to the intrinsic geometry that is generally different fromthe initial geometry because of cell shape changes, cell interca-lation, and related processes. For this reason, developmentalevents as intricate as the inversion process of the green alga
Volvox can be modelled quantitatively using a Hookean shelltheory [38, 39]. By contrast, large deformations of many abiological material are not in general described well by neo-Hookean constitutive equations, although other families ofhyperelastic constitutive equations predict behavior in quan-titative agreement with experimental data for brain and fattissues [40–42]. However, and in spite of the ubiquity of theseelastic models, in particular in the modelling of the folds of the cerebrum [15–19], it was pointed out very recently that thefolding of the cerebellum is fundamentally inconsistent withthe differential-growth hypothesis [43]: in the cerebellum, theoscillations of the thicknesses of the core and the growingcortex are out-of-phase, while elastic bilayer instabilities leadto in-phase oscillations [43]. All of this emphasizes the needfor a deeper understanding of how continuum models relate toproperties of structures at the cell level. By explicitly show-ing how both nonlinear and nonlocal elastic terms arise in thecontinuum limit of a simple discrete model and impact onits behavior, the present analysis has taken a first step in thisdirection.
ACKNOWLEDGMENTS
We thank Matej Krajnc for discussions of Ref. [10] at anearly stage of this work, and are grateful for support fromthe Engineering and Physical Sciences Research Council (Es-tablished Career Fellowship EP/M017982/1, REG; DoctoralPrize Fellowship, PAH), Wellcome Trust Investigator Award207510/Z/17/Z, the Schlumberger Chair Fund, and MagdaleneCollege, Cambridge (Nevile Research Fellowship, PAH).
Appendix A: Derivation of Eq. (12)
In this appendix, we relate φ and ψ by deriving Eq. (12)used in the main text. At the same time, we verify that thisexpansion is indeed consistent at all orders. Taylor expandingthe right-hand side of Eq. (11), ψ ( s + ks ) − ψ ( s ) = φ ( s ) + (cid:169)(cid:173)(cid:171) k − (cid:213) j = φ ( s + js ) (cid:170)(cid:174)(cid:172) + φ ( s + ks ) = k φ ( s ) + ∞ (cid:213) n = s n n ! φ ( n ) ( s ) k n + k − (cid:213) j = j n = k φ ( s ) + ∞ (cid:213) n = s n n ! φ ( n ) ( s ) k n + n + n (cid:213) j = (− ) j (cid:18) n + j (cid:19) B j ( k − ) n + − j , (A1)where we have used Faulhaber’s formula [44] to expand the sum of powers of integers, and where B = , B = − , . . . denotethe Bernoulli numbers (of the first kind) [32]. Expanding ( k − ) n + − j using the binomial theorem and simplifying the binomialcoefficients, ψ ( s + ks ) − ψ ( s ) = k φ ( s ) + ∞ (cid:213) n = φ ( n ) ( s ) n ! (cid:96) n A ( k , n ) , (A2)where we have introduced A ( k , n ) = k n + (− ) n + n ! n (cid:213) j = n + − j (cid:213) i = (− ) i B j k i i ! j ! ( n + − i − j ) ! = n + (cid:213) i = a i ( n ) k i , (A3)wherein the coefficients a , a , . . . , a n + depend on n . We notice in particular that a = , a n = , a n + = n + , (A4)of which the last two are obtained by direct computation, and the first one follows using an identity of Bernoulli numbers [32], n (cid:213) j = (cid:18) n + j (cid:19) B j = n = , , . . . . (A5)Moreover, for i = , , . . . , n − a i ( n ) = (− ) n + − i n ! i ! n + − i (cid:213) j = B j j ! ( n + − i − j ) ! . (A6)Accordingly, A ( k , n ) = k n + n + + (− ) n + n ! n − (cid:213) i = (− ) i k i i ! n + − i (cid:213) j = B j j ! ( n + − i − j ) ! . (A7)Upon inverting the order of summation, Eq. (A2) becomes ψ ( s + ks ) − ψ ( s ) = k φ ( s ) + ∞ (cid:213) n = φ ( n ) ( s ) n ! (cid:96) n (cid:18) k n + n + (cid:19) + ∞ (cid:213) i = ∞ (cid:213) n = i + (− ) n + − i φ ( n ) ( s ) (cid:96) n k i i ! n + − i (cid:213) j = B j j ! ( n + − i − j ) ! . (A8)Finally, upon relabelling indices in the first summation and changing variables n (cid:55)−→ m = n + − i in the last summation, ψ ( s + ks ) − ψ ( s ) = ∞ (cid:213) i = k i i ! (cid:96) i (cid:96) φ ( i − ) ( s ) + ∞ (cid:213) m = (− ) m (cid:96) m − φ ( i − + m ) ( s ) (cid:169)(cid:173)(cid:171) m (cid:213) j = B j j ! ( m − j ) ! (cid:170)(cid:174)(cid:172) . (A9)But, rearranging Eq. (A5), m (cid:213) j = B j j ! ( m − j ) ! = B m m ! for m = , , . . . . (A10)Since B n = n >
1, and using B =
1, we finallyobtain ψ ( s + ks ) − ψ ( s ) = ∞ (cid:213) i = k i i ! (cid:96) i (cid:40) ∞ (cid:213) m = B m ( m ) ! φ ( i − + m ) ( s ) (cid:96) m − (cid:41) . (A11)Comparing this to the Taylor expansion of the left-hand side, ψ ( s + ks ) − ψ ( s ) = ∞ (cid:213) i = k i i ! (cid:96) i ψ ( i ) ( s ) , (A12)we deduce that the expansion is consistent at all orders, with,in particular, ψ (cid:48) ( s ) = ∞ (cid:213) m = ψ m φ ( m ) ( s ) (cid:96) m − where ψ m = B m ( m ) ! , (A13)which is Eq. (12). As noted in the main text, we are not awareof a closed form for the inverted series that expresses φ as afunction of the derivatives of ψ . Formally, inverting Eq. (12)gives φ ( s ) = ∞ (cid:213) m = φ m ψ ( m + ) ( s ) (cid:96) m + , (A14)where the coefficients φ , φ , . . . are determined recursivelyby φ ψ = m (cid:213) j = φ j ψ m − j = m = , , . . . . (A15) In agreement with Eq. (13), we find φ = , φ = − , φ = , . . . . (A16) Appendix B: Eigenmodes of the buckled epithelium
Eigenmodes of the epithelium are non-zero solutions of thegoverning equation (21) with µ =
0. They thus obey ....ψ = Ξ (cid:220) ψ − ∆Ξ (cid:219) ψ (cid:220) ψ + (cid:219) ψ (cid:220) ψ, (B1) . . . . . ∆ ∗ ( Ξ )
10 20 30 4044 . . . . Ξ∆ ∗ no s o l u ti on s ∆ D FIG. 5. Eigenmodes of a buckled epithelium: plot of relative end-to-end shortening D against ∆ . Parameter value: Ξ =
20. Noeigenmodes were found for ∆ < ∆ ∗ ( Ξ ) . On the dashed part of thebranch, E <
2. Continuation failed at the point marked × . Inset: plotof ∆ ∗ against Ξ (solid line) and power-law fit (dashed line). ψ ( ) = ψ ( ) = , (cid:220) ψ ( ) = (cid:220) ψ ( ) = . (B2)To find eigenmodes numerically, we remove the trivial, zerosolution by imposing a non-zero compression D and varyingthis compression until a solution with µ = ∆ (cid:62) ∆ ∗ , but findno solutions if ∆ < ∆ ∗ , for some value ∆ ∗ depending on Ξ (Fig. 5). Some of these solutions have energy E < D > ∆ ∗ against Ξ (Fig. 5, inset) suggests that ∆ ∗ approaches a constant value as Ξ grows large. We observe that the numericaldata are well approximated by a power-law ∆ ∗ = c + c Ξ − / ,where c ≈ . c ≈ . ∆ and hence δ for these eigenmodesbeckon a comment on the formal range of validity of the contin-uum model: stability of the underlying discrete model requires α, β (cid:62) δ (cid:54) (cid:96) . While the asymptotic expan-sion leading to the geometric relation (13) was an expansionin the large parameter (cid:96) , it did not involve δ . By contrast, theexpansion of the Lagrangian (20), which did involve δ , wasan expansion in a different large parameter, Ξ . Hence largevalues of δ (cid:46) (cid:96) are indeed in the formal range of validity ofthe continuum model provided that Ξ is large enough. [1] T. Lecuit and P.-F. Lenne, “Cell surface mechanics and thecontrol of cell shape, tissue patterns and morphogenesis,” Nat.Rev. Mol. Cell Biol. , 633–644 (2007).[2] G.M. Odell, G. Oster, P. Alberch, and B. Burnside, “Themechanical basis of morphogenesis,” Dev. Biol. , 446–462(1981).[3] M. Rauzi, A. Hočevar Brezavšček, P. Ziherl, and M. Leptin,“Physical models of mesoderm invagination in Drosophila em-bryo,” Biophys. J. , 3–10 (2013).[4] A. G. Fletcher, F. Cooper, and R. E. Baker, “Mechanocellularmodels of epithelial morphogenesis,” Phil. Trans. R. Soc. B ,20150519 (2016).[5] J. Derganc, S. Svetina, and B. Žekš, “Equilibrium mechanics ofmonolayered epithelium,” J. Theor. Biol. , 333–339 (2009).[6] A. Hočevar Brezavšček, M. Rauzi, M. Leptin, and P. Ziherl, “Amodel of epithelial invagination driven by collective mechanicsof identical cells,” Biophys. J. , 1069–1077 (2012).[7] M. Krajnc, N. Štorgel, A. Hočevar Brezavšček, and P. Ziherl,“A tension-based model of flat and corrugated simple epithelia,”Soft Matter , 8368–8377 (2013).[8] E. Hannezo, J. Prost, and J.-F. Joanny, “Theory of epithelialsheet morphology in three dimensions,” Proc. Nat. Acad. Sci.USA , 27–32 (2014).[9] M. Rauzi, U. Krzic, T. E. Saunders, M. Krajnc, P. Ziherl, L. Huf-nagel, and M. Leptin, “Embryo-scale tissue mechanics during Drosophila gastrulation movements,” Nat. Commun. , 8677(2015).[10] M. Krajnc and P. Ziherl, “Theory of epithelial elasticity,” Phys.Rev. E , 052713 (2015).[11] N. Štorgel, M. Krajnc, P. Mrak, J. Štrus, and P. Ziherl, “Quanti-tative morphology of epithelial folds,” Biophys. J. , 269–277(2016).[12] M. Krajnc, S. Dasgupta, P. Ziherl, and J. Prost, “Fluidizationof epithelial sheets by active cell rearrangements,” Phys. Rev. E , 022409 (2018).[13] J. Howard, S. W. Grill, and J. S. Blois, “Turing’s next steps: themechanochemical basis of morphogenesis,” Nat. Rev. Mol. CellBiol. , 392–398 (2011).[14] K. Dasbiswas, E. Hannezo, and N. S. Gov, “Theory of epithe-lial cell shape transitions induced by mechanoactive chemicalgradients,” Biophys. J. , 968–977 (2018).[15] D. P. Richman, R. M. Stewart, J. W. Hutchinson, and V. S.Caviness, Jr., “Mechanical model of brain convolutional devel-opment,” Science , 18–21 (1975).[16] T. Tallinen, J. S. Biggins, and L. Mahadevan, “Surface sulci in squeezed soft solids,” Phys. Rev. Lett. , 024302 (2013).[17] O. V. Manyuhina, D. Mayett, and J. M. Schwarz, “Elasticinstabilities in a layered cerebral cortex: a revised axonal tensionmodel for cortex folding,” New J. Phys. , 123058 (2014).[18] S. Budday, P. Steinmann, and E. Kuhl, “Secondary instabilitiesmodulate cortical complexity in the mammalian brain,” Philos.Mag. , 3244–3256 (2015).[19] T. Tallinen, J. Y. Chung, F. Rousseau, N. Girard, J. Lefève, andL. Mahadevan, “On the growth and form of cortical convolu-tions,” Nat. Phys. , 588–593 (2016).[20] E. Lejeune, A. Javili, J. Weickenmeier, E. Kuhl, and C. Linder,“Tri-layer wrinkling as a mechanism for anchoring center initi-ation in the developing cerebellum,” Soft Matter , 5613–5620(2016).[21] E. Hannezo, J. Prost, and J.-F. Joanny, “Instabilities of mono-layered epithelia: Shape and structure of villi and crypts,” Phys.Rev. Lett. , 078104 (2011).[22] T. Savin, N. A. Kurpios, A. E. Shyer, P. Florescu, H. Liang,L. Mahadevan, and C. J. Tabin, “On the growth and form of thegut,” Nature , 57–62 (2011).[23] A. E. Shyer, T. Tallinen, N. L. Nerurkar, Z. Wei, E. S. Gil, D. L.Kaplan, C. J. Tabin, and L. Mahadevan, “Villification: how thegut gets its villi,” Science , 212–218 (2013).[24] L. Hufnagel, A. A. Teleman, H. Rouault, S. M. Cohen, andB. I. Shraiman, “On the mechanism of wing size determinationin fly development,” Proc. Nat. Acad. Sci. USA , 3835–3840(2007).[25] R. Farhadifar, J.-C. Röper, B. Aigouy, S. Eaton, and F. Jülicher,“The influence of cell mechanics, cell-cell interactions, andproliferation on epithelial packing,” Curr. Biol. , 2095–2104(2007).[26] S. Hilgenfeldt, S. Erisken, and R. W. Carthew, “Physical mod-eling of cell geometric order in an epithelial tissue,” Proc. Nat.Acad. Sci. USA , 907–911 (2008).[27] D. Bi, J. H. Lopez, J. M. Schwarz, and M. L. Manning, “Adensity-independent rigidity transition in biological tissues,”Nat. Phys. , 1074–1079 (2015).[28] D. Bi, X. Yang, M. C. Marchetti, and M. L. Manning, “Motility-driven glass and jamming transitions in biological tissues,” Phys.Rev. X , 021011 (2016).[29] L. D. Landau and E. M. Lifshitz, Theory of Elasticity , 2nded., Course of Theoretical Physics, Vol. 7 (Pergamon, Oxford,England, 1970) Chap. 21, pp. 97–100.[30] B. Audoly and Y. Pomeau,
Elasticity and Geometry (OxfordUniversity Press, Oxford, England, 2010) App. A, pp. 546–557; Chap. 6.4, pp. 184–192 and App. D, pp. 571–581.[31] In this respect, the present analysis differs from that of Ref. [10],where the energy was expressed as an integral with respectto the arclength S in the deformed configuration, related infact to s by (cid:96) d S = (cid:96) d s . In Ref. [10], a weaker global areaconservation constraint for the cell sheet had therefore to beimposed separately.[32] J. Spanier and K. B. Oldham, An Atlas of Functions (Hemi-sphere, Washington, DC, USA, 1987) Chap. 4, pp. 35–38.[33] These and other expansions in this paper were carried out withthe help of Mathematica (Wolfram, Inc.) to facilitate the ma-nipulation of complicated algebraic expressions.[34] A non-trivial term arises already at order O (cid:0) Ξ − (cid:1) in the secondpart of Eq. (20) that imposes the condition of fixed displacement,stemming from the corrections at order O (cid:0) (cid:96) − (cid:1) in Eqs. (16).Hence, already at this order that does not even resolve the effectof non-zero differential tension δ , the governing equation differsfrom Euler’s Elastica equation.[35] A. Libai and J. G. Simmonds, The Nonlinear Theory of Elas-tic Shells , 2nd ed. (Cambridge University Press, Cambridge,England, 1998) Chap. III.J, pp. 36–38.[36] In Ref. [10], the energy was expressed in terms of ψ and (cid:96) ,which were then both varied, while they are in fact related by (cid:96) sec φ ( ψ ) = const. for matching cells. If the lateral sides arenot required to match up exactly, corrections to the final termin Eq. (2) would have to be introduced, however, to ensurea consistent description of the adhesion between neighboringcells. Notwithstanding this, minimising Eq. (6) with respect to λ leads to e ( φ ) = (cid:96) √ sec φ + δ tan φ , with e (cid:48)(cid:48) ( ) = (cid:0) − δ (cid:1) (cid:96) / δ > √ Auto-07p: Continuation and Bifurca-tion Software for Ordinary Differential Equations , Tech. Rep.(Concordia University, Montreal, Canada, 2012).[38] S. Höhn, A. R. Honerkamp-Smith, P. A. Haas, P. Khuc Trong,and R. E. Goldstein, “Dynamics of a