Nonlinear and nonlocal elasticity in coarse-grained differential-tension models of epithelia

The shapes of epithelial tissues result from a complex interplay of contractile forces in the cytoskeleta of the cells in the tissue and adhesion forces between them. A host of discrete, cell-based models describe these forces by assigning different surface tensions to the apical, basal, and lateral sides of the cells. These differential-tension models have been used to describe the deformations of epithelia in different living systems, but the underlying continuum mechanics at the scale of the epithelium are still unclear. Here, we derive a continuum theory for a simple differential-tension model of a two-dimensional epithelial monolayer and study the buckling of this epithelium under imposed compression. The analysis reveals how the cell-level properties encoded in the differential-tension model lead to linear and nonlinear elastic as well as nonlocal, nonelastic behavior at the continuum level.


I. INTRODUCTION
Intercellular adhesion proteins and cortical actin networks are well established as regulators of cell surface mechanics, and hence of the deformations of epithelia during morphogenesis [1].Ever since the seminal work of Odell et al. [2], these cellular components have therefore underlain mathematical models of epithelia [3,4].One large class of such models are differential-tension models [5][6][7][8][9][10][11][12], in which cell polarity, cell-cell adhesion properties and the actomyosin network induce different surface tensions in different sides of the discrete cells.Coupling the mechanical models describing epithelial deformations to models of the intracellular biochemistry is a key challenge in the field [13], but some progress has recently been made by coupling the differential-tension model of Ref. [8] to the diffusion of a 'mechanogen' that induces contractility [14].
While the differential-tension models can quantitatively reproduce the morphology of epithelial folds in many different living systems [11], it is likely that, in general multilayered epithelia, the formation of epithelial folds must be ascribed to a combination of these intra-epithelial stresses and differential growth of different parts of the tissue.Models based on the latter only have for example been invoked to describe, at the scale of the epithelium, the formation of cortical convolutions in the brain [15][16][17][18][19][20] and of the intestinal villi [21][22][23], the 'fication' of which lends itself to the pun that gave Ref. [23] its title.The coarse-grained limit of the differential-tension models at this scale is less well-studied, however, and this continuum limit is the topic of this paper.
We shall focus on the most basic setup of these differentialtension 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 .The internal state of the cells induces different surface tensions Γ a , Γ b , and Γ in the apical, basal, and lateral sides of the cells, respectively.The energy of a single cell therefore reads where the factor of 1/2 has been introduced for mere convenience [7,10].The theory can be extended to incorporate additional physics such as a basement membrane [11] or a confining vitelline membrane [6].In a full three-dimensional setup, this leads to the study of the shapes of prism-shaped cells [8,12].Here, we shall restrict to the two-dimensional setup [7,10] of an epithelium consisting of isosceles trapezoidal 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 (Fig. 1).Since there are two such lateral sides, the cell energy (1) reduces to per unit extent in the third dimension [7,10].A different twodimensional limit is obtained by averaging over the thickness of the cell sheet and describing in-plane deformations only.Such models, termed area-and perimeter-elasticity models, have been studied extensively [24][25][26][27][28].
The simplest problem in the mechanics of elastic rods is their Euler buckling under applied forces [29]; it is therefore meet to ask how the buckling behavior of an active material such as this model epithelium differs from that of an elastic rod.This problem was considered in Ref. [10], where the continuum limit of Eq. ( 2) was mapped to Euler's Elastica equation [30].These calculations, complementing simulations of the discrete model in Ref. [7], were phrased in terms of spontaneous buckling of the epithelium, but an additional compressive or extensile force is required to produce these deformations, making the analysis of Ref. [10] more appropriate to the present setup of buckling under imposed forces.Moreover, the analysis of Ref. [10] is not completely consistent with the discrete model, since it does not impose the condition that the trapezoidal cells match up exactly along their lateral sides.In this paper, we perform a consistent asymptotic expansion of the discrete geometry of this model, revealing nonlinear and nonlocal elastic terms in the continuum limit.We then analyze the buckling behavior of the continuum model under imposed compression analytically and numerically.

A. Single-cell Energy
To obtain the energy of a single cell, we express the sidelengths of the trapezoidal cells in terms of their mean base K, their height L, and the angle 2φ that their lateral sides make with each other (Fig. 2a): Incompressibility implies the cell area conservation constraint Upon eliminating K and L using these relations, the energy of a single cell is expressed, from Eq. ( 2) and as a function of L and φ, as

Non-dimensionalization
We non-dimensionalize this expression by scaling lengths with the square root of the cell area and thus define the nondimensional length of the lateral sides of the trapezoidal cells, λ = L /A 1/2 c .We further set = L/A 1/2 c and κ = K/A 1/2 c .Finally, following Ref.[10], we introduce the parameters We note that 0 is the (uniform and non-dimensionalized) thickness of the epithelium in the flat configuration.Hence s 0 = 1/ 0 is the (non-dimensionalized) width of a single cell (i.e. its arclength in the flat configuration).The nondimensional energy e = E/Γ A 1/2 c of a single cell is therefore Without loss of generality, we assume that δ > 0 in what follows, so, for a single cell, φ < 0 is energetically favorable.

Transition to Constricted Cells
The transition to constricted triangular cells is a geometric singularity in the discrete model.These triangular cells arize as limiting cases of the trapezoidal cells when L a = 0 or L b = 0. Using Eq. (3), the conditions L a , L b 0 reduce to

B. Energy of an Epithelium
In the continuum limit, we take φ to be a function of the arclength s of the midline of the undeformed, flat epithelium.Summing over all the cells, we obtain the non-dimensional energy E of the epithelium, By imposing an energy density equal to 1/s 0 = 0 times the (non-dimensional) energy of a single cell and integrating with respect to the reference arclength in this manner, we have imposed local cell area conservation [31].
The boundary conditions are most naturally expressed in terms of the angle ψ of the deformed midline of the epithelium below the horizontal (Fig. 2b).We therefore express the energy in terms of ψ.This is usefully done in the scaling limit 0 1 of a columnar epithelium, as explained next.

Asymptotic Expansion
We make two further scaling assumptions: These scalings correspond to the regime ∼ 0 and φ 1, where the cells deform but slightly from their equilibrium configuration.In this limit, Eq. (3c) implies that λ ∼ ∼ 0 .Further, area conservation (3d) requires that κ ∼ 1/ 0 , and thus, from Eqs. (3a,3b), we must have φ κ/λ ∼ 1/ 2 0 .The second scaling thus corresponds to the largest deformations allowed.We therefore introduce the parameter We are now set up to relate φ and ψ, for which purpose we use the geometric relation valid for any positive integer k, as sketched in Fig. 2b.In Appendix A, we show that, with our scaling assumptions, the continuum limit of this relation is wherein B 0 = 1, B 1 = − 1 2 , . . .are the Bernoulli numbers (of the first kind) [32].The next step is to invert this series, to express φ in terms of the derivatives of ψ.While we are not aware of any explicit expression for the coefficients of the inverted series, it is straightforward to invert the series orderby-order by substituting back and forth, and thus obtain where dashes denote differentiation with respect to s.In Ref. [10], only the first term of this expansion was obtained.Inclusion of the second term will enable us to analyze the buckling behavior of the epithelium in what follows.

Shape Equations for the Buckled Epithelium
We describe the shape of the buckled epithelium by the coordinates x(s), y(s) of the centreline of the epithelium, defined by the axes in Fig. 2b.To derive the continuum equations describing the centreline, we begin by projecting the discrete geometry onto the axes, Using κ(s) = ( 0 Λ) −1 sec φ(s) and expanding these equations order-by-order in inverse powers of 0 using Eq. ( 13), we obtain, after a considerable amount of algebra [33], where Integrating these differential equations yields the shape of the buckled epithelium.Deviations from the 'standard' values

Derivation of the Governing Equation
We shall seek to describe buckled configurations of an epithelium of undeformed length 2Σ s 0 .We change variables by introducing σ = s/Σ, use dots to denote differentiation with respect to σ, and define Since Ξ = Σ/s 0 , the number of cells in the epithelium is simply N = 2Ξ.We shall seek buckled solutions with clamped boundary conditions and a prescribed relative compression D, so that where the coordinates are now expressed relative to the scaled arclength σ.We shall restrict to symmetrically buckled configurations for which ψ(σ) = −ψ(2−σ).The second condition above is then satisfied.We may further reduce the solution to the range 0 σ 1, with the condition of prescribed compression reading x(1) − x(0) = 1 − D. To minimize the energy of the epithelium at this imposed displacement, we therefore consider the Lagrangian where the Lagrange multiplier µ imposes the displacement condition and has the interpretation of a horizontal, compressive force.Upon substituting for φ using Eq. ( 13), expanding in inverse powers of Ξ, discarding terms that vanish upon integration, and integrating by parts, we find To analyze the dependence of the energy on the differential tension δ, we must go beyond lowest order [34].We therefore truncate the expansion at fourth order to obtain a description of the epithelium in the spirit of a Landau theory.Not only do nonlinear elastic terms arise at this order of truncation, but nonlocal terms appear, too: the theory is not elastic [35], since the energy depends 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 the trapezoidal cells are required to match up exactly [36].After a considerable amount of algebra, we find wherein ∆ = δΛ 2 , subject to the boundary conditions and the integral condition The last condition imposes the fixed end-to-end shortening of the epithelium.These equations have a trivial solution ψ = 0, Λ = (1−D) −1 , corresponding to the compressed but unbuckled state of the epithelium.We note that, although Eq. ( 21) only depends on δ and Λ through their agglomerate ∆, a separate dependence on Λ arises in condition (22b).Minimising the energy of buckled solutions of Eqs.(21,22) with respect to Λ finally determines Λ.

III. BUCKLING ANALYSIS
In this section, we analyze the buckling behavior of the epithelium, first determining the threshold for buckling analytically and then discussing the post-buckling behavior using a weakly nonlinear analysis of the governing equations.
The buckling analysis naturally divides into two parts: we first seek buckled configurations of small amplitude for each value of Λ, and then minimize the energy of these configurations with respect to Λ.

A. Solution of the buckling problem
The form of the trivial solution and of condition (22b) suggest that the appropriate small parameter for the first part of the analysis is We therefore expand and write It is important to note that, while the governing equations derived above are only valid in the limit Ξ 1, this parameter is not an asymptotic parameter for the buckling analysis.It will however be useful to introduce Next, we solve Eq. ( 21), subject to the boundary and integral conditions (22), order-by-order.

Calculation of the amplitudes
The constants Ψ 0 ,Ψ 1 ,Ψ 2 left undetermined by the above calculation are obtained by expanding both sides of the integral condition (22b).Solving order-by-order, we obtain where we have chosen Ψ 0 > 0 without loss of generality.The result Ψ 1 = 0 is to be expected for a supercritical bifurcation; in fact, in a standard elastic buckling problem, one would have ψ 1 ≡ 0; here, a non-zero ψ 1 is required because of the symmetry breaking resulting from the term proportional to ψ 3 in the Lagrangian (20).Further, we obtain

B. Minimization of the Energy
In the second part of the buckling analysis, we determine the buckling threshold, and then analyze the post-buckling behavior.Substituting for Λ using Eq. ( 23) in the energy term in Eq. ( 20) and expanding, the energy of the buckled configuration is wherein the first bracketed term is the energy of the trivial solution ψ = 0, Λ = (1 − D) −1 .Accordingly, buckled config-urations become energetically favorable if In particular, the buckling threshold is independent of the differential tension δ.
We are left to determine the value of Λ that minimizes the energy of the buckled configuration.This is equivalent with relating ε to the excess compression d = D − D * > 0. We therefore write ε = ε 0 d 1/2 + O(d), and obtain an expansion of the energy in d 1, where The energy is thus minimized when . Substituting this result into the expression for µ 2 , we finally obtain This is the main result of our asymptotic analysis of the buckling: the force required to compress the epithelium decreases with increasing differential tension δ > 0.

IV. POST-BUCKLING BEHAVIOR
While asymptotic analysis can describe the deformations of the epithelium just beyond the buckling threshold, larger compressions must be studied numerically.We solve the governing equation ( 21), complemented by the boundary and integral conditions (22), numerically using the boundary-valueproblem solvers bvp4c, bvp5c of M (The MathWorks, Inc.) and the continuation package [37].

A. Transition to Constricted Cells
For the numerical solution, we fix Ξ and δ, and obtain solutions for different values of Λ.By interpolation, we determine the value of Λ that minimizes the energy (Fig. 3a).Thence we obtain the corresponding value of the compressive force µ (Fig. 3b), in agreement with the asymptotic results of the previous section.This also validates our numerical implementation of the system.
There is, however, one extra constraint that has not been incorporated into the continuum equations: the constraint, related to the transition to constricted cells that we have briefly discussed before, that the lateral sides of the trapezoidal cells cannot self-intersect.At the level of the continuum description, this constraint translates to the condition that the apical and basal surfaces of the epithelium cannot self-intersect.The apical and basal surfaces of the cell sheet are described by the curves where = 0 Λ cos φ is the local thickness of the cell sheet.It is important to note that Eqs. ( 41) are exact equations since the Kirchhoff 'hypothesis' of the analysis of slender elastic structures, the asymptotic result [30] that the normal to the undeformed midline remains normal in the deformed configuration, is an exact result in the discrete model that underlies our continuum theory.(For this same reason, an analogous analysis for an elastic object beyond asymptotically small deformations would, rather more intricately, require solving for the 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 critical compression D .The numerical solutions also reveal that self-intersections first arise at σ = 0, when x + = 0 there.Expanding this condition using ψ(0) = 0 and φ(0) = 0, of which the latter follows from Eq. ( 13) by symmetry, where f is defined as in Eq. (16a).We note that, with this condition, an explicit dependence on both 0 and Ξ has arisen for the first time in our analysis.

Estimating the critical compression D
The numerical data in Fig. 3a,b suggest that the asymptotic results of the previous section can approximate the buckling behavior of the epithelium well up to compressions as large as D .We therefore use our asymptotic results to estimate the critical compression D .For this purpose, we treat r = Ξ/ 2 0 as an O(1) quantity.Then, using whence, to lowest order in ξ,  This approximation is not itself an asymptotic result, yet, for small enough values of r, it compares well to numerical estimates of D obtained by a bisection search (Fig. 3c).The numerical results also show that, at fixed Ξ, D decreases with increasing δ, with relative variations of about 10% for the range of δ under consideration.
We also find numerically that, at large values of r, steric interactions between different parts of the cell sheet become important before D reaches D .A detailed analysis of these interactions 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 as δ is increased, so this possibility 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 not energy minimizers, but, for completeness, we discuss these solutions in Appendix B.

B. Buckled shapes for D > D
As D is increased beyond D , we might expect fans of constricted cells to expand around the trough and (later) the crest of the buckled shape, but deriving the equations describing these fans and solving for these shapes is beyond the scope of the present discussion.
Here, we note simply that, for D > D , buckled shapes without self-intersections can be found.Two scenarios can be envisaged a priori, depending on the sign of ∂ x(0)/∂Λ at the energy minimum (Fig. 4a,b): if ∂ x(0)/∂Λ > 0, buckled solutions without self-intersection arise as Λ is increased; if ∂ x(0)/∂Λ > 0, such shapes are found as Λ is decreased.Interestingly, both of these possibilities arise in the system: there exists a critical value Ξ = Ξ * (δ) at which ∂ x(0)/∂Λ = 0.The first possibility occurs in the case Ξ < Ξ * , and the second one in the case Ξ > Ξ * (Fig. 4c).
If Ξ < Ξ * , the buckling amplitude decreases for these solutions as D is increased beyond D : Λ tends to its value (1 − D) −1 in the unbuckled, compressed configuration, while µ decreases (Fig. 4c,d).By contrast, if Ξ > Ξ * , the buckling amplitude is increased: as D increases and Λ decreases, µ increases faster than in the self-intersecting configurations at the energy minimum (Fig. 4c,d).
While these qualitative considerations cannot capture the exact mechanics of the fans of constricted cells near the through and crest of the buckled shape, we expect them to give a qualitative indication of the buckling behavior as D is increased just beyond D .For larger values of D, there are more intricate possibilities: there are in general two values of Λ such that x(0) = 0, on either side of the energy minimum.One of these solutions defines the branch shown in Fig. 4c, but the second solution may become energetically favorable over the first one as D is increased sufficiently.This actually happens on the branch with Ξ = Ξ 1 < Ξ * in Fig. 4c at D ≈ 0.44, but, for the second solution, different parts of the cell sheet start to touch before this value of D is reached and hence we do not pursue this further here.

V. CONCLUSION
In this paper, we have derived, by taking a rigorous asymptotic limit, the continuum limit of a simple discrete differentialtension model of a two-dimensional epithelium.If the expansion is carried to high enough order for the differential tension between the apical and basal sides of the epithelium to arise in the energy, nonelastic terms that are nonlocal in the strains appear.This is the key lesson to be drawn from taking the continuum limit.We have gone on to use this continuum model to study the buckling of the epithelium under imposed confine-ment, showing how, post-buckling, the compressive force is reduced with increasing differential tension.A second buckling transition occurs when constricted cells start to form near the troughs and crests of the buckled shape; we have discussed the behavior close to this transition qualitatively.Taking the analysis of the buckling behavior of epithelium in this continuum framework beyond the transition to constricted cells is the key challenge for future work on this problem.
Possible extensions of the continuum framework include mimicking the setup of studies of the discrete model [6,9,11] by coupling the epithelium to an elastic substrate or incorporating fixed-volume constraints for a closed one-dimensional epithelium.The question of how to extend the continuum model to describe a two-dimensional epithelium also remains open.In particular, how are the deviations from elasticity affected by the increase of the dimensionality of the system?
Cell sheet deformations during development commonly feature large geometric deformations, but the elastic deformations can remain small provided that the deformed geometry remains close to the intrinsic geometry that is generally different from the initial geometry because of cell shape changes, cell intercalation, and related processes.For this reason, developmental events as intricate as the inversion process of the green alga Volvox can be modelled quantitatively using a Hookean shell theory [38,39].By contrast, large deformations of many a biological material are not in general described well by neo-Hookean constitutive equations, although other families of hyperelastic constitutive equations predict behavior in quantitative agreement with experimental data for brain and fat tissues [40][41][42].However, and in spite of the ubiquity of elastic models, in particular in the modelling of the folds of the cerebrum [15][16][17][18][19], it was pointed out very recently that the folding of the cerebellum is fundamentally inconsistent with the differential-growth hypothesis [43]: in the cerebellum, the oscillations of the thicknesses of the core and the growing cortex are out-of-phase, while elastic bilayer instabilities lead to in-phase oscillations [43].All of this emphasizes the need for a deeper understanding of how continuum models relate to properties of structures at the cell level.By explicitly showing how both nonlinear and nonlocal elastic terms arise in the continuum limit of a simple discrete model and impact on its behavior, the present analysis has taken a first step in this direction.To find eigenmodes numerically, we remove the trivial, zero solution by imposing a non-zero compression D and varying this compression until a solution with µ = 0 is found.Numerically, we obtain eigenmodes if ∆ ∆ * , but find no solutions if ∆ < ∆ * , for some value ∆ * depending on Ξ (Fig. 5).Some of these solutions have energy E < 2 (Fig. 5), lower than the energy of the uncompressed, flat solution; these are spontaneous buckled modes that arise in the absence of external forces, but, as is apparent from the corresponding values D > 1 (Fig. 5), these solutions are unphysical.Plotting ∆ * against Ξ (Fig. 5, inset) suggests that ∆ * approaches a constant value as Ξ grows large.observe that the numerical data are well approximated by a power-law ∆ * = c 1 + c 2 Ξ −5/4 , where c 1 ≈ 3.96, c 2 ≈ 19.5 (Fig. 5, inset).
The 'large' values of ∆ and hence δ for these eigenmodes beckon a comment on the formal range of validity of the continuum model: stability of the underlying discrete model requires α, β 0 [7], and hence δ 2 0 .While the asymptotic expansion leading to the geometric relation (13) was an expansion in the large parameter 0 , it did not involve δ.By contrast, the expansion of the Lagrangian (20), which did involve δ, was an expansion in a different large parameter, Ξ. Hence large values of δ 2 0 are indeed in the formal range of validity of the continuum model provided that Ξ is large enough.

FIG. 1 .
FIG. 1. Model epithelium.Isosceles trapezoidal cells of apical and basal bases L a and L b are connected along their lateral sides, which have length L .The continuous line and shaded area provide a cartoon of the continuum limit.

FIG. 2 .
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 , the lateral sides being at an angle 2φ to each other.(b) Definition of the tangent angle ψ of the midline below the horizontal.The geometry of contiguous cells defines the relation between φ and ψ, as expressed in Eq. (11).

FIG. 3 .
FIG. 3. Numerical buckling results.(a) Plot of Λ against relative compression D. Above a critical compression D (buckled shape for D = D shown), solution shapes at the energy minimum begin to self-intersect (dotted line for D > D ).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 (solid line and dotted line for D > D ) in agreement with asymptotic results (dashed line).Inset: zoomed plot of µ against D close to D * .Parameter values for numerical calculations: Ξ = 20, δ = 1, 0 = √ 10.(c) Critical compression D against r = Ξ/ 2 0 , for different values of δ, at fixed 0 = √ 10, and approximation (44) thereof.Insets show buckled shapes at D = D and δ = 1, for different values of Ξ.