InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Computer and Information Science » Numerical Analysis and Scientific Computing » "Perusal of the Finite Element Method", book edited by Radostina Petrova, ISBN 978-953-51-2820-5, Print ISBN 978-953-51-2819-9, Published: December 14, 2016 under CC BY 3.0 license. © The Author(s).

Chapter 6

Improved Reduced Order Mechanical Model and Finite Element Analysis of Three-Dimensional Deformations of Epithelial Tissues

By Ara S. Avetisyan, Asatur Zh. Khurshudyan and Sergey K. Ohanyan
DOI: 10.5772/65027

Article top


Ventral furrow formation in Drosophila according to [3] (upper) and [4] (lower).
Figure 1. Ventral furrow formation in Drosophila according to [3] (upper) and [4] (lower).
Apical constriction of ventral furrow at different times [10].
Figure 2. Apical constriction of ventral furrow at different times [10].
Ventral furrow formation in Drosophila according to [11].
Figure 3. Ventral furrow formation in Drosophila according to [11].
Epithelial cell aspect ratio as a bistable phenomenon [27]: αl is the contractile force and Λa is the actin-myosin tension.
Figure 4. Epithelial cell aspect ratio as a bistable phenomenon [27]: αl is the contractile force and Λa is the actin-myosin tension.
Diagrams of a single-cell layer: cross section (left) and 3D shape (right).
Figure 5. Diagrams of a single-cell layer: cross section (left) and 3D shape (right).
Single-cell model (left) and adherens junction model (right): diagonals of the top hexagon imitate actin fibres, red areas between cells imitate junction bonds.
Figure 6. Single-cell model (left) and adherens junction model (right): diagonals of the top hexagon imitate actin fibres, red areas between cells imitate junction bonds.
Apical constriction of initially flat tissue. It consists from 5050 elements or cells and has 529,245 DOFs.
Figure 7. Apical constriction of initially flat tissue. It consists from 5050 elements or cells and has 529,245 DOFs.
von Mises stress distribution in deformed configurations.
Figure 8. von Mises stress distribution in deformed configurations.
Apical constriction of initially cylindrical tissue. It consists from 5975 elements or cells and has 621,375 DOFs.
Figure 9. Apical constriction of initially cylindrical tissue. It consists from 5975 elements or cells and has 621,375 DOFs.
von Mises stress distribution in deformed configurations.
Figure 10. von Mises stress distribution in deformed configurations.
Apical constriction of initially spherical tissue. It consists from 3126 elements or cells and has 325,245 DOFs.
Figure 11. Apical constriction of initially spherical tissue. It consists from 3126 elements or cells and has 325,245 DOFs.
von Mises stress distribution in deformed configurations.
Figure 12. von Mises stress distribution in deformed configurations.

Improved Reduced Order Mechanical Model and Finite Element Analysis of Three-Dimensional Deformations of Epithelial Tissues

Ara S. Avetisyan, Asatur Zh. Khurshudyan and Sergey K. Ohanyan
Show details


In this chapter, we analyse non-uniform bending of single-layer cell tissues—epithelia, surrounding organs throughout the body. Dimensionally reduced model is suggested, which is equivalent to membranes with bending stiffness: the total elastic energy of the tissue is a combination of stretching and bending energies. The energy, suggested in this chapter, is a piecewise function, the branches of which correspond to a specific deformation regime: compression, pure bending and stretch.

Under specific assumptions on a single cell, such as all cells are identical, elastic, homogeneous, have the shape of a hexagonal prism, are linked to each other by pure elastic junctions, and so on, tissues of different shapes are analysed using finite element discretization. The deformation of initially flat, cylindrical and spherical tissues is visualized to reproduce similar deformations in real epithelial tissues observed during experiments.

Keywords: morphogenesis, embryonic tissue, adherens junctions, relaxation, quasi-convex energy densities, isometric deformations, short maps, tissue mechanics, Γ-convergence, lower semi-continuity, hexagonal elements

1. Introduction

Morphogenesis, the evolutionary process of shape and structure development of organisms or their parts, is driven by certain types of cell shape changes or, in other words, deformations. One type of such deformations is the apical constriction—visible shrinkage of the apical side of cells, leading to bending of epithelia—cell sheets, which surround organs throughout the body. First occurring at early stages of embryogenesis, the apical constriction results initially flat epithelia to be bent obtaining three-dimensional form, which, depending on physiological context and morphogenetic stage, leads to different consequences. It has been first hypothesized in [1] that in various developmental systems the apical constriction may drive the bending of epithelia. Using bulky mechanical construction, the hypothesis of Rhumbler is first practically tested in [2]. The testing mechanism consists of 13 identical bars (cell walls) that are kept apart by means of stiff tubes connecting their centres (cell kernels), and are held in a row by rubber bands connecting their ends (cell membrane). In the first stage, the rubber bands at both sides are stretched equally and the mechanism is in straight equilibrium. Shortening the rubber bands at, for instance, the upper side uniformly in each segment (cell), the mechanism is bent on the side of the greater tension, as it could be expected. Thus, Rhumbler’s hypothesis could take place.

A more illustrative, computer model for verifying Rhumbler’s hypothesis is suggested in [3], where bars and bands are replaced by virtual cells. Prescribing certain mechanical properties to the cell membrane and cytoplasmic components, and assuming that as a result of deformations the volume of a single cell is preserved, which is implicitly done by Rhumbler and Lewis, a model of cuboidal epithelia folding is suggested, based on the local behaviour of individual cells. The model is demonstrated on the example of ventral furrow formation in Drosophila. A sequence of cells in the form of a cylindrical shell, representing the cross section of ventral furrow, is deformed such that apical constriction occurs in its lower row of cells. Increasing the applied stresses, different steps of the furrow formation are illustrated. Figure 1 expresses how much Odell’s model is close to the real microscopic picture [4].


Figure 1.

Ventral furrow formation in Drosophila according to [3] (upper) and [4] (lower).

Much is known about the causes of the apical constriction, but some issues still remain unexplored [58]. The most studied causes include contraction of actin filaments—fibrous network localized at the cortex of the cell, by interacting with motor protein myosin—the most known converter of chemical energy into mechanical work. Under influence of myosin, actin fibres contract leading to shrinkage of the cell in the area of their localization [9, 10]. The contraction size and principal direction of the fibres, that is, the microscopical deformation of a single cell, mainly depend on the tissue type: actin-myosin network contraction may deform columnar (occurring, for instance, in digestive tract and female reproductive system) or cuboid (resp. in kidney tubules) cells into trapezoidal-, wedge- or bottle-shaped cells. The deformation size strongly depends also on the emplacement of deforming cells: cells with different placement are constricted differently, so that macroscopically one observes localized wrinkles (see Figure 2).


Figure 2.

Apical constriction of ventral furrow at different times [10].

Because of the mechanical nature of cell shape change, in particular, apical constriction, its theoretical study first of all must rely on mechanical principles and constitutive laws. In early mechanical models, the tissue is modelled as a continuum material, so that the position and behaviour of individual cells are unimportant, that is, only macroscopic deformation of the tissue is studied. Moreover, the height of the cells (viz. the thickness of the tissue) is supposed to be negligible with respect to other measures, such that the deformation of the tissue can be described in terms of its mid-surface. In such models, the actin network is not explicitly accounted and the contraction forces are modelled as a force term acting on the outer surface of the epithelium. The actin network is explicitly accounted in [11], where thin elastic shell model based on linear Cauchy relations is derived to describe apical constriction in initially non-flat epithelia. The model is tested to simulate apical constriction in initially flat, cylindrical and spherical tissues. Ventral furrow formation in Drosophila is simulated (see Figure 3), and it becomes evident from comparison of Figures 1 and 3 that the model, being three-dimensional, is in good correspondence with that from [3]: Figure 1 corresponds to cross section of cylindrical shell from Figure 3. However, it does not involve the behaviour of individual cells and describes only macroscopic deformation. For further introduction into other mechanical models, we refer to [1126].


Figure 3.

Ventral furrow formation in Drosophila according to [11].

The physical and geometrical characteristics of individual cells are taken into account by mechanical model suggested in [27] describing three-dimensional deformations of cell sheet. Both actin-myosin network contractile tensions and cell-cell adhesion stresses are involved to contribute in epithelial tissue bending. Prescribing specific mechanical characteristics to cytoplasmic components and to the kernel, the effective energy of a cell, which is assumed to be a hexagonal prism with constant volume, is expressed in terms of its basal length. As a result of simulations, it is particularly established that when adhesion stresses, distributed on lateral sides of the cell, are large enough, the effective energy has two local minima, that is, there are two equilibrium cell shapes. More particular, depending on position, cell shape may be discontinuously transited from squamous to columnar forms. Figure 4 shows that there is a whole domain in the plane of lateral adhesion and actin-myosin tensions for which the effective energy has two local minima.


Figure 4.

Epithelial cell aspect ratio as a bistable phenomenon [27]: αl is the contractile force and Λa is the actin-myosin tension.

In mathematical terms, it means that the total energy of the tissue is not lower semi-continuous, which implies that there is no convergence of its minimizers, and, therefore, it becomes impossible to derive real deformation of the tissue [28, 29]. To be more illustrative, consider a membrane, that is, a tissue without bending stiffness, which is compressed by boundary stresses. The only way that the membrane can accommodate a compressed state is due to wrinkling [30]. Wrinkles may occur more and more finely, so the limiting deformation will be smooth, but evidently it will not minimize the membrane energy.

To overcome such difficulties, the membrane energy density is usually substituted by its quasi-convex envelope [28], which ensures the lower semi-continuity of its total energy. In calculus of variations, there are several ways to construct or even compute convex envelopes for functionals, depending on first- and second-order derivatives [28, 3136].

Another possibility for incorporating the lack of lower semi-continuity of membrane energy functionals is accounting the bending stiffness of the membrane and adding to the energy the bending contribution [37]. So, it is established in [33, 34] that the total energy of a membrane with some bending stiffness and thickness h, that is, the functional


of deformation r consisting of membrane (Em) and bending (Eb) contributions is lower semi-continuous for arbitrary Em as far as Eb is lower semi-continuous. Explicit forms for Em and Eb as well as sufficient conditions on material constants for which Eb (and therefore E) is lower semi-continuous are derived ibid. Such materials are called stable.

However, the most rigorous and general approach seems to be the derivation of low-order energies from general three-dimensional non-linear elasticity by means of Γ-convergence. The concept was originally developed by De Giorgi, see [29]. It turns out that the Γ-limit of three-dimensional elastic energy functional, when the thickness of the body goes to zero, is always lower semi-continuous and since most commonly used two-dimensional energy functionals, that is, membrane energy, pure bending (Kirchhoff) energy and higher-order terms (von Kármán–like energies), are already derived as such limits (see Section 3), it becomes possible to analyze real deformations of particular tissues.

By suggested improvements, it is supposed to get rid of disadvantages of previous models and take into account the fact that the total energy of the tissue can have two local minima. Thus, we are intended

  • to account the height of each individual cell, so the model is fully three-dimensional,

  • to pick up strongly localized internal forces to model the phenomenon of apical constriction more precisely,

  • to introduce stretching and bending contributions into the total elastic energy of the epithelial tissue, thus making the model more realistic and convenient not only for qualitative but also for quantitative analysis,

  • to use the quasi-convexification of the total energy, so the real three-dimensional deformations of the epithelial tissues can be identified as minimizers of the total energy by a gradient flow technique.

The former allows making simplified toolboxes for analysing three-dimensional deformations of epithelial tissues. The rough diagram of the improved model looks like in Figure 5, so that the characteristics of each individual cell are important, as well as the actin belt is explicitly involved in the model.


Figure 5.

Diagrams of a single-cell layer: cross section (left) and 3D shape (right).

The chapter is organized as follows: In Section 2, some preliminary definitions and main notations are brought to make the chapter independent of outer sources. In Section 3, known results from Γ-convergent energies are shortly described. In Section 4, the quasi-convex envelope of membrane energy densities of general form is derived explicitly and sufficient conditions for their convexity are derived. Moreover, a specific energy is suggested for different deformation regimes that fully satisfies the needs of basic mechanical theory suggestion. Finally, in Section 5 main results of finite element analysis of tissues of particular initial shapes (flat, cylindrical, spherical) are discussed. Overall conclusion completes the body of the chapter.

2. Preliminaries

We begin with the definition of concepts used throughout the paper, which and much more on this topic can be found in [28, 29].

There are different concepts of convexity in higher dimensions, such as quasi-, poly-, rank-one and separately convexity. Here, we use all but poly-convexity, so it is of no use to enlarge the chapter by its definition.

Definition 1 (quasi-convexity). A function media/ueq1.png is said to be quasi-convex if


for every bounded open set Ω, for every media/ueq2.png and for every media/ueq3.png

Definition 2 (rank-one convexity). A function media/ueq4.png is said to be rank-one convex if


for every θ[0,1], media/ueq5.png with rank F=1.

Definition 3 (separately convexity). A function media/ueq6.png is said to be separately convex if


is convex for every i=1,,N and for every fixed media/ueq7.png

In other words, a function media/ueq8.png is separately convex if it is convex on every line parallel to a coordinate axis.

In general, convexity implies quasi-convexity, which implies rank-one convexity, and finally rank-one convexity implies separately convexity.

The quasi-convexity notion is generalized by Meyers [38] for functionals depending on higher-order derivatives. We refer to [32, 35, 36] for practical use of the definition.

Definition 4 (k-quasi convexity). A function media/ueq9.png is said to be kquasi-convex if


for every open-bounded set Ω, for every media/ueq10.png and for every media/ueq11.png

In particular, 2-quasi-convex functions, which we use here, are the restriction of quasi-convex functions to symmetric matrices [39] (in the case of k-quasi convexity, see [40]).

We widely use the quasi-convex envelope of W(G) or apparent energy density, Wq(G), which can be defined in terms of deformations of the form r(x)=Gx+r0(x) with media/ueq12.png constant, media/ueq13.png and r0(x)=0 on the boundary Ω [41, 42]:


It is evident that for any media/inline14.png and when W=0, Wq=0 as well.

Deformation is always denoted by media/inline15.png with open-bounded media/inline16.png We introduce the first and second fundamental forms of a deformation r


respectively, in which media/ueq17.png is the unit of the outer normal to Ω. The first fundamental form quantitatively characterizes the stretch or compression of the membrane, while the second fundamental form characterizes its curvature.

Definition 5. Energy density function media/inline18.png is called frame indifferent if for all rotations R3SO(3)

If for all R2SO(2),

then it is called isotropic.

In other words, frame indifferent energy densities are invariant against three-dimensional rotations before a deformation is applied and that isotropic energy densities are invariant against two-dimensional rotations after a deformation is applied.

Definition 6 (short (δ-short) deformation). A deformation media/inline19.png is called short (δ-short) if I(r)Id2×2 (δId2×2) a.e. in Ω.

Definition 7 (isometric deformation). A deformation media/inline20.png is called isometric if I(r)=Id2×2 a.e. on Ω.

Short deformations result compressive stresses, while isometric deformations do not allow stretching or compression. The case IId2×2 corresponds to tensile or, in general, non-compressive stresses. It is known that for non-compressive stresses, the membrane energy density function is already quasi-convex [41, 42].

3. Γ-convergence and Γ-limits of three-dimensional non-linear elastic energy

Rigorous derivation of two-dimensional energy functionals for thin bodies from three-dimensional functionals is of particular interest. There are several ways for dimension reduction, such as asymptotic expansion, Γ-convergence technique, and so on (for general survey, see [43]). In [44], a hierarchy of plate models is derived as Γ-limit of three-dimensional elastic energy functional when the thickness of the body h0. It turns out that the scaling of external force in h plays an important role for Γ-limit derivation. Suppose a body occupies the domain Ω×[h2,h2], and the applied forces media/inline21.png active on the outer surface of the body satisfy




denotes the rescaled elastic energy of the body, with


The total energy of the body will be


The Γ-limit


for β0 will provide constitutive mechanical models to analyse deformed three-dimensional shape of initially flat thin bodies.

We combine the results of [4448] in Theorem 1.

Theorem 1. Let the stored energy density function is frame indifferent and satisfy


and in a neighbourhood of SO(3), WC2.


(i) (Membrane theory [45]). Suppose α=β=0 and


Then media/inline25.png is a β-minimizing sequence, that is, if


then media/inline26.png (for a subsequence). The limiting deformation r¯ is independent of x3 and minimizes


among all media/inline27.png. For derivation of Wq from given W see Section 4 subsequently.

(ii) (Constrained membrane theory [46]). Suppose that 0<α<53 and β=α. Then media/ueq31.png and every β-minimizing sequence r(h) has a subsequence with media/ueq25.png in . The limit r¯ is independent of x3 and minimizes


among all short deformations media/inline27.png.

(iii) (Non-linear bending theory [47]). Suppose α=β=2. Then |infεh|Ch2 and if r(h) is a β-minimizing sequence then there is strong convergence r(h)r¯ in media/ueq28.png (for a subsequence). The limiting deformation r¯ is isometric, independent of x3, belongs to media/inline31.png and minimizes


among all isometric deformations media/inline32.png, which belong to media/inline33.png The non-linear strain satisfies


where 2symG=GT+G, amin is the solution of


In all cases, there is convergence of energy, that is,


For α>0, the limiting deformation r¯ is not only isometric but is close to a rigid motion. For that reason in [44], for deformation media/inline34.png the deformation


with R˜(h)SO(3) and media/inline35.png is introduced and the averaged in-plane and out-of-plane displacements


are rescaled by



Then, we have

Theorem 2. (iv) (Linearized isometry constraint). Suppose 2<α<3 and set β=2α2, γ=2(α2), 2δ=γ. If 2<α<52, suppose in addition that Ω is simply connected. Then, Chβinfεh0. If r(h) is a β-minimizing sequence then there exist constant R(h)SO(3) and c(h) ∈ ℝ3 such that


Moreover, the pair (v¯,R¯) minimizes the functional


subject to det2v=0.

(v) (von Kármán theory). Suppose that α=3 and set β=4, γ=2, δ=1. Then media/inline36.png and for a subsequence of a β-minimizing sequence, (1) and (2) hold and the limit triple (u¯,v¯,R¯) minimizes von Kármán functional


(vi) (Linearized von Kármán theory). Suppose α>3 and set β=2α2, γ=α1 and δ=α2. Then Chβinfεh0 and for a sequence of a β-minimizing sequence. (2) holds with u¯=0 and the pair (v¯,R¯) minimizes the linearized von Kármán functional


Moreover, the non-linear strain satisfies


In all cases, we have convergence of the rescaled energy media/inline42.png to the minimum of the limit functional media/inline43.png. Moreover, for f30 we have R¯33=1 or 1.

Remark 1. (i) Practically, Theorems 1 and 2 establish Γ-convergence and give the Γ-limits of three-dimensional non-linear energy functional in the range β[0,53)[2,). The range β[53,2) remains unexplored.

(ii) Even in the case β=53, it is still possible to find a sequence rk(h) for which the functional media/ueq35.png is bounded when h0. The limiting deformations r are composed of finitely many affine isometries separated by sharp folds and are called origami deformations [46]. It is proved ibid that the class of origami deformations is the closure of the class of short deformations with respect to the uniform convergence.

(iii) In particular, when α=β=2, for isotropic homogeneous bodies we obtain


in which λ and μ are Lamé constants, coinciding with plate energy derived by Kirchhoff much earlier in [49].

4. Relaxation of epithelial elastic energy by Pipkin’s procedure and by adding the bending contribution

As was mentioned in Section 1, there are two ways to relax the epithelial elastic energy. The first way is to relax the stretching energy for deformation regimes corresponding to compressive stresses. For that reason, Pipkin‘s procedure [41] seems to be the most common tool.

Suppose that the stretching energy density function W=W(r) is frame indifferent, isotropic and positive except the cases


where it turns to zero, that is, attains its minimum. On elastic energies with two minima, see [50]. This assumption is motivated by bistability of cell shapes during apical constriction (see Figure 4).

We restrict attention to deformations r with


and denote Wq(r):=W˜q(σ1,σ2). Here |σj|=λj, j{1;2}, are the principal stretches, the squares of which are the eigenvalues of the first fundamental form. Therefore


In view of frame indifference of W, and thence, W˜, we conclude that W˜q(±σ1,±σ2)=W˜q(λ1,λ2).

From quasi-convexity of Wq, in particular, follows its rank-one convexity, therefore since the matrix


is rank one, we have

therefore W˜q is convex with respect to λ1. Repeating the argument when


we see that W˜q(λ1,λ2) is separately convex.

Under assumptions made, we have W˜q(λ1,λ2)=0 if λ1=λ2=1 and λ1=λ2=δ. Since δ>1 and 0W˜q are convex along lines parallel to λ1=0 and λ2=0, therefore W˜q=0 on edges connecting the vertexes (±δ,±δ) in the σ1σ2 plane. Repeating the argument, we conclude that W˜q=0 inside the square Ω0={(λ1,λ2);λ1δ,λ2δ}. So, for δ-short maps Wq0.

Since W˜q(λ1,λ2) is even and convex in λ2, it will attain its minimum with respect to λ2 only when λ2=0. Denote by W˜0(λ1)=minλ2W˜(λ1,λ2) when λ1>δ. Suppose that λ0(λ1) is the largest value of λ2 for which W˜q(λ1,λ0)=W˜0(λ1), λ1>δ. The convexity of W˜q in λ2 implies that W˜q(λ1,λ2)=W˜0(λ1) for λ2λ0(λ1), λ1>δ. In the same way, using the symmetry of W˜q in its arguments, we will arrive at W˜q(λ1,λ2)=W˜0(λ2) for λ1λ0(λ2), λ2>δ. Evidently, for the rest of deformation regimes,W is convex.

If we denote by Ω1={(λ1,λ2);λ1>δ,λ2λ0(λ1)} and Ω2={(λ1,λ2);λ1λ0(λ2),λ2>δ}, the ranges of principal stretches, corresponding to uniaxial tensions in directions of those stretches, respectively [51], we finally arrive at


Besides opportunity of explicit relaxation, in [41] the derivation of criteria for convexity and quasi-convexity is described. In view of rank-one convexity, W˜q satisfies Legendre-Hadamard inequality from which it follows


For convexity of W˜q in corresponding domain instead of the last two equalities, the stronger conditions


must be satisfied. Above


The other way is the adding of pure bending contribution to the elastic energy of the tissue. Since actin-myosin network contraction leads to compressive or non-stretching stresses, we have to incorporate the elastic energy mainly for such deformations. According to Section 3 for non-stretching stresses, we have

  1. (short deformations and compressive stresses)

    Γlimh01hβE(h)=0in the rangeβ[0,53),

  2. (isometric deformations, neither stretch nor compressive stresses)


The deformation regime IId corresponding to tensile stresses leads to Em[r], which is already lower semi-continuous. In the case of uniaxial tension, there is no bending and the energy expression will look like the middle rows of Wq (42). Thence, the total energy might be given as follows:

E[r]=Ωfrdx+(0,IδId,12ΩW˜0(λ1(r))dx,for uni-axial tension in the direction ofλ1,12ΩW˜0(λ2(r))dx,for uni-axial tension in the direction ofλ2,124ΩQ2(II)dx,I=Id,12ΩWq(r)dx,otherwise.

We assume that the tissue is homogeneous and isotropic, so according to Remark 1


Furthermore, the membrane energy density function we take in the following form [52, 53]:


which is valid for large deformations of incompressible hyper-elastic membranes. Above, I1=trI is the first strain invariant, γ1>0 and γ2>0 are material parameters, evaluated experimentally in [53]. For a neo–Hookean material 2γ1=μ and 2γ2=κ, in which μ is the shear modulus and κ is the bulk modulus. W˜q and W˜0 are determined according to the procedure described above. For that reason, I1 must be represented in terms of the principal stretches λ1 and λ2 as follows:


in which the incompressibility condition λ1λ2λ3=1 is taken into account.

For any fixed λ1>1 (uniaxial tension in the first principal direction), the corresponding W˜(λ1,λ2) attains its minimum with respect to λ2 at the maximal positive root of


Since in the uniaxial tension regime λ12+λ22>1, we finally get



The correspondent tension is defined by


Since W˜ is obviously symmetric in its arguments, in the case of uniaxial tension in the second principal direction, we would have




It is evident, that when we substitute λ1=1 in the first case and λ2=1 in the second case, we arrive at undeformed configuration of the membrane, so the relaxed energy is fully consistent.

Forces driving tissue deformation are strongly localized and in general are compressive. Force f must possess those properties.

5. Finite element analysis of three-dimensional shape of some characteristic deformations of tissues during morphogenesis

In this section, we summarize main results of finite element analysis of a single layer tissue model, the elastic energy of which is given by


in which Em and Eb are defined in Eqs. (49)(51). Since the discretization procedure is standard, we omit the details and bring only the main results.

In Figure 6, we bring the model of a single-cell (element) and cell-cell junction (in red). All structures (plate and shell) considered in this section entirely consist of such cell groups. In all tissues considered below, the height of a cell h=4r and h=200d, in which r is its side and d is the diameter of actin fibres. The ratio of Young‘s moduli of a cell and actin fibres is Eact=3Ecell, and of actin fibres and cell-cell links Eact=1.5Elink.

We consider

  1. initially flat rectangular plate (Figure 7 (left)),

  2. cylindrical shell (Figure 9 (left)),

  3. spherical shell (Figure 11 (left)).


Figure 6.

Single-cell model (left) and adherens junction model (right): diagonals of the top hexagon imitate actin fibres, red areas between cells imitate junction bonds.


Figure 7.

Apical constriction of initially flat tissue. It consists from 5050 elements or cells and has 529,245 DOFs.


Figure 8.

von Mises stress distribution in deformed configurations.


Figure 9.

Apical constriction of initially cylindrical tissue. It consists from 5975 elements or cells and has 621,375 DOFs.


Figure 10.

von Mises stress distribution in deformed configurations.


Figure 11.

Apical constriction of initially spherical tissue. It consists from 3126 elements or cells and has 325,245 DOFs.

Elements of the middle part of the rectangular tissue are compressed in apical sides to imitate apical constriction in cells. Increasing the compressing stresses, the tissue is bent and a blaster-shaped pattern is formed as shown in Figure 7 (right). The quantitative picture of the stresses arising in the tissue is drawn in Figure 8.


Figure 12.

von Mises stress distribution in deformed configurations.

Next, we consider a cylindrical shell to imitate ventral furrow (generally all tubular patterns) formation. Elements of the top part of the cylinder are constrained in the apical sides by compressing the links standing for apical fibres (see Figure 9). Increasing the compressing stresses, the ventral furrow formation is simulated similar to stages presented in Figure 1. The quantitative picture of the stresses arising in the tissue is drawn in Figure 10.

Finally, a spherical shell is simulated. Cells at the top of a semi-sphere are constrained in apical sides and depending on values of compressing stresses various stages of blastopore formation in archenteron can be described (see Figure 11 (right)). The quantitative picture of the stresses arising in the tissue is drawn in Figure 12.

6. Conclusions

Analysis based on lower semi-continuous energy functionals reveals real three-dimensional deformations of soft epithelial tissues. Having different forms for different deformation regimes, such as compressive stresses (short deformations), uniaxial tensions collinear to principal directions of the first fundamental form, no stretching or pure bending stresses (isometric deformations) and stretching stresses (large deformations), the resulting total energy is lower semi-continuous, so the existence of its minimizers, that is, real deformations, is ensured. Particular energy density functions are chosen and the explicit form of the total energy functional is obtained, thereby the discretization is made easy.

On the basis of obtained energy functional, a three-dimensional discretized model of epithelial tissues undergoing combined stretching and bending deformations is constructed. Discretization elements correspond to single cells forming the tissue. Actin fibres and cell-cell adhesion links, mainly contributing on the tissue energy, are explicitly embedded in elements. Deformations characteristic to specific embryonic tissues (ventral furrow, neutral tube, neurosphere) observed earlier are described quantitatively increasing contractile stresses in fibres.


We thankfully dedicate the chapter to the blessed memory of our good fellow and colleague, a candidate of physical and mathematical sciences, Hamlet V. Hovhannisyan (1956–2016), who unexpectedly died before he could realize his best scientific ideas.


The theoretical part of the chapter is investigated under the guidance of Doctor, Professor Benedikt Wirth, Institute for Computational and Applied Mathematics, University of Münster, whom we are heartily thankful. The work of As. Kh. and S. O. was made possible in part by a research grant from the Armenian National Science and Education Fund (ANSEF) based in New York, NY, USA.


Inner (dot) product
λ,μLamé coefficients
λ1,λ2Principal stretches
EhRescaled three-dimensional total energy
μ(Ω)Measure of bounded open set Ω
Gradient (nabla) operator
ΩMidsurface of the epithelium
rank FRank of the matrix F
Tensor product
distUsual distance in three-dimensional Euclidean space
dDiameter of actin fibres
EThree-dimensional elastic energy
EbBending energy
EhRescaled three-dimensional elastic energy
EmMembrane energy
EactYoung’s modulus of actin fibres
EcellYoung’s modulus of a cell
ElinkYoung’s modulus of cell-cell links
hThickness of the epithelium or a single cell
rSide of a single cell
WMembrane energy density function
WqThe quasi-convex envelope of W
IIThe second fundamental form associated with deformation r
IThe first fundamental form associated with deformation r
nUnit normal vector
rDeformation acting from R2 to R3
Id2×22×2 identity matrix
SO(3)The group of all rotations about the origin of three-dimensional Euclidean space


1 - Rhumbler L., Zur Mechanik des Gastrulationsvorganges insbesondere der Invagination Eine entwickelungsmechanische Studie. Archiv für Entwicklungsmechanik der Organismen, 1902, vol. 14, issue 3, pp. 401–476.
2 - Lewis W. H., Mechanics of invagination. The Anatomical Record, 1947, vol. 97, issue 2, pp. 139–156.
3 - Odell G. M., Oster G., Alberch P., Burnside B., The mechanical basis of morphogenesis: I. Epithelial folding and invagination. Developmental Biology, 1981, vol. 85, issue 2, pp. 446–462.
4 - Polyakov O., Bing H., Michael S., Joshua W. S., Matthias K., Wieschaus E., Passive mechanical forces control cell-shape change during Drosophila ventral furrow formation. Biophysical Journal, 2014, vol. 107, issue 4, pp. 998–1010.
5 - Keller R., Developmental biology. Physical biology returns to morphogenesis. Science, 2012, vol. 338, pp. 201–203.
6 - Mammoto T., Ingber D. E., Mechanical control of tissue and organ development. Development, 2010, vol. 137, pp. 1407–1420.
7 - Martin A. C., Goldstein B., Apical constriction: themes and variations on a cellular mechanism driving morphogenesis. Development, 2014, vol. 141, pp. 1987–1998.
8 - Sawyer J. M., Harrell J. R., Shemer G., Sullivan-Brown J., Roh-Johnson M., Goldstein B., Apical constriction: a cell shape change that can drive morphogenesis. Developmental Biology, 2010, vol. 341, pp. 5–19.
9 - Martin A. C., Pulsation and stabilization: contractile forces that underlie morphogenesis. Developmental Biology, 2010, vol. 341, pp. 114–25.
10 - Martin A. C., Kaschube B., Wieschaus E. F., Pulse contractions of an actin-myosin network drive apical constriction. Nature, 2009, vol. 457, issue 22, pp. 495–501.
11 - Jones G. W., Chapman S. J., Modelling apical constriction in epithelia using elastic shell theory. Biomechanics and Modeling in Mechanobiology, 2010, vol. 9, issue 3, pp. 247–261.
12 - Ciarletta P., Ben Amar M., Labouesse M., Continuum model of epithelial morphogenesis during Caenorhabditis elegans embryonic elongation. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2009, vol. 367, pp. 3379–3400.
13 - von Dassow M., Davidson L., Variation and robustness of the mechanics of gastrulation: the role of tissue mechanical properties during morphogenesis. Birth Defects Research Part C: Embryo Today, 2007, vol. 81, issue 4, pp. 253–269.
14 - Du. X., Osterfield M., Shvartsman S. Y., Computational analysis of three-dimensional epithelial morphogenesis using vertex models. Physical Biology, 2014, vol. 11, issue 6, 15 p.
15 - Fletcher A. G., Osterfield M., Baker R. E., Shvartsman S. Y., Vertex models of epithelial morphogenesis. Biophysical Journal, 2014, vol. 106, issue 11, pp. 2291–2304.
16 - Imai M., Furusawa K., Mizutani T, Kawabata A., Haga H., Three-dimensional morphogenesis of MDCK cells induced by cellular contractile forces on a viscous substrate. Scientific Reports, 2015, vol. 5, 14208.
17 - Ingber D. E., Mechanical control of tissue morphogenesis during embryological development. International Journal of Developmental Biology, 2006, vol.50, pp. 255–266.
18 - Inoue Y., Suzuki M., Watanabe T., Yasue N., Tateo I., Adachi T., Ueno N., Mechanical roles of apical constriction, cell elongation, and cell migration during neural tube formation in Xenopus. Biomechanics and Modeling in Mechanobiology, 2016, DOI 10.1007/s10237-016-0794-1.
19 - Jauvert S., Peyroux R., Richefeu V., A mechanical model for cell motility and tissue morphogenesis. Computer Methods in Biomechanics and Biomedical Engineering, 2013, vol. 16, pp. 13–14.
20 - Lecuit T., Lenne P.-F., Cell surface mechanics and the control of cell shape, tissue patterns and morphogenesis. Nature Reviews Molecular Cell Biology, 2007, vol. 8, pp. 633–644.
21 - Murisic N., et al., From discrete to continuum models of three-dimensional deformations in epithelial sheets. Biophysical Journal, 2015, vol. 109, issue 1, pp. 154–163.
22 - Patwari P., Lee R. T., Mechanical control of tissue morphogenesis. Circulation Research, 2008, vol. 103, pp. 234–243.
23 - Rauzi M., Krzic U., Saunders T. E., Krajnc M., Ziherl P., Hufnagel L., Leptin M., Embryo-scale tissue mechanics during Drosophila gastrulation movements. Nature Communications, 2015, vol. 6, Article number: 8677, doi: 10.1038/ncomms9677.
24 - Siedlik M. J., Nelson C. M., Mechanics of tissue morphogenesis. In “Cell and Matrix Mechanics” edited by Kaunas R., Zemel A., CRC Press, 2014.
25 - Vaughan B. L. Jr., Baker R. E., Kay D., Maini P. K., A modified Oster-Murray-Harris mechanical model of morphogenesis. SIAM Journal of Applied Mathematics, 2013, vol. 73, issue 6, pp. 2124–2142.
26 - Wyczalkowski M. A., Chen Z., Filas B. A., Varner V. D., Taber L. A., Computational models for mechanics of morphogenesis. Birth Defects Research Part C: Embryo Today: Reviews, 2012. vol. 96, issue 2, pp. 132–152.
27 - Hannezo E., Prost J., Joanny J.-F., Theory of epithelial sheet morphology in three dimensions. PNAS, vol. 111, issue 1, pp. 27–32.
28 - Dacorogna B., Direct Methods in the Calculus of Variations. Springer, 2009.
29 - Dal Maso G., An Introduction to Γ-Convergence. Birkhäuser, Boston, 1993.
30 - Steigman D. J., Pipkin A. C., Finite deformations of wrinkled membranes. Quarterly Journal of Mechanics and Applied Mathematics, 1989, vol. 42, pp. 427–440.
31 - Kohn R. V., Strang G., Optimal design and relaxation of variational problems I, II, III. Communications on Pure and Applied Mathematics, 1986, vol. 39, pp. 1–25 (I), 139–182 (II), 353-377 (III).
32 - Hilgers M. G., Pipkin A. C., Elastic sheets with bending stiffness. Quarterly Journal of Mechanics and Applied Mathematics, 1992, vol. 45, issue 1, pp. 57–75.
33 - Hilgers M. G., Pipkin A. C., Bending energy of highly elastic membranes. Quarterly of Applied Mathematics, 1992, vol. L, issue 2, pp. 389–400.
34 - Hilgers M. G., Pipkin A. C., Bending energy of highly elastic membranes II. Quarterly of Applied Mathematics, 1996, vol. LIV, issue 2, pp. 307–316.
35 - Hilgers M. G., Pipkin A. C., Energy–minimizing deformations of elastic sheets with bending stiffness. Journal of Elasticity, 1993, vol. 31, pp. 125–139.
36 - Hilgers M. G., Pipkin A. C., The Graves condition for variational problems of arbitrary order. IMA Journal of Applied Mathematics, 1992, vol. 48, pp. 265–269.
37 - Schmidt B., Fraternali F., Universal formulae for the limiting elastic energy of membrane networks. Journal of the Mechanics and Physics of Solids, 2012, vol. 60, pp. 172–180.
38 - Meyers N. G., Quasi-convexity and lower semi-continuity of multiple variational integrals of any order. Transactions of AMS, 1965, vol. 119, pp. 125–149.
39 - Dal Maso G., Fonseca I., Leoni G., Morini M., Higher-order quasi-convexity reduces to quasi-convexity. Archive for Rational Mechanics and Analysis, 2004, vol. 171, pp. 55–81.
40 - Cagnetti F., k-Quasi-convexity reduces to quasi-convexity. Proceedings of the Royal Society of Edinburgh: Section A Mathematics, 2011, vol. 141, issue 4, pp. 673–708.
41 - Pipkin A. C., The relaxed energy density for isotropic elastic membranes. IMA Journal of Applied Mathematics, 1986, vol. 36, pp. 85–99.
42 - Pipkin A. C., Relaxed energy densities for large deformations of membranes. IMA Journal of Applied Mathematics, 1994, vol. 52, pp. 297–308.
43 - Ciarlet P. G., Mathematical Elasticity, in 3 vol. Elsevier, 1993.
44 - Friesecke G., James R. D., Muller S., A hierarchy of plate models derived from nonlinear elasticity by Gamma-convergence. Archive for Rational Mechanics and Analysis, 2006, vol. 180, pp. 183–236.
45 - Le Dret H., Raoult A., The nonlinear membrane model as variational limit of nonlinear three-dimensional elasticity. Journal de Mathématiques Pures et Appliquées, 1995, vol. 74, pp. 549–578.
46 - Conti S., Maggi F., Confining thin elastic sheets and folding paper. Archive for Rational Mechanics and Analysis, 2008, vol. 187, pp. 1–48.
47 - Friesecke G., James R. D., Muller S., A theorem on geometric rigidity and the derivation of nonlinear plate theory from three-dimensional elasticity. Communications on Pure and Applied Mathematics, 2002, vol. 55, issue 11, pp. 1461–1506.
48 - Friesecke G., James R. D., Mora M. G., Muller S., Derivation of nonlinear bending theory for shells from three-dimensional nonlinear elasticity by Γ-convergence. Comptes Rendus de l‘Académie des Sciences. Series I Mathematics, 2003, vol. 336, pp. 697–702.
49 - Kirchhoff G., Über das Gleichgeweicht und die Bewegung einer elastischen Scheibe. Journal für die reine und angewandte Mathematik, 1850, vol. 40, pp. 51–88.
50 - Pipkin A. C., Elastic materials with two preferred states. Quarterly Journal of Mechanics and Applied Mathematics, 1991, vol. 44, issue 1, pp. 1–15.
51 - Cesana P., Plucinsky P., Bhattacharya K., Effective behavior of nematic elastomer membranes. Archive for Rational Mechanics and Analysis, 2015, vol. 218, issue 2, pp. 863–905.
52 - Chagnon G., Rebouah M., Favier D., Hyperelastic energy densities for soft biological tissues: a review. Journal of Elasticity, 2015, vol. 120, pp. 129–160.
53 - Raghavan M.L., Vorp D. A., Toward a biomechanical tool to evaluate rupture potential of abdominal aortic aneurysm: Identification of a finite strain constitutive model and evaluation of its applicability. Journal of Biomechanics, 2000, vol. 33, issue 4, pp. 475–482.