## Abstract

Villi are ubiquitous structures in the intestine of all vertebrates, originating from the embryonic development of the epithelial mucosa. Their morphogenesis has similar stages in living organisms but different forming mechanisms. In this work, we model the emergence of the bi-dimensional undulated patterns in the intestinal mucosa from which villi start to elongate. The embryonic mucosa is modelled as a growing thick-walled cylinder, and its mechanical behaviour is described using an hyperelastic constitutive model, which also accounts for the anisotropic characteristics of the reinforcing fibres at the microstructural level. The occurrence of surface undulations is investigated using a linear stability analysis based on the theory of incremental deformations superimposed on a finite deformation. The Stroh formulation of the incremental boundary value problem is derived, and a numerical solution procedure is implemented for calculating the growth thresholds of instability. The numerical results are finally discussed with respect to different growth and materials properties. In conclusion, we demonstrate that the emergence of intestinal villi in embryos is triggered by a differential growth between the mucosa and the mesenchymal tissues. The proposed model quantifies how both the geometrical and the mechanical properties of the mucosa drive the formation of previllous structures in embryos.

## 1. Introduction

The inner surface of the intestinal tissue is covered by a columnar epithelium exhibiting invaginations, known as the crypts of Lieberkuhn, and finger-like projections, called villi [1]. Intestinal villi contain the majority of differentiated absorptive cells, possessing a peculiar morphology which increases the surface area of the intestine as well as its capacity to absorb liquids and nutrients from food. In fact, an increased surface area decreases the average distance travelled by nutrient molecules, thus improving the effectiveness of the diffusion process. Villi are ubiquitous in the intestine of all vertebrates and their morphogenesis follows similar processes in all the organisms. In fact, villi originate from the embryonic development of a soft tissue, the mucosa, that covers the most internal part of organs and cavities exposed to an external environment, as airways, urogenital and gastrointestinal tract. The mucosa is made of three layers: the epithelium, the lamina propria and the muscularis mucosae, as depicted in figure 1. Epithelial cells constitute the epithelial layer and cover the external part of the mucosa. The lamina propria is the layer that is most specific to the anatomical origin. In particular, the lamina propria of the intestinal mucosa includes blood vessels, lymph nodes and connective tissue composed of cells and an extracellular matrix made of ground substances (fluids) and fibres (collagen and elastin). Finally, the lamina muscularis mucosae is a continuous thin sheet of smooth muscle cells.

The morphogenetic processes leading to the formation of the intestinal villi have been experimentally studied in various living organisms, highlighting similar development stages but different initiation mechanisms. In general, three formation stages have been identified: elevation of previllous ridges, delineation of villus bases and outgrowth of the definitive villi [2]. In human embryos, the formation of villi follows the development of mucosal structures, which are characteristic of the anatomical origin of the tissue. In the upper part of the duodenum invaginations are observed during the mucosa growth, in the remainder of the duodenum vacuoles are present, villi are first seen as separate thickening of the mucosa in the upper part of the ileum while they are preceded by circumferential folding in its lower part [3]. In chick embryo samples [4,5], the formation of previllous ridges starts with the delineation of circumferential folds followed by a longitudinal zigzag patterns, emerging on the mucosal wall around 11–15 days after incubation [6]. Between 14 and 17 days, the base of each villus starts to broaden, showing a transition from a zigzag to a regular bi-dimensional undulation pattern, which forms the base structure for villi development. In fact, a cell population starts growing on each crest of the undulated pattern. The mature villi emerge from the elongation process before 25 days after incubation, forming the regular pattern shown in figure 2.

During this morphogenetic process, the growing mucosa is geometrically confined by its surrounding tissues, the muscularis propria and the serosa. It is well established in morpho-elastic theories that geometrical incompatibilities during growth can give rise to residual stresses [7–9]. A number of morpho-elastic models have been proposed to study the mechanical and geometrical aspects of circumferential buckling in tubular organs, considered as isotropic hyperelastic tissues [10–14]. Only recently, Ciarletta & Ben Amar [15] proposed a variational method for studying the segmentation instability in intestinal tissues accounting also for the anisotropic distribution of collagen fibres inside the tissue. Anyway, none of these models considered the simultaneous occurrence of circumferential and longitudinal instabilities.

In this work, we aim to model the formation of the bi-dimensional undulated patterns in the embryonic intestinal mucosa, which characterize the initial stage of villi morphogenesis. In §2, we give a geometrical description of the model, using the multiplicative decomposition of the deformation gradient for describing the volumetric growth of the elastic mucosa. The mechanical behaviour of the mucosa will be defined by the hyperelastic constitutive equations presented in §3, where the anisotropy of the tissue is accounted in terms of the stiffness and the orientation of the collagen fibres distributed under the epithelial layer. In §4, we will derive the equilibrium equations for the basic axial-symmetric solution, and we perform a linear stability analysis for studying the emergence of the bi-dimensional patterns. The Stroh formulation of the incremental elastic equations will be given and the boundary value problem (BVP) will be solved using a particular numerical scheme. Results relative to different growth and materials properties will be presented in §5, paying attention to the links between elastic parameters, growth and instability thresholds.

## 2. Definition of the geometrical model

In the following, we refer to (*R*, * Θ*,

*Z*) as the cylindrical coordinates in the configuration with orthonormal basis vectors (

**E**

*,*

_{R}**E**

*,*

_{Θ}**E**

*) and to (*

_{Z}*r*,

*,*

*θ**z*) as the cylindrical coordinates in the actual configuration with orthonormal basis vectors (

**e**

*,*

_{r}**e**

*,*

_{θ}**e**

*). The mucosa is described as a thick-walled tube, as sketched in figure 3, whose geometry in the fixed reference configuration is defined by 2.1 where*

_{z}*R*and

_{i}*R*

_{0}are the internal and external radii, respectively, and is the axial length of the cylinder. The morpho-elastic deformation is identified by the mapping

*, defined as 2.2 where*

*χ***X**,

**x**are the position vectors in the reference and actual configurations and , respectively. The external radius is fixed during the entire process (

*r*

_{0}=

*r*(

*R*

_{0}) =

*R*

_{0}), since the outer layers, in particular the muscularis propria and the serosa, are much stiffer than the mucosa and can be considered a rigid constraint during mucosal growth. Moreover, the free annular surfaces cannot slide longitudinally and the mucosal internal surface is free of external traction, because the inner intestinal pressure in embryos is negligible. As proposed by Rodriguez

*et al.*[16], the geometric deformation gradient

**F**= Grad

**x**=

*∂*(

*χ***X**)/

*∂*

**X**, associated with the volumetric growth of the mucosa, can be split into two components, as follows: 2.3 where

**F**

_{g}represents pure volumetric growth and

**F**

_{e}is the elastic deformation tensor. The growth tensor represents a linear map of the tangent space and defines a natural grown state, which is not a configuration because the geometry of the tissue is not necessarily compatible with the constraint, as depicted in figure 4. The elastic deformation restores the compatibility of the tissue deformation, while residual stresses arise. For the sake of simplicity, we consider an homogeneous growth tensor in the form: 2.4 where the operator diag indicates that the growth tensor is diagonal. The growth rates in the radial and longitudinal directions are expressed as

*g*and

_{r}*g*, respectively, and they are assumed constant and positive definite. Finally, as the cells and extracellular matrix constituting the mucosa are prevalently composed by water, the tissue can be considered incompressible. Such an incompressibility constraint reads: 2.5

_{z}In §3, we introduce the nonlinear constitutive equations of the material based on the microstructural properties of the mucosal layer.

## 3. Constitutive equations

Recalling that the characteristic viscoelastic time of soft tissues (seconds) is much smaller than their characteristic time of growth (days), we consider the mucosa having a hyperelastic behaviour, so that a strain energy density function can be defined. The mucosa is modelled as an one-layered, homogeneous, incompressible tissue, composed of a cross-ply continuous distribution of collagen and elastin fibres (anisotropic component), immersed into a homogeneous ground substance (isotropic component). We assume that collagen and fibres are distributed along the two principal directions **m**_{α} and **m**_{−α}, defined as
3.1
where * α* is the cross-ply fibre angle with respect to the longitudinal direction

**e**

*. Hence, neglecting the mutual interaction between fibres and ground substance, we can use the additive decomposition proposed by Holzapfel & Ogden [17], and we can express the strain energy function as a sum of two terms 3.2 where*

_{z}

*Ψ*_{Iso}is the isotropic component, a scalar function of the right Cauchy–Green tensor

**C**

*; and*

_{e}

*Ψ*_{Aniso}is the anisotropic component, which also depends on the orientation of the fibres through the vectors

**m**

_{±α}; and

*p*is a Lagrange multiplier for enforcing incompressibility. Let us assume a Neo-Hookean behaviour of the ground substance, so that the isotropic contribution

*Ψ*_{Iso}(

**C**

*) in equation (3.2) reads: 3.3 where*

_{e}*is the shear modulus,*

*μ**are the principal stretches, with and*

*λ*_{l}*I*

_{1}= tr[

**C**

*] is the first invariant of*

_{e}**C**

*. Let us now define the structural tensors*

_{e}**M**

_{±α}=

**m**

*⊗*

_{±α}**m**

*, where the symbol ⊗ indicates the dyadic product between two vectors ((*

_{±α}**a**⊗

**b**)

*=*

_{kj}*a*). The anisotropic strain energy function

_{k}b_{j}

*Ψ*_{Aniso}(

**C**

*,*

_{e}**m**

_{±α}) can be defined as [18] 3.4 where

*k*

_{1}is the anisotropic stiffness of the material, and

*I*

_{4}is a structural pseudo-invariant, which reads: 3.5 where the symbol indicates the so-called saturation operator between two tensors (

**A**:

**B**= tr[

**A**

^{T}

**B**] =

*A*, and

_{kj}B_{kj}

*λ*_{±α}=

*|*

**Fm**

_{±α}

*|*are the fibre stretches along the directions ±

*. We also note that the strain energy in equation (3.4) is polyconvex and physically consistent both with the compression and the extension of fibres. From constitutive arguments, the nominal stress tensor*

*α***S**of the mucosa can be expressed as 3.6 while the Cauchy stress tensor

*reads: 3.7 where*

*σ*

*ψ*_{1}=

*∂*/

*Ψ**∂I*

_{1},

*ψ*_{4}=

*∂*/

*Ψ**∂I*

_{4}and with

**B**

*indicating the left Cauchy–Green tensor. In the absence of body forces, the equilibrium equations for the mucosa read: 3.8 where Div and div are the divergence operators in the reference and actual configurations, respectively. Finally, we can define the following boundary conditions: 3.9 where*

_{e}**N**and

**n**are the outer normal units in the reference and actual configurations, respectively. In §4, we derive the basic axial-symmetric solution of equations (3.8) and (3.9) and we perform its linear stability analysis.

## 4. Axial-symmetric solution and linear stability analysis

Growth instabilities are usually investigated using the method of incremental deformations superposed on finite deformations. As the basis of this method, the fundamental idea is to consider an infinitesimal perturbation of a basic elastic solution in order to perform a linear stability analysis. Let us first solve the equilibrium equations in equations (3.8) and (3.9) for calculating the basic axial-symmetric solution of the associated BVP, as well as the residual stress distribution inside the tissue. Second, we define an incremental deformation on the finite solution, and we perform a bifurcation analysis for solving the associated BVP and identifying the instability thresholds.

### 4.1. Basic axial-symmetric solution

Let us now derive the axial-symmetric solution for the growing mucosa, taking into account the residual stresses arising from the geometrical constraint. Considering the incompressibility constraint in equation (2.5), the spatial components of the axial-symmetric deformation read:
4.1
where is imposed by the geometrical constraint at *R*_{0}. According to equation (4.1), the principal stretches * λ_{r}*,

*and*

*λ*_{θ}*are given by 4.2 where the comma denotes the partial differentiation. Accordingly, the basic elastic deformation gradient rewrites 4.3 so that the left Cauchy–Green tensor*

*λ*_{z}**B**

*is defined as Since we are considering an incompressible, hyperelastic solid, the principal components of the Cauchy stress*

_{e}*from equation (3.7) read: 4.4 where*

*σ**Σ*

*,*

_{k}*k*= {

*r*,

*,*

*θ**z*} are the diagonal components of

**in equation (3.7). Considering the axial-symmetric solution in equation (4.1), the only non-trivial equilibrium equation of equation (3.8) is given by 4.5 and the only unknown is the Lagrange multiplier**

*Σ**p*. Imposing the boundary condition in equation (3.9), ensuring no stress on the internal surface,

*p*can be obtained solving equation (4.5) as 4.6 Substituting equation (4.6) in equation (4.4), we obtain the spatial distribution of the residual stresses inside the mucosa as a function of the growth rates. In the following, we will perform a linear stability analysis introducing an infinitesimal perturbation on the axial-symmetric solution.

### 4.2. Incremental deformation

Let us consider an infinitesimal perturbation of the basic grown position **x**^{(0)} = *χ*^{(0)}(**X**), defined as
4.7
where so that the term *χ*^{(1)}(**x**^{(0)}) can be regarded as a first-order incremental displacement with respect to the axial-symmetric configuration.

Considering the spatial displacement gradient **Γ** = grad

*χ*^{(1)}(

**x**

^{(0)}) associated with the incremental deformation, the perturbed elastic deformation gradient reads: 4.8 where is the increment of the basic elastic deformation gradient as depicted in figure 5. Now, we want to rewrite the equilibrium equations in the perturbed configuration The nominal stress in equation (4.8) is given after perturbation by 4.9 where the zero-order term

**S**

^{(0)}is the nominal stress calculated for the axial-symmetric configuration. Let be the associated increment of the nominal stress tensor . By differentiating the constitutive equation in equation (3.6), it reads: 4.10

In equation (4.10), the term is the increment in *p* and is the fourth-order tensor of the elastic moduli [19], defined as
4.11
Since the zero-order term in equation (4.9) satisfies the first of equation (3.8), the incremental equilibrium equations rewrite
4.12

Now, let us define the push-forward of in the perturbed configuration as 4.13

where is the push-forward of also known as the tensor of instantaneous moduli 4.14

From equations (4.3) and (4.14), the non-zero components of are
4.15
where * Ψ_{k}* =

*∂*/

*Ψ**∂*,

*λ*_{k}*=*

*Ψ*_{kj}*∂*

^{2}

*/*

*Ψ**∂*and with

*λ*_{k}∂*λ*_{j}*k*and

*j*running over

*r*,

*and*

*θ**z*. Using the well-known properties of Piola transformations, the equilibrium equations in equation (4.12) can be written in the spatial form as 4.16 where is given by equation (4.13). Introducing the stress components in equation (4.16), we can obtain the three incremental equilibrium equations as 4.17

Let us now assume the incremental deformation in the following form: *χ*^{(1)}(**x**^{(0)}):
4.18
where (*u*, *v*, *w*) are scalar functions representing the incremental displacements. According to equation (4.18), the displacement gradient **Γ** reads:
4.19
Since we are dealing with an incompressible tissue, the incremental incompressibility condition reads:
4.20

Equations (4.16) and (4.20) define a system of four partial differential equations, where the four unknowns are the three incremental displacements in equation (4.18) and the increment *q* of the Lagrange multiplier. Once derived the incremental equilibrium equations, we need to impose the necessary boundary conditions for solving the system of equations (4.16) and (4.20) and defining an adjacent equilibrium to the axial-symmetric deformation.

### 4.3. Stroh formulation of the boundary value problem and numerical solution

In order to solve equations (4.16) and (4.20), we transform the system of partial differential equations into a system of ordinary differential equations of the first order. This strategy is also known as the Stroh formulation of the elastic problem [20]. First of all, let us consider incremental displacements in the form:
4.21
where *m* (respectively, *k _{z}* = 2

*n*) is the circumferential (respectively, longitudinal) mode of the perturbation, with

*π*/L*m*and

*n*positive integers, and ‘i’ is the imaginary unit. Using the perturbation defined in equation (4.21), the deformed mucosa is depicted in figure 6, showing the characteristic bi-dimensional undulated pattern at the inner surface emerging at the initial stages of intestinal villi formation.

For physical compatibility, the incremental stress components have the similar form:
4.22
where *S*_{0kj} is a function of *r*, with indices {*k*, *j*} running over {*r*, * θ*,

*z*}. Substituting equations (4.15), (4.21) and (4.22) in the incremental constitutive equations for given by equation (4.13), we obtain 4.23

Similarly, from the constitutive equations for the incremental components and we can explicate *V′* and *W’* as
4.24
and
4.25
Furthermore, the incremental incompressibility condition in equation (4.20) yields
4.26
Now, we first replace the components in equation (4.17) using equation (4.22) and setting *s*_{0kj} = (*irS*_{0kj}) for (*k*, *j*) = {*r*, * θ*,

*z*} [21]. Using equations (4.13), (4.15) and (4.21), the incremental equilibrium equations take the following expressions: 4.27 4.28 and 4.29 where prime denotes the derivative respect to the variable

*r*. Let us define the displacement-traction vector

*as 4.30*

*η*Substituting equations (4.23)–(4.26) in equations (4.27)–(4.29), we eliminate the unknown *Q* and we obtain three ordinary differential equations of first order that depend only on (*U*, *V*, *W*, *s*_{0rr}, *s*_{0rθ}, *s*_{0rz}). Accordingly, equations (4.24)–(4.29) can be written in a more compact formulation as follows:
4.31
where the Stroh matrix **G** has the form:
4.32
In particular, the four blocks of **G** read:
4.33
4.34
where is the adjugate (transpose conjugate) of **G**_{1} and
4.35

with
4.36
In order to solve equation (4.31), we must impose the boundary conditions ensuring the vanishing of the perturbation at the external radius *r* = *R*_{0}, being
4.37

Finally, the stress-free boundary conditions at the internal radius *r* = *r _{i}* impose
4.38

In the following, we will present the numerical procedure for solving equation (4.31) with the boundary conditions in equations (4.37) and (4.38).

### 4.4. Numerical resolution procedure

The numerical solution is calculated using the determinant method, as in [15,22]. The aim is to find the values of a control parameter for which there is a solution of equation (4.31) with mixed boundary conditions given by equations (4.37) and (4.38). Since the growth has been considered homogeneous, the growth rates *g _{r}* and

*g*can be considered as the control parameters of the elastic instability for the mucosa. The numerical procedure used to obtain these bifurcation parameters is implemented as follows. We express the solution

_{z}**of equation (4.30), as a linear combination of three scalar functions**

*η*

*η*_{1},

*η*_{2},

*η*_{3}, being 4.39 where

*ν*_{1},

*ν*_{2},

*ν*_{3}are constant coefficients. The scalar functions

*η**,*

_{k}*k*= {1,2,3}, in equation (4.39) are three linearly independent copies of the solution, obtained numerically integrating the system of equation (4.31) between

*r*=

*R*

_{0}and

*r*=

*r*and imposing three linearly independent sets

_{i}

*η*_{0k}of initial conditions, expressed as 4.40 with arbitrary incremental stress components (

*s*

_{0rr})

*, (*

_{k}*s*

_{0rθ})

*, (*

_{k}*s*

_{0rz})

*at*

_{k}*r*=

*R*

_{0}. Then, we must impose that the solution in equation (4.39) satisfies the other boundary conditions in equation (4.38) at

*r*=

*r*, which rewrites: 4.41 where (

_{i}*s*

_{0jl})

*,*

_{k}*jl*= {

*r*,

*,*

*θ**z*},

*k*= {1, 2, 3} are the incremental stresses numerically calculated at

*r*=

*r*from the initial value

_{i}

*η*_{0k}. As shown in figure 7, we calculate the bifurcation threshold

*g*(

_{τ}*H*),

*= {*

*τ**r*,

*z*} with the help of two cycles of iteration: we first iterate on the aspect ratio

*H*=

*R*

_{0}/

*R*, then, in a second cycle, we iterate on the bifurcation parameter until the stop condition in equation (4.41) is reached.

_{i}## 5. Results

The numerical results obtained from the linear stability analysis of the growing intestinal mucosa are presented in the following. First, we consider the mucosa as an isotropic material, with the aim of investigating the effect of the volumetric growth on the onset of instability. The volume increase is considered resulting from both isotropic (*g _{r}* =

*g*) and anisotropic (

_{z}*g*≠

_{r}*g*) growth processes. The curves of marginal stability depict the growth thresholds at which a bifurcation occurs as a function of the aspect ratio

_{z}*H*and are shown for different perturbation modes. Second, the mucosa is considered as a fibre-reinforced tissue according to equation (3.2). The numerical results are shown for investigating the effect of the material anisotropy on the bifurcation thresholds for

*g*and

_{r}*g*, investigating the effects of the stiffness and the orientation of the reinforcing fibres on the onset of instability.

_{z}### 5.1. Isotropic behaviour of the mucosa

Let us first consider the mucosa as an isotropic material (i.e. setting *k*_{1} = 0 in equation (3.2)). The marginal stability curves obtained assuming an isotropic growth process (*g _{r}* =

*g*) with

_{z}*m*=

*k*are plotted in figure 8 within the range of 1 <

_{z}*H*≤ 2. An increase in the value of the perturbation modes results in a decrease in the growth thresholds, highlighting the occurrence of a surface instability at very short wavelengths.

The same instability mechanism occurs when considering anisotropic growth processes. The marginal stability curves are plotted in figure 9, referring to a volume increase completely due to a radial (*g _{z}* = 1, figure 9

*a*) or to a longitudinal (

*g*= 1, figure 9

_{r}*b*) growth process.

In figure 10, the marginal stability curves are shown for different circumferential perturbation modes *m*, at a fixed perturbation mode in the longitudinal direction, *k _{z}* = 10. In the same way, fixing the circumferential perturbation mode at

*m*= 10, the instability thresholds are depicted in figure 11, for different longitudinal modes

*k*.

_{z}These results confirm the occurrence of a surface instability mechanism: the growth thresholds for high values of *m*, *k _{z}* collapse to a single master curve in the case of thick tissues (

*H*> 1.5), while a large variability on the perturbation mode appears as the aspect ratio decreases. Even if the instability is predicted for (

*m*,

*k*) →

_{z}*∞*, we recall that the biological system will select a finite wavelength because of the existence of boundary energies penalizing the increase in the surface area of the mucosa. Although such correction of the wavelength can be calculated following the method proposed by Ben Amar & Ciarletta [23], it will be neglected here for the sake of simplicity.

Finally, it is useful to compare the instability thresholds for the three different growth processes investigated in terms of the total volume increase, given by In figure 12, we depict the marginal stability curves for *J*_{g} in the cases *g _{r}* =

*g*,

_{z}*g*= 1 and

_{r}*g*= 1, predicting that the bi-dimensional surface undulations occur first when the mucosal growth is uniquely longitudinal.

_{z}### 5.2. Anisotropic behaviour of the mucosa

In this paragraph, we take into account the material anisotropy of the mucosa, and the numerical results are calculated for different cross-ply fibre angles * α* (figure 13) and for different material anisotropy ratios

*k*

_{1}/

*(figure 14).*

*μ*The curves plotted in figure 13 show how increasing the cross-ply angle * α* increases the instability thresholds, both for the radial and longitudinal growth processes. Hence, a surface instability first occurs when the fibre orientation angle is

*= 0.*

*α*Moreover, the curves of marginal stability in figure 14 show that the presence of a material anisotropy is a stabilizing factor, as increasing *k*_{1}/* μ* increases the growth thresholds of the instability. Finally, we underline that the results in figure 13 (respectively, figure 14) are shown for 1 <

*H*< 1.25 (respectively, 1.5) as the equation system became stiff for thicker tissues. Although other numerical techniques exist to overcome this difficulty (see [24], for a review), we preferred to avoid further complications here, because the mechanical effects of fibre reinforcements on instability are clearly identifiable in the observed range.

## 6. Discussion and conclusion

In this paper, we have proposed a morpho-elastic model of the embryonic mucosa in order to investigate the occurrence of undulated bi-dimensional patterns during the initial stages of villi morphogenesis. In §2, the intestinal mucosa has been modelled as a thick-walled cylinder with an outer spatial confinement. Its deformation gradient has been decomposed into a homogeneous growth component and an incompressible elastic deformation, in order to account for the spatial distribution of residual stresses arising during the constrained growth process. The mechanical behaviour of the mucosa has been described using the hyperelastic model presented in §3. Although similar geometrical approaches have been used to study the occurrence of folding instability in tubular organs [11,12,25], a novelty of our model is to consider the material anisotropy of the tissue. In fact, the nonlinear strain energy function includes a polyconvex energy term that depends on both the orientation and the stiffness of the collagen and elastin fibres lying under the epithelial layer. In §4, we have derived the equilibrium equations for the basic axial-symmetric solution of the growing mucosa, and we have performed a linear stability analysis using the theory of small deformations superposed on the axial-symmetric solution. Another novel aspect is the definition of a bi-dimensional perturbation for describing the undulated morphology at the free surface of the mucosa, as depicted in figure 6. The resulting incremental BVP has been derived using the Stroh formulation, and it has been solved with the help of a particular numerical scheme. The numerical results have been presented in §5. The growth rate thresholds for the onset of instability have been depicted as a function of the aspect ratio *H* =*R*_{0}/*R _{i}*, considering isotropic and anisotropic growth processes in figures 8 and 9, respectively. The marginal stability curves show that a short-wavelength bi-dimensional undulation occurs. An increase in the perturbation modes results in decreasing instability thresholds, showing the occurrence of a surface instability on the internal surface of the mucosal wall, as depicted in figures 10 and 11. The calculated volume increase thresholds are smaller when considering an anisotropic growth process, with the most unstable scenario being the mucosa growing only along the longitudinal direction, as shown in figure 12. In §5.2, we investigated the role played by the material anisotropy of the mucosa on the onset of the surface instability. The results in figures 13 and 14 show that an increase in both the cross-ply angle and the stiffness of the reinforcing fibres provokes an increase in the growth instability thresholds, highlighting that the material anisotropy is a stabilizing effect.

Finally, we aim to compare the growth thresholds obtained here for the occurrence of undulated bi-dimensional patterns with those reported in [25] for prismatic deformations, representing both circumferential folding and longitudinal segmentation. This comparison shows how the instability thresholds for the bi-dimensional perturbations are smaller than those required for folding and segmentation, both for isotropic (figure 15*a*) and anisotropic (figure 15*b*) growth processes. In this way, we also prove that a bi-dimensional undulation always occurs first for a thin mucosa (i.e. *H* < 1.15), while for thicker tubes the growth threshold is roughly the same as the one reported for circumferential folding. Using the experimental curves reported by Sbarbati [26], we have calculated the aspect ratio of the mucosa in mouse embryos within the range *H* = 1.75–1.93 (duodenum) and *H* = 1.51–1.72 (large intestine), while the external radius is about *R*_{0} = 20−40 μm between 12 and 16 days after incubation. The predictions of our model using such values of aspect ratio are consistent with the experimental observations of circumferential folding preceding the villi elongation mouse embryos.

Moreover, we can make a comparison of our theoretical predictions by using the experimental data of Sbarbati & Strakee [27], who reported the volume of the mouse embryonic mucosa within the period 11–19 days after incubation. Considering an initial aspect ratio *R*_{0}/*R _{i}* = 1.7 from the geometrical data in [26], the critical volume increase for instability is about 1.3 according with the results sketched in figure 15. This critical value is reached during experiments between days 13 and 15, which corresponds to the period when the undulated pattern is first detected, therefore validating our theoretical predictions. It is worth mentioning that even if the linear stability analysis fixes the threshold and the wavenumber of the unstable wrinkling modes, the emerging three-dimensional pattern at the tissue surface will be driven by nonlinear effects. This aspect is beyond the scope of this paper, but future work will address the study of the nonlinear development of villi morphology by minimizing the incremental strain energy inside the mucosa at orders in

*higher than two. Weakly nonlinear elastic theories [28] and numerical post-buckling techniques [29] will be used for this purpose.*

*ɛ*In conclusion, the proposed morpho-elastic model confirms that the emergence of intestinal villi in embryos is triggered by a differential volumetric growth between the epithelial mucosa and the mesenchymal tissues. We have demonstrated that villi morphogenesis can start directly from a bi-dimensional undulation of the mucosa or indirectly after a circumferential folding, later followed by a zigzag pattern, as reported by experimental observations in embryos of different species. Our theoretical analysis predicts that the selected previllous structure on the mucosa is mainly driven by the initial aspect ratio of the tubular tissue, with thinner tubes not requiring any preceding mucosal folding. The proposed morpho-elastic model highlights that both the geometrical and the mechanical properties of the mucosa strongly influence the formation of previllous structures in embryos, providing useful suggestions for interpreting the dynamics of villi morphogenesis in living organisms. The comprehension of the interplay between growth and elasticity in the development of embryonic mucosa might also be important for understanding altered villi structures driven by pathological conditions, such as the marked flattening of the mucosal surface in coeliac disease [30].

## Acknowledgements

Partial funding by the European Community grant FP7 ERG-256605 is gratefully acknowledged.

- Received February 2, 2013.
- Accepted February 21, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.