Elephant trunk wrinkles: A mathematical model of function and form
Yang Liu Alain Goriely L. Angela Mihai





Elephant trunk wrinkles:
 A mathematical model of function and form
    Yang Liu[Mathematical Institute, University of Oxford, Woodstock Road, Oxford, OX2 6GG, UK, Department of Mechanics, School of Mechanical Engineering, Tianjin University, Tianjin 300354, China]
  Alain Goriely[Mathematical Institute, University of Oxford, Woodstock Road, Oxford, OX2 6GG, UK]
  L. Angela Mihai[Corresponding author: L.A. Mihai, Email: , School of Mathematics, Cardiff University, Senghennydd Road, Cardiff, CF24 4AG, UK]

    
==========================================================================================================================================================================================================================================================================================================================================================================================================================================================


§ ABSTRACT 
12pt
A remarkable feature of the elephant trunk is the pronounced wrinkling that enables its great flexibility. Here, we devise a general mathematical model that accounts for characteristic skin wrinkles formed during morphogenesis in elephant trunk. Using physically realistic parameters and operating within the theoretical framework of nonlinear morphoelasticity, we elucidate analytically and numerically the effect of skin thickness, relative stiffness and differential growth on the physiological pattern of transverse wrinkles distributed along the trunk. We conclude that, since the skin and muscle components have similar material properties, geometric parameters, such as curvature, play important roles. In particular, our model predicts that, in the proximal region close to the skull, where curvature is lower, fewer wrinkles form and sooner than in the distal narrower region where more wrinkles develop. Similarly, less wrinkling is found on the ventral side, which is flatter, compared to the dorsal side. In summary, the mechanical compatibility between the skin and the muscle enables them to grow seamlessly, while the wrinkled skin acts as a protective barrier that is both thicker and more flexible than the unwrinkled skin. 
6pt
6pt
Keywords: hyperelastic solids, nonlinear deformation, instabilities, wrinkling, bilayer systems, mathematical modeling.
6pt
Mathematics Subject Classification: 74B20, 74G10, 74G60, 9210.
 




“Elephants have wrinkles, wrinkles, wrinkles,

Elephants have wrinkles, wrinkles everywhere.

On their trunks, on their ears, on their knees,

On their toes, no one knows, no one knows

Why-I-I-I-Ih!”– Nursery Rhyme









§ INTRODUCTION


The elephant trunk (Latin proboscis, Ancient Greek πρ oγ oσκíç or proboskís) is an iconic, highly versatile limb containing approximately 100,000 radial and lateral muscles fibers, connected to the elephant's head by an opening in the skull and controlled by a proboscis nerve <cit.>. Due to its exceptional characteristics, it has attracted much interest from the research community motivated by biomimicking its properties and movements <cit.>. 





A salient feature contributing to elephant trunk's agility is its rich pattern of wrinkles (see Figure <ref>), which form and develop throughout elephant's life <cit.>. During fetal growth, transverse trunk wrinkles and furrows appear gradually <cit.>. In adult elephants, the number of transverse wrinkles differs longitudinally, with more wrinkles formed distally (at the tip) than proximally (near the skull), and more dorsal (top) than ventral wrinkles. The number of wrinkles further changes with trunk-lateralisation. MicroCT-imaging indicates that the outer trunk skin has a relatively constant thickness, while the inner skin is thinner within folds than between folds. Thus the generation of wrinkles in elephant trunk is primarily a result of differential growth appearing during development <cit.>.

Wrinkling instability is a ubiquitous mechanism involving deformations across the scales. A general theory of small strain superposed on large strain deformations for homogeneous isotropic hyperelastic materials was developed by <cit.>. Particular cases of infinitesimal deformations superposed on finite extension or compression were analyzed in  <cit.>. The stability of a solid circular cylinder subject to finite extension and torsion was treated in <cit.>. For tubes of Mooney-Rivlin material with arbitrary length and thickness, by applying the Stroh formalism  <cit.>, an explicit formulation for the Euler buckling load with its first nonlinear corrections was derived in <cit.>. A similar approach was employed in <cit.> to calculate the nonlinear buckling load of a compressible elastic cylinder. Surface wrinkles formed by straightening of a circular sector or by bending of a cylindrical sector were examined in <cit.> and <cit.>, respectively. Wrinkling due to bending of an inflated cylindrical shell was discussed in <cit.>. An extension of Flügge's formalism <cit.> applied to the buckling of thin-walled cylinders of nonlinear elastic material was presented in <cit.>. The instability of a thick-walled hyperelastic tube subject to external pressure and axial compression was considered in <cit.>. Other forms of instability, like necking or bulging, also arise during inflation of elastic tubes. These phenomena were studied extensively by <cit.>. The influence of residual stresses on the stability of circular cylinders, and in particular on stretch-induced localized bulging and necking, was addressed by <cit.>. Post-buckling modes for a core-shell cylindrical system of neo-Hookean material under axial compression with a perfectly bonded interface were investigated experimentally, theoretically and computationally in <cit.>. The influence of geometrical incompatibility on pattern selection of growing bilayer tubes was modelled in <cit.>. A cylinder with shear modulus arbitrarily varying along the radial direction which buckled and wrinkled under axial compression was examined semi-analytically in <cit.>. The onset of wrinkling in an anisotropically growing fibre-reinforced tube was analyzed in <cit.>. Wrinkling patterns in a cylindrical bilayer where only the outer layer grew were explored in <cit.>. Radially distributed wrinkles induced by a thin elastic film growing on a soft elastic substrate were presented in <cit.>.

In this paper, building on the rich methodology for elastic instabilities, we construct a mathematical model that accounts for transverse skin wrinkling in elephant trunk. Working within the theoretical framework of large strain morphoelasticity <cit.>, we address the following key question: What is the effect of the relative thickness, stiffness, and growth of the skin and muscle substrate on wrinkles formation along the elephant trunk? In Section <ref>, we model the trunk as a tubular cylindrical bilayer where the skin forms the thin outer layer, the muscle constitutes the thick inner layer, and a perfect bond exists between these two constituents. To account for large elastic strains, we describe each layer as an incompressible homogeneous isotropic nonlinear hyperelastic material. Assuming axisymmetric deformations, in Section <ref>, we present extensive numerical results for skin wrinkling when some model parameters change while others are fixed. In Section <ref>, we treat analytically and numerically different limiting cases. From our comprehensive analysis, we conclude that relative growth, geometry and material properties, together with loading conditions, compete to generate the characteristic pattern of transverse wrinkles in elephant trunks. While this study is motivated by a specific application scenario, our fundamental results extend further and are relevant to other applications as well.




§ PROBLEM FORMULATION


We model (a segment of) the elephant trunk as a cylindrical tube bilayer composed of two concentric homogeneous isotropic incompressible hyperelastic tubes. The skin forms the outer shell, while the muscle substrate constitutes the inner thicker core. We denote the reference state by ℬ_0, and set R_0, R_1, h_0 and l_0 the inner and outer radii of the core, the uniform radial shell thickness, and the axial length of the cylindrical system, respectively, where R_0<R_1 and h_0≪ R_1-R_0. Assuming that both the shell and the core grow until a deformed state is attained, we designate ℬ as the current configuration where the cylindrical geometry is maintained and the geometrical parameters become r_0, r_1, h, and l, respectively, with r_0<r_1 and h≪ r_1-r_0. We further assume that, at one end, the cylindrical system is fixed (homogeneous Dirichlet boundary condition), while at the other end, the system is constrained elastically (Robin boundary condition <cit.>) by a linear (Hookean) spring of stiffness K (Hooke's constant). Our bilayer system is represented schematically in Figure <ref>.





Within the usual system of cylindrical polar coordinates, we denote by 𝐗=(R,Θ,Z) and 𝐱=(r,θ,z) a material point in the reference and current configuration, respectively, where

    R_0⩽ R⩽ R_1+h_0,     0⩽Θ⩽ 2π,     0⩽ Z⩽ l_0,
    
    r_0⩽ r⩽ r_1+h,     0⩽θ⩽ 2π,     0⩽ z⩽ l.

For both the shell and the core, the deformation gradient 𝐅 from the reference to the current configuration takes the form <cit.>

    𝐅=𝐀𝐆,

where A is the elastic deformation tensor and G is the growth tensor.

Given the strain-energy function W(𝐀) of an incompressible homogeneous isotropic hyperelastic material, the corresponding nominal stress tensor 𝐒 and Cauchy stress tensor σ are, respectively,

    𝐒=J𝐆^-1∂ W/∂𝐀-pJ𝐆^-1𝐀^-1        σ=𝐀∂ W/∂𝐀-p𝐈,

where J=𝐅=𝐆 represents the volume change due to biological growth, p denotes the Lagrange multiplier associated with the incompressibility constraint for the isochoric elastic deformation 
A=1, and 𝐈 is the second-order identity tensor.

In the current configuration, the equilibrium equation for a quasi-static problem without body forces takes the form

    ÷ σ=0.

For the axisymmetric deformation considered here, the only non-trivial equation is

    dσ_11/d r+σ_11-σ_22/r=0,

where the indices 1,2,3, correspond to r-, θ-, z-directions, respectively.

In the subsequent analysis, we add the subscripts s and m to distinguish between physical quantities associated with the skin shell and the muscle core, respectively. For example, W_s represents the strain-energy function for the skin and W_m for the muscle. Physical quantities that are valid for both these components will be written without their subscripts. 

The primary deformation from ℬ_0 to ℬ is governed by the equations

    r=r(R),    θ=Θ,     z=λ_zZ,

where the longitudinal stretch ratio λ_z is constant. The deformation gradient then simplifies as follows,

    𝐅=d r/d R𝐞_r⊗𝐞_r+r/R𝐞_θ⊗𝐞_θ+λ_z𝐞_z⊗𝐞_z,

where {𝐞_r,𝐞_θ,𝐞_z} is the usual orthonormal basis for the coordinate system.

Denoting by {λ_i}_i=1, 2, 3 the principal stretch ratios in the radial, azimuthal and axial direction, respectively, by equation (<ref>), we obtain

    λ_1=dr/dR,    λ_2=r/R,    λ_3=λ_z.


We further consider the following growth tensors for the skin and muscle, respectively:

    {   𝐆_s = 𝐞_r⊗𝐞_r+𝐞_θ⊗𝐞_θ+γ g𝐞_z⊗𝐞_z,
       𝐆_m= 𝐞_r⊗𝐞_r+𝐞_θ⊗𝐞_θ+g𝐞_z⊗𝐞_z,
    .

where g≥1 is their common growth factor and γ≥1 is the differential (or relative) growth factor between the skin (outer layer) and the muscle (inner layer) in the axial (length) direction. Hence,

    𝐆_s=𝐆_γ𝐆_m,
            𝐆_γ = 𝐞_r⊗𝐞_r+𝐞_θ⊗𝐞_θ+γ𝐞_z⊗𝐞_z.

From equation (<ref>), we deduce that  

    {   𝐀_s = α_s1𝐞_r⊗𝐞_r+α_s2𝐞_θ⊗𝐞_θ+α_s3𝐞_z⊗𝐞_z,
       𝐀_m= α_m1𝐞_r⊗𝐞_r+α_m1𝐞_θ⊗𝐞_θ+α_m3𝐞_z⊗𝐞_z,
    .

and, by equation (<ref>), we have

    α_m3=γα_s3=λ_zg^-1.


From the incompressibility conditions A_s=α_s1α_s2α_s3=1 and A_m=α_m1α_m2α_m3, together with the multiplicative decomposition (<ref>), we infer that 

    {   rdr/dR = λ_z^-1γ gR,   r_1⩽ r⩽ r_1+h,
       rdr/dR = λ_z^-1gR,      r_0⩽ r⩽ r_1.
    .

By integration, we obtain

    {   r^2 = λ_z^-1γ g(R^2-R_1^2)+r_1^2,     r_1⩽ r⩽ r_1+h,
       r^2 = λ_z^-1g(R^2-R_0^2)+r_0^2,    r_0⩽ r⩽ r_1.
        .

The geometry of the system in the current configuration ℬ is described by

    {   r_1 = [λ_z^-1g(R_1^2-R_0^2)+r_0^2]^1/2,
       h = {λ_z^-1gγ[(R_1+h_0)^2-R_1^2]+r_1^2}^1/2-r_1,
       l = λ_z l_0,
    .

where the deformed inner radius r_0 remains to be determined.

Expressing the strain-energy function in terms of the three principal elastic stretches as W(α_1,α_2,α_3), the non-zero components of Cauchy stress tensor take the following form,

    {   σ_sii = α_siW_si-p_s,
       σ_mii= α_miW_mi-p_m,
    .

where W_si=∂ W_s/∂α_i and W_mi=∂ W_m/∂α_i, i=1,2,3 (there is no summation for repeated indices in the above equations). Assuming that both the outer and inner surfaces are free and the interface between the shell and the core remains continuous (perfectly bonded), we have

    σ_m11|_r=r_0=0,    σ_s11|_r=r_1+h=0,    (σ_s11-σ_m11)|_r=r_1=0.


Since α_1α_2α_3=1, we can define the following function depending only on two variables, 

    w(α,α_3)=W(α_1,α_2,α_3),

where α_2=α and α_1=α^-1α_3^-1. Applying the chain rule, we obtain

    w_1=∂ w/∂α=α^-1(α_2W_2-α_1W_1),     
    w_2=∂ w/∂α_3=α_3^-1(α_3W_3-α_1W_1).


From (<ref>)-(<ref>), it follows that

    {   σ_s11 = ∫_α_h^α_sw_s1/1-α^2α_s3dα,
       σ_m11= ∫_α_r_1^α_mw_s1/1-α^2α_m3dα+∫_α_h^α_r_1w_s1/1-α^2α_s3dα,
       0  = ∫_α_r_1^α_r_0w_s1/1-α^2α_m3dα+∫_α_h^α_r_1w_s1/1-α^2α_s3dα,
    .

where 

    α_r_0=r_0/R_0,    α_r_1=r_1/R_1,    α_h=r_1+h/R_1+h_0.

In view of (<ref>), the above elastic stretches are connected by 

    {   α_r_1^2=α_r_0^2R_0^2+g(R_1^2-R_0^2)R_1^2,
       α_h^2=α_r_0^2R_0^2+g(R_1^2-R_0^2)+gγ[(R_1+h_0)^2-R_1^2](R_1+h_0)^2.
        .


The associated Lagrange multipliers can be determined from the following identities:

    {   p_s=α_s1W_s1-σ_s11,
       p_m=α_m1W_m1-σ_m11.
    .


The resultant axial force is equal to

    N=    2π(∫_r_0^r_1σ_m33rdr+∫_r_1^r_1+hσ_m33rdr)
    =    π∫_r_0^r_1[2(α_m3W_m3-α_m1W_m1)-(α_m2W_m2-α_m1W_m1)]rdr
       +π∫_r_1^r_1+h[2(α_s3W_s3-α_s1W_s1)-(α_s2W_s2-α_s1W_s1)]rdr
    
    =    π R_0^2(λ_zα_r_0^2-g)∫_α_r_1^α_r_02α_m1w_m2-α w_m1/(1-α^2α_m3)(g-α^2α_m3)αdα
       +π R_1^2(λ_zα_r_1^2-gγ)∫_α_r_h^α_r_12α_s1w_s2-α w_s1/(1-α^2α_s3)(gγ-α^2α_s3)αdα.


Since the axial displacement is restricted at one end while the other end is attached to a spring of stiffness K, we can express the axial force as

    N=K(l-l_0)=K(λ_z-1)l_0.

Then, for a given value of g, we can derive the axial extension λ_z and the deformed inner radius r_0 from the system of equations: 

    N(r_0,λ_z)=K(λ_z-1)l_0,    σ_m11(r_0,λ_z)=0,

as illustrated in Figure <ref>.

For definiteness of our problem, we consider the incompressible neo-Hookean-type functions:

    {   W_s=μ_s/2(λ_1^2+λ_2^2+λ_3^2-3),
       W_m=μ_m/2(λ_1^2+λ_2^2+λ_3^2-3),
    .

where μ_s>0 and μ_m>0 are shear moduli at infinitesimal strain. The respective non-dimensionalized axial Cauchy stress components are indicated in Figure <ref>.

It is useful to introduce the following dimensionless quantities:

    R_0=R_0/l_0,    R_1=R_1/l_0,    h_0=h_0/l_0,    ζ=h_0/R_1,
       N=N/πμ_m[(R_1+h_0)^2-R_0^2],    σ=σ/μ_m,    K=K/πμ_m(R_1+h_0+R_0),    β=μ_s/μ_m.

Then, by equation (<ref>), the dimensionless resultant axial force is

    N=-K(λ_z-1)/R_1+h_0-R_0.

Since the skin shell is stiffer than the muscle core, surface wrinkles will be generated at some critical growth. For simplicity, henceforth, we drop the over-hat from our notation. 








§ LINEAR BIFURCATION ANALYSIS


In this section, we derive the critical condition for transverse wrinkling (see Figure <ref>) using the Stroh formalism <cit.>. We then construct a robust numerical scheme using the surface impedance matrix method <cit.> to solve the eigenvalue problem arising from wrinkling instability. 








 §.§ Incremental theory


To obtain the incremental equations for the stability analysis, we denote the perturbed state for the bilayer system by ℬ, with the associated position vector 𝐱. The relation between the position vectors in ℬ and ℬ is

    𝐱(r,z)=𝐱+u,

where u=u_1(r,z)𝐞_r+u_3(r,z)𝐞_z is the incremental displacement when θ is fixed. We focus on the axial instability <cit.> since major wrinkles in elephant trunk tend to develop in the longitudinal direction <cit.>.

The deformation gradient from the reference configuration ℬ_0 to the perturbed configuration ℬ can be expressed as follows,

    𝐅=∂𝐱/∂𝐗=(𝐈+η)𝐅,

where η=𝐮 stands for the displacement gradient given by

    η=[
        [      u_1,1 u_1,2-u_2r      u_1,3;      u_2,1 u_1r+u_2,2      u_2,3;      u_3,1      u_3,2      u_3,3 ]]    =[
        [ u_1,1     0 u_1,3;     0  u_1r     0; u_3,1     0 u_3,3 ]],

with u_i,1=∂ u_1/∂ r, u_i,2=∂ u_1/(r∂θ), u_i,3=∂ u_i/∂ z, i=1,2,3. 

Assuming that the growth tensor 𝐆 is constant, we have

    𝐀𝐆=𝐅=(𝐈+η)𝐅=(𝐈+η)𝐀𝐆,

and therefore,

    𝐀=(𝐈+η)𝐀.


Recalling that the elastic deformation is isochoric, incompressibility implies (𝐈+η)=1, which in its linearized form specialises to 

    trη=u_1,1 +u_1/r+u_3,3 =0,

where `tr' is the trace operator. 

The nominal stress tensor in ℬ takes the form 

    𝐒=J𝐆^-1∂ W/∂𝐀-pJ̃𝐆^-1𝐀^-1,

with p the associated Lagrange multiplier. 

We define the Lagrange multiplier increment 

    p=p-p

and introduce the incremental stress tensor

    χ^T=J^-1(-)^T,

where `T' denotes the transpose operator. 

The incremental equilibrium equation reads

    divχ^T=0.


As the magnitude of each component of η is small, we can expand χ in terms of η, retaining the linear terms, i.e.,

    χ_ij=𝒜_jilkη_kl+pη_ji-pδ_ji+𝒪(|η_ij|^2),   i,j,k,l=1,2,3,

where the summation convention for repeated indices is applied and 𝒜=(𝒜_jilk)_i,j,k,l=1,2,3 denotes the first-order instantaneous modulus tensor with the following nonzero entries:

    {   𝒜_iijj=α_iα_j W_iW_j,
       𝒜_ijij=α_iW_i-α_jW_j/α_i^2-α_j^2α_i^2,   i≠ j and α_i≠α_j,
       𝒜_ijij=𝒜_iiii-𝒜_jjjj+α_iW_i/2,   i≠ j and α_i=α_j,
       𝒜_ijji=𝒜_ijij-α_iW_i,   i≠ j.
        .

Note that there is no summation for repeated indices in the above expressions and the tensor 𝒜 has pairwise symmetry, 𝒜_ijkl=𝒜_klij. 

By (<ref>), the nonzero incremental stress components are

    {   χ_11=𝒜_1111η_11+𝒜_1122η_22+𝒜_1133η_33+p η_11-p, 
       χ_22=𝒜_2211η_11+𝒜_2222η_22+𝒜_2233η_33+p η_22-p,
       χ_33=𝒜_3311η_11+𝒜_3322η_22+𝒜_3333η_33+p η_33-p,
       χ_13=𝒜_3131η_13+𝒜_3113η_31+p η_31 ,
       χ_31=𝒜_1313η_31+𝒜_1331η_13+p η_13.
    .

Then equation (<ref>) reduces to 

    {   (rχ_11)_,1+rχ_13,3-χ_22=0,
       (rχ_31)_,1+rχ_33,3=0.
        .


We emphasize that the above general formulae are valid for both the skin and muscle layers. To specify a layer, the related subscript will be added to indicate its affiliation. 

The corresponding boundary and interface conditions are, respectively: 

    χ_s𝐞_r|_r=r_1+h= 0,    (χ_s-χ_m)𝐞_r|_r=r_1= 0

and

    (u_s1-u_m1)|_r=r_1=0,    (u_s3-u_m3)|_r=r_1=0.


The sliding conditions at the ends are

    𝐞_r·χ_s𝐞_z|_z=0,l=0,    𝐞_r·χ_m𝐞_z|_z=0,l=0.





 §.§ Stroh formulation


Employing the Stroh method <cit.>, we look for periodic solutions of the form:

    {   u_1(r,z)=U(r)cos kz,     u_3(r,z)=W(r)sin kz,
       χ_11(r,z)=T_11(r)cos kz,    χ_31(r,z)=T_31(r)sin kz,
       .

where k represents the wave number in the axial direction. In view of equations (<ref>), (<ref>), denoting the dimensionless wave number by n, the sliding conditions (<ref>) yield

    k=nπ/l,   n=1,2,3,⋯.


We further define a displacement-traction vector

    Γ=[𝐔(r),𝐓(r)]^T,    𝐔(r)=[U(r),W(r)]^T    𝐓(r)=[T_11(r),T_31(r)]^T.

From (<ref>), (<ref>), and (<ref>)_1,5, we obtain 

    Γ'(r)=1/r𝐍(r)Γ(r),

with the prime denoting differentiation with respect to r and the Stroh block matrix taking the form

    𝐍=[[    𝐍_1    𝐍_2;    𝐍_3 -𝐍_1^T ]],

where 𝐍_i (i=1,2,3) are 2×2 sub-matrices, such that 𝐍_2 and 𝐍_3 are symmetric. These sub-matrices can be expressed as follows,

    𝐍_1=[[                 -1                -rk; rk(𝒜_1331+p)𝒜_1313                  0 ]],    𝐍_2=[[       0       0;       0 1𝒜_1313 ]],    𝐍_3=[[ t_11 t_12; t_12 t_22 ]],

where 

    {   t_11=𝒜_1111-2𝒜_1122+𝒜_2222+r^2k^2𝒜_3131-r^2k^2/𝒜_1313(𝒜_1331+p)^2+2p,
       t_12=rk(𝒜_1111-𝒜_1122-𝒜_1133+𝒜_2233+p),
       t_22=r^2k^2(𝒜_1111-2𝒜_1133+𝒜_3333+2p).
    .





 §.§ The surface impedance matrix method


Next, we apply the surface impedance matrix method <cit.>, and introduce the conditional impedance matrix 𝐙(r) which satisfies

    𝐓=𝐙(r)𝐔.

From the relations

    𝐔'=1/r(𝐍_1𝐔+𝐍_2𝐓),    𝐓'=1/r(𝐍_3𝐔-𝐍_1^T𝐓),

we obtain the Riccati equation

    𝐙'=1/r(𝐍_3-𝐍_1^T𝐙-𝐙𝐍_1-𝐙𝐍_2𝐙).


We then apply the above general expressions to the bilayer system. Fur the skin layer, we use the boundary condition 𝐙(r_1+h)=0 determined from (<ref>) to integrate (<ref>) from r_1+h to r_1 and find 𝐙_s(r_1). Applying the same procedure to the muscle layer yields 𝐙_m(r_1). The displacement and traction continuities at interface results in (𝐙_m(r_1)-𝐙_s(r_1))𝐔_s(r_1)=0. For surface wrinkling, the existence of a non-trivial solution of 𝐔 finally leads to the bifurcation condition

    Φ(β,γ,n,g,R_0,R_1,h_0,K)=0,

where

    Φ(β,γ,n,g,R_0,R_1,h_0,K)=(𝐙_m(r_1)-𝐙_s(r_1)).





 §.§ Numerical examples


In this section, the interplay between different model parameters is captured by numerical examples where some model parameters vary while others are kept fixed. For the elephant trunk, in <cit.>, the Young moduli of skin and muscle tissue are estimated as E_s=3μ_s≤ 1190±120 kPa and E_m=3μ_m≈ 938 kPa, respectively, and the skin is found at most 1.27 times stiffer than the muscle.
Nonetheless, we should mention that those moduli were measured by deforming wrinkled skin samples, first at small strain then at large strains, and therefore, the large strain  measurements would be closer to those for the unwrinkled skin, which are not available. Additionally, these quantities will be different in the unborn elephant, where wrinkles first form <cit.>, than in calf and in adult elephants where wrinkles continue to evolve. Further measurements are needed. Meanwhile, our relevant plots in this paper are those where the modulus ratio is β=1.27. Notwithstanding these estimates, our modelling approach is valid more generally, and we include a wide range of examples which are valuable in their own right and can be useful to this and other applications as well.





In Figure <ref>, we show bifurcation curves identified from the solutions of (<ref>) for different parameter values. These plots suggest that the relative growth factor γ decreases as a function of the wave number n when the pre-growth factor g increases and also when the relative modulus β increases. For each curve, the local minimum gives the first bifurcation point with the critical relative growth γ_cr and the critical wave number n_cr. In particular, when β=1.27, the relative growth factor γ∼ 2 appears to be optimal.





For β  large, as illustrated in Figure <ref> shows bifurcation curves,  all closed curves are shrinking for increasing β, indicating a higher modulus ratio delays the critical relative growth γ_cr. Additionally, the bifurcation condition (<ref>) has no solution for some ranges of n values. This feature is quite different from the previously known result that a stiffer film triggers an earlier instability.  Another inference is that, for a large enough β, the bifurcation curve may vanish. To study this feature, we plot in Figure <ref> the function Φ, given by equation (<ref>), with the differential growth factor γ increasing while n=12 and the other parameters as in Figure <ref>. The dots where the black arrow crosses the curves in Figure <ref> highlight the minimum points, not the bifurcation points (note that the bifurcation point satisfies the bifurcation condition Φ=0). Meanwhile, the arrow indicates that  all curves increase with increasing β, until a critical value  β_max^(n=12), where the function Φ is tangential to the line Φ=0. When β exceeds this value,  the bifurcation condition has no solution.  In general, for  n=n_cr, the critical value is β_max, which gives the global maximum value of β for generating wrinkles. This optimal value can be identified by solving the simultaneous equations:

    Φ=0,    ∂Φ/∂ n=0,    ∂Φ/∂γ=0.

where the bifurcation function Φ is described by  equation (<ref>). For the elephant trunk, the relevant results are those where β is close to 1.





Figure <ref> displays the monotonically increasing β_max as a function of K. When K→0, corresponding to a free boundary, we find β_max≈0.6466. With a free end, if the skin modulus is higher or comparable to that for the muscle, so that β⩾1, surface wrinkles are impossible, even when compression is generated by differential growth, i.e., when γ>1. The underlying mechanism, in this case, is that growth happens through free trunk elongation. In the other limiting case, where K→∞, the spring constraint becomes a fixed end (Dirichlet boundary) condition, which is usually adopted in the context of growth-induced wrinkling  <cit.>. In this scenario, regardless of the value of β, surface wrinkles will form at a critical growth.





To investigate the effect of different parameters on the critical state, we set β=1.27. Figure <ref> shows how the critical relative growth factor and wave number vary with respect to the skin thickness h_0. It suggests that a thicker skin always retards surface wrinkling and, at the same time, reduces the number of wrinkles. This is consistent with other reported results for growth-induced wrinkling <cit.>. Since R_1=0.25, a thicker skin layer will give rise to a lower curvature. As a result, we have that curvature increase will cause an earlier instability if the muscle thickness is fixed. 





Next, we comment on the plots in Figure <ref> showing the dependence of γ_cr and n_cr on the radius of the bilayer tube R_1+h_0. In particular, we keep the thickness ratio ζ=h_0/R_1 constant, meaning that both h_0 and R_1 will increase when R_1+h_0 increases. It can be seen from Figure <ref> that the critical differential growth γ_cr is non-monotonic but varies only in a small interval [2.099,2.199]. Meanwhile, the critical wave number n_cr drops uniformly. This implies that, when the thickness ratio ζ is fixed, the critical differential growth does not alter much as the bilayer tube becomes thicker or thinner. However, the wave number changes rapidly. This may explain why there are more wrinkles in the distal narrower part of the elephant trunk than in the proximal wider part near the skull (see Figure <ref>).

For planar film/substrate bilayers, surface wrinkles are usually not affected by the thickness of the substrate which can be seen as a half-space <cit.>. To further check whether this is true for curved structures, we display the variations of γ_cr and n_cr as functions of R_1 by fixing all other parameters. As expected, R_1 alone has no effect on the wrinkled pattern but it does affect the critical differential growth γ_cr. We find that a thicker muscle substrate will cause an earlier instability. Remember that, since h_0 is constant in this case, a higher R_1 corresponds to a lower curvature. Therefore, we summarize by stating that curvature increase can delay surface wrinkles if the thickness of the skin layer is fixed. Furthermore, we discover that the inner radius R_0 has no influence on both the critical differential growth and critical wave number, thus we do not show the associated results here.





Figure <ref> illustrates the effect of the spring stiffness K on the critical state, namely that the critical differential growth γ_cr decreases as K increases, while the critical wave number n_cr is a two-step function of K.  





Finally, Figure <ref> shows the critical relative growth factor γ_cr and wave number n_cr when the relative modulus β varies. In this case, we set K=20, instead of K=1 adopted in previous examples, since this latter value is too small to generate wrinkles when β is large. From Figure <ref>, we have β_max≈4.4746 when K=1. However, if K=20, then the associated value becomes β_max≈73.5354. So wrinkles are only possible provided that β<β_max, as seen from Figure <ref>. 





We observe from Figure <ref> that the dependence of γ_cr on β is non-monotonic, which is different from the case when K→∞ <cit.>. By carefully analyzing the deformation, we find that the bifurcation is induced by the resultant axial force N. By equation (<ref>), N is proportional to both K and λ_z. When β increases, corresponding to a stiffer outer layer, a larger axial force N is required to generate an instability. Furthermore, Figure <ref> shows the critical stretch λ_z^cr for bifurcation as a function of β. As an increasing function is obtained, more force is necessary to trigger surface wrinkling for higher β. In Figure <ref>, the critical wave number n_cr reduces with increasing β, consistent with the other results.




§ ASYMPTOTIC ANALYSIS


In this section, we present approximate analytical results for the primary deformation and the bifurcation condition.




 §.§ Primary deformation


To make analytical progress, we set λ_z=1, corresponding to the case with Dirichlet boundary conditions at both ends of the cylindrical system. This is also the limiting case when the elastic spring is infinitely stiff, i.e., K→∞. Given that the inner surface of the tubular system is free, the equation for the elastic hoop stretch α_r_0 is

    g/α_r_0^2-g/α_r_1^2+βγ^2g/α_r_1^2-βγ^2g/α_h^2+logα_r_1^2/α_r_0^2+βγlogα_h^2/α_r_1^2=0,

where α_1^2 and α_h^2 depend on α_r_0^2 through the relations (<ref>). Assuming 

    α_r_0^2∼ g,    α_r_1^2∼ g,    α_h^2∼ g,

we have the Taylor expansion

    logα_r_1^2/α_r_0^2≈α_r_1^2/α_r_0^2-1,    logα_h^2/α_r_1^2≈α_h^2/α_r_1^2-1.

Then  (<ref>) is approximated as follows,

    g/α_r_0^2-g/α_r_1^2+βγ^2g/α_r_1^2-βγ^2g/α_h^2+α_r_1^2/α_r_0^2-1+βγ(α_h^2/α_r_1^2-1)=0,

or equivalently,

    α_h^2/g(α_r_1^2/g-α_r_0^2/g)(α_r_1^2/g+1)+βγα_r_0^2/g(α_h^2/g-α_r_1^2/g)(α_h^2/g+γ)=0.

Expressing α_1^2 and α_h^2 in terms of α_r_0^2, as in (<ref>), yields a cubic equation in x=α_r_0^2/g, which admits only one real solution and can be solved directly. In particular, if γ=1, then α_r_0=α_r_1=α_h=g^1/2.

When R_0∼ h_0^1/2, we can approximate α_r_0∼ g^1/2 by seeking a solution of the form

    α_r_0≈ x_0+x_1h_0+x_2h_0^2+x_3h_0^3+⋯,

where x_i, i=0,1,2,⋯, are unknowns to be determined. After inserting the above formula into the cubic equation for α_r_0^2/g, then expanding in powers of h_0 and collecting the coefficients of h_0^i, i=0,1,2,⋯, we obtain infinitely many algebraic equations. Therefore we are able to solve for x_0, x_1, x_2 and find an asymptotic solution for r_0=α_r_0 R_0 of the form

    r_0≈ g^1/2R_0[1+βγ h_02R_1(1+8 R_0^2R_1^2)(γ ^2-1)+3
       βγ(γ^2-1)^2-2(4γ^2-γ+3)(γ-1)8 R_1^2].

Note that r_0→ g^1/2R_0 as γ→ 1, which is the exact solution when there is no differential growth, i.e., the skin and the muscle grow simultaneously at the same rate. 





Figure <ref> illustrates the exact and asymptotic solutions for r_0 as functions of the differential growth factor γ, when the pre-growth factor g and muscle outer radius R_1 are given, the inner radius is R_0=0.1, the skin thickness is h_0=0.01, and the relative stiffness is β=1.27, as estimated in elephant trunk.




 §.§ WKB approximation


To derive analytical formulae for the critical stretch and critical wave number, we employ the WKB (Wentzel–Kramers–Brillouin) method <cit.>, which has been found useful in the context of wrinkling in curved and graded structures <cit.>. 

First, we see from Figure <ref> that the boundary spring stiffness K has a marginal influence on the critical wave number n_cr. Therefore, we can assume K→∞, which is equivalent to λ_z=1 (i.e., Dirichlet boundary conditions are imposed at both ends of the tubular system). 

Second, we consider γ=1 so there is no differential growth, so that the pre-growth factor g as the control parameter. As the main purpose of the asymptotic analysis is to gain analytical insight into the effect of different geometrical and material parameters on the critical state, this is a reasonable simplification.

Then the primary deformation can be analytically characterized, as follows. Starting from the original eigenvalue problem arising from equations (<ref>), we introduce a stream function 

    ψ(r,z)=f(r)sin(kz).

Then the incremental displacements are 

    u(r,z)=1/r∂ψ/∂ z,     w(r,z)=-1/r∂ψ/∂ r,

and the linearized incremental incompressibility condition is automatically satisfied. 

In terms of the stream function, the incremental equations and boundary conditions become, respectively,

    r^3𝒜_1313f””+2r^2(r𝒜_1313'-𝒜_1313)f”'+a_2(r)f”+a_1(r)f'+a_0f=0,

and

    {   ϱ_s11(r_1+h)=0,  ϱ_s31(r_1+h)=0,  ϱ_m11(r_0)=0,  ϱ_m31(r_0)=0,
       ϱ_s11(r_1)=ϱ_m11(r_1),  ϱ_s31(r_1)=ϱ_m31(r_1),   f_s(r_1)=f_m(r_1),   f^'_s(r_1)=f^'_m(r_1),
        .

where

    {   ϱ_11=r^2𝒜_1313f”'+(r^2𝒜_1313'-r𝒜_1313)f”+b_1(r)f'+b_0(r)f=0,
       ϱ_31=r𝒜_1313f”-𝒜_1313f'+k^2r(𝒜_1313+p)f=0,
    .

with the coefficients {a_i}_i=0,1,2 and {b_i}_i=0,1 taking the following forms: 

    {   a_0=k^2r(𝒜_2222-𝒜_1111+2𝒜_1133-2𝒜_2233+𝒜_1111'-𝒜_1122'-𝒜_1331'-𝒜_1133'+𝒜_2233')
       +k^4r^3𝒜_3131+k^2r^3𝒜_1331”+k^2r^3p”,
       
            a_1=k^2r^2(𝒜_1111-2𝒜_1331-2𝒜_1133+𝒜_3333)+k^2r^3(2𝒜_1133'-𝒜_1111'+2𝒜_1331'-𝒜_3333')
       -3𝒜_1313+3r𝒜_1313'-r^2𝒜_1313”,
       
            a_2=k^2r^3(2𝒜_1133-𝒜_1111+2𝒜_1331-𝒜_3333)+3r𝒜_1313-3r^2𝒜_1313'+r^3𝒜_1313',
       
            b_0=k^2r(𝒜_1111-𝒜_1122-𝒜_1133+𝒜_2233+p+𝒜_1331'+p'),
       
            b_1=k^2r(2𝒜_1133-𝒜_1111+𝒜_1331-𝒜_3333-p')+𝒜_1313-r𝒜_1313'.
        .


For the WKB-type solution to  (<ref>), we have 

    f=exp(∫_r_c^r𝒮(τ)dτ),

where 𝒮(r) is a function to be determined, the lower limit of integration is r_c=r_1 for the skin and r_c=r_0 for the muscle. 

To solve the eigenvalue problem with variable coefficients, arising from the mathematical model for transverse wrinkles when the wave number n=k/π is large, we look for solutions of the form

    𝒮(r)=n 𝒮_0(r)+𝒮_1(r)+1/n𝒮_2(r)+⋯,

where 𝒮_i (i=1,2,3,⋯) are unknown functions. 

Substituting the forms (<ref>) and (<ref>) into  (<ref>) then equating the coefficient of n to zero provides an algebraic equation for 𝒮_0. Solving this equation directly yields four independent solutions. Similarly, we are able to derive 𝒮_1 and 𝒮_2 in a systematic manner. 

For the neo-Hookean material and the restricted growth case considered here, namely, γ=1 and λ_z=1, we obtain the concise formulae given by

    𝒮_0^(1,2)=±π,  𝒮_0^(3,4)=±π g^-1/2,  𝒮_1^(1,2,3,4)=1/2r,  𝒮_2^(1,2)=±3/8π r^2,   𝒮_2^(3,4)=±3 g^1/2/8π r^2.


As a result, we express the general solutions as follows

    {   f_m(r)=∑_i=1^4 C_iexp(∫_r_0^r𝒮^(i)_m(τ)dτ),  r_0⩽ r ⩽ r_1,
       f_s(r)=∑_i=1^4 C_i+4exp(∫_r_1^r𝒮^(i)_s(τ)dτ),   r_1⩽ r ⩽ r_1+h.
    .

Substituting these into the boundary and interface conditions (<ref>), we obtain

    𝐌𝐂=0,

where 

    𝐂=[C_1,C_2,C_3,C_4,C_5,C_6,C_7,C_8]^T,

and

    𝐌=[
    [  M_11  M_12;  M_21 M_22, ]],

with the sub-matrices

    𝐌_11=[[         0         0         0         0;         0         0         0         0; 𝔄_m1(r_0) 𝔄_m2(r_0) 𝔄_m3(r_0) 𝔄_m4(r_0); 𝔅_m1(r_0) 𝔅_m2(r_0) 𝔅_m3(r_0) 𝔅_m4(r_0) ]],


    𝐌_12=[[ E_s1𝔄_s1(r_1+h) E_s2𝔄_s2(r_1+h) E_s3𝔄_s3(r_1+h) E_s4𝔄_s4(r_1+h); E_s1𝔅_s1(r_1+h) E_s2𝔅_s2(r_1+h)     𝔅_s3(r_1+h) E_s4𝔅_s4(r_1+h);               0               0               0               0;               0               0               0               0 ]],


    𝐌_21=[[     -E_m1𝔄_m1(r_1)     -E_m2𝔄_m2(r_1)     -E_m3𝔄_m3(r_1)     -E_m4𝔄_m4(r_1);     -E_m1𝔅_m1(r_1)     -E_m2𝔅_m2(r_1)     -E_m3𝔅_m3(r_1)     -E_m4𝔅_m4(r_1);              -E_m1              -E_m2              -E_m3              -E_m4; -E_m1𝒮_m^(1)'(r_1) -E_m2𝒮_m^(2)'(r_1) -E_m3𝒮_m^(3)'(r_1) -E_m4𝒮_m^(4)'(r_1) ]],


    𝐌_22=[[     𝔄_s1(r_1)     𝔄_s2(r_1)     𝔄_s3(r_1)     𝔄_s4(r_1);     𝔅_s1(r_1)     𝔅_s2(r_1)     𝔅_s3(r_1)     𝔅_s4(r_1);             1             1             1             1; 𝒮_s^(1)'(r_1) 𝒮_s^(2)'(r_1)  𝒮_s^(3'(r_1) 𝒮_s^(4)'(r_1);               ]],

and the components

    .
       E_mi=exp(∫_r_0^r_1𝒮_m^(i)dr),   E_si=exp(∫_r_1^r_1+h𝒮_s^(i)dr),
       𝔄_i=b_0+b_1𝒮^(i)'+r(r𝒜_1313'-𝒜_1313)(𝒮^(i)')^2+r^2𝒜_1313(𝒮^(i)')^3
       +r(r𝒜_1313'-𝒜_1313)𝒮_s^(i)''+3r^2𝒜_1313𝒮_s^(i)'𝒮_s^(i)''+r^2𝒜_1313𝒮_s^(i)''',
    
       𝔅_i=rk^2(𝒜_1313+p)-𝒜_1313𝒮^(i)'+r𝒜_1313(𝒮^(i)')^2+r𝒜_1313𝒮^(i)'',
    }  i=1,2,3,4.

Pursuing a non-trivial solution yields

    𝐌=0.


It can further to be deduced from (<ref>) that E_m1 and E_m3 are exponentially large. On the other hand, as the skin is very thin, i.e., h is also small, E_si, i=1,2,3,4, no longer possess the exponentially large or small nature. In view of these facts, the terms proportional to E_m1E_m3 are dominant such that 

    𝐌/E_m1E_m3=Ψ+.

We therefore obtain an approximate but explicit bifurcation condition:

    Ψ(g,n,h_0,β,R_0,R_1)=0.






Before proceeding further, we specify the parameter values to be used in the subsequent calculations. To this end, we plot in Figure <ref> the exact bifurcation curves, in the absence of relative growth (i.e., setting γ=1), when the radius R_1 is fixed, R_0=0.1, h_0=0.01, and β=1.27. We find that each curve has two local minima: a first one which increases as R_1 increases, with its critical wave number decreasing, and a second one which is independent of R_1, and thus remains the same for all the curves. For example, when R_1=0.25, the  first local minimum is the global minimum, attained at the wave number n=2, while when R_1=0.35, the second local minimum is the global minimum attained at n=17. This suggest a possible mode transition as the aspect ratio of the bilayer tube increases. Similar mode transitions have also been reported for compressed tubes <cit.> and film/substrate structures <cit.>. 

For the current problem, we can determine the critical value for the geometric parameter R_1 where the mode transition occurs by applying the following  two-step procedure: 

    
  (I) First, setting R_1=0.35 say, we determine g_min and n_min for the global minimum of the bifurcation curve by solving simultaneously the equations
    
    Φ(g,n)=0,    ∂Φ(g,n)/∂ n=0.

    
  (II) Second, substituting g=g_min in the bifurcation condition Φ, we identify R_1=R_1^(m) and the associated small wave number n^(m) by solving the following two equations  

    Φ(R_1,n)=0,    ∂Φ(R_1,n)∂ n=0.
 
    
By applying the above procedure when R_0=0.1, h_0=0.01, β=1.27 and γ=1, as in Figure <ref>, we obtain R_1^(m)≈ 0.2989.
    
This interesting mode transition deserved to be treated in detail in a separate study. As our attention here is focused on wrinkling, i.e., on the case when the wave number is large, we shall select R_1=0.4 in all subsequent examples and thus avoid the above mode transition. 





Figure <ref> compares the exact and asymptotic bifurcation curves for the pre-growth factor g as a function of the wave number n, without differential growth (i.e., we set γ=1), when the skin thickness h_0 is given, the muscle inner radius is R_0=0.1, its outer radius is R_1=0.25, and their modulus ratio is β=1.27 (small) or β=50 (large). Excellent agreement is found for all illustrated cases. 

Next, we concentrate on two distinct limiting cases, namely when β∼ 1 and when β→∞, respectively.




  §.§.§ The limiting case when the shear modulus ratio is close to unity
 

In elephant trunk, the skin and muscle substrate have comparable stiffness <cit.>, thus β∼ 1. In this case, it can be seen from Figure <ref> that the critical pre-growth factor g_cr is far from 1, even when the skin layer is very thin, i.e., 0<h_0≪ 1. In this case, we have two small parameters h_0 and 1/n, and approximate the critical growth factor and critical wave number as follows, <cit.>

    g_cr=g_0+g_1h_0+g_2h_0^2+⋯

and

    n_cr=h_0^-1(n_0+n_1h_0+n_2h_0^2+⋯).

Inserting the above forms into the simultaneous equations given by the bifurcation condition (<ref>) and the equation ∂Ψ/∂ n=0, we find that the two leading order equations for g_0 and n_0 can only be solved numerically. Once they are solved, all higher order unknowns are then derived recursively. Although we do not include here the lengthy expressions for the leading-order equations and the recursive relations, we point out that g_0 and n_0 are both independent of R_0 and R_1, and are only related to β, indicating that the critical pre-growth g_cr is dominated by β only. This is useful in explaining Figure <ref> where γ_cr belongs in the small interval [2.099,2.199]. Similarly, n_cr h_0 is governed by β, and thus the critical wave number n_cr decreases as h_0 increases, which explains the result in Figure <ref>. 

The three-term approximations (<ref>)-(<ref>) for the critical pre-growth factor g_cr and critical wave number n_cr when R_0=0.1, R_1=0.4, and ζ=0.0025 (h_0=0.01) are compared with the exact solutions in Figure <ref>. It can be seen from these plots that, for β close to 1, an excellent agreement is found, validating our asymptotic solutions. 








  §.§.§ The limiting case when the shear modulus ratio is large
 

When the modulus ratio β is large, we observe from Figure <ref> that the critical pre-growth factor is close to 1. So there are in total four small parameters, namely, g_cr-1, 1/β, 1/n, and h_0. Again, we assume that there is no differential growth, i.e., γ=1. To derive analytical formulae for the critical stretch and critical wave number, we set the scalings <cit.>

    g_cr-1∼𝒪(β^-2/3),     n_cr∼𝒪(β^2/3),     h_0∼𝒪(β^-1),

then pursue an asymptotic solution in terms of the small parameter 1/β. In a similar manner as in the previous subsection, we obtain the following explicit forms for the critical pre-growth factor 

    g_cr   =1+3^2/3β^-2/3+7/203^4/3β^-4/3+h_0^2/R_1^23^1/3β^2/3-h_0/2R_1β^-1
       +1/23^2/3β^-5/3+1009/2800 3β^-2
       +79 h_0^2/20 R_1^2-h_0^4/R_1^4β^2+𝒪(β^-7/3),

and the critical wave number

    n_cr   =1/π h_0[3^1/3β^-1/3-1/203^3β^-1+h_0^2/R_1^2β-h_0/2R_13^-1/3β^-2/3.
       .+47 h_0^2/40R_1^23^-1/3β^1/3-2 h_0^4/R_1^43^-1/3β^7/3+1229/29003^2/3β^-5/3+𝒪(β^-1)].







From equations (<ref>)-(<ref>) we infer that the inner radius R_0 is not involved up to the truncated order. Meanwhile, the muscle radius R_1 and the skin thickness h_0 affect the critical pre-growth factor trough higher order terms, indicating that their influence is relatively minor compared to that of the relative modulus β. In addition, both h_0 and β feature in the leading-order term in  (<ref>). Recalling the notation ζ=h_0/R_1, we can rewrite  

    g_cr   =1+3^2/3β^-2/3+7/203^4/3β^-4/3+ζ^2 3^1/3β^2/3-1/2ζβ^-1
       +1/23^2/3β^-5/3+1009/2800 3β^-2+79/20ζ^2-ζ^4β^2+𝒪(β^-7/3),

and 

    n_cr   =1/π h_0[3^1/3β^-1/3-1/203^3β^-1+ζ^2β-ζ/23^-1/3β^-2/3.
       .+47/40ζ^2 3^-1/3β^1/3-2ζ^4 3^-1/3β^7/3+1229/29003^2/3β^-5/3+𝒪(β^-1)].

It can be seen from the above equations that g_cr and n_crπ h_0 depend h_0 and R_1 through ζ. In particular, when β and ζ are fixed, the critical pre-growth factor g_cr is constant. 

Figure <ref> shows the exact and asymptotic solutions for the critical pre-growth factor g_cr given by equation (<ref>) and the critical wave number n_cr given by equation (<ref>), as functions of the modulus ratio β, when R_0=0.1, R_1=0.4, and ζ=0.04 (h_0=0.01), in the absence of differential growth (i.e., when γ=1).  

We also remark that, in equation (<ref>) for transverse wrinkles distributed along the tube length, the leading-order term 3^1/3β^-1/3 is the same as those for the circumferential wrinkles in growing bilayer tubes 
<cit.>, core/shell cylinders <cit.>, and surface wrinkles in planar film/substrate bilayers <cit.>. This analsyis further confirms the universal scaling between the modulus ratio and number of wrinkles.




§ CONCLUSION


This work stems from the need to understand wrinkling in elephant trunk where physical measurements suggest that skin and muscle substrate have comparable stiffness. Due to their mechanical compatibility, the two components can grow together seamlessly, and the wrinkled skin acts as a protective barrier that is at the same time thicker and more flexible than the unwrinkled skin. Moreover, geometric parameters, such as curvature, play key roles in the formation of transverse wrinkles. In particular, our model predicts that fewer wrinkles form and earlier in the proximal region where curvature is lower than in the distal region which has larger curvature, as observed in elephant trunks. Similarly, the dorsal side presents more wrinkling than the ventral side, since curvature is higher dorsally than ventrally.

While our theoretical and numerical results describe mathematically how transverse wrinkles form, our  investigation extends beyond that and can be useful to many other applications. Follow-up models should take into account additional deformations and loading conditions, e.g., bending, unbending, and inflation under internal pressure, as in a water filled elephant trunk, for example. More sophisticated hyperelastic strain-energy functions for skin and muscle tissues could also be considered. 

The elephant trunk is a rich source of inspiration for bio-mechanical devices, but its physiology is yet to be understood. We hope that our analysis will stimulate further quantitative studies of elephant trunk and elephant skin more generally. 





  
Acknowledgment
We gratefully acknowledge the UKRI Horizon Europe Guarantee MSCA (Marie Skłodowska-Curie Actions) Postdoctoral Fellowship to Yang Liu (EPSRC Grant No. EP/Y030559/1). Yang Liu further acknowledges the financial support from the National Natural Science Foundation of China (Project No. 12072227). 




  
Data availability statement
There are no additional data associated with this article.




  
Conflict of  interest

The authors have no competing interests to declare that are relevant to the content of this article.



9

[Ayer and Mariappa, 1950]Ayer:1950:AM Ayer AA, Mariappa D. 1950. External characters of three foetuses of the Indian elephant. Proceedings of the Indian Academy of Sciences B 31(4), 193–209 (doi: 10.1007/BF03050577).

[Ben Amar and Goriely, 2005]BenAmar2005JMPS 
    Ben Amar M., Goriely A., 2005. 
    Growth and instability in elastic tissues.
    Journal of the Mechanics and Physics of Solids 53, 2284–2319 (doi: 10.1016/j.jmps.2005.04.008).

[Biryukov, 1985]Biryukov1985 
    Biryukov, S.V., 1985. 
    Impedance method in the theory of elastic surface waves.
    Sov. Phys. Acoust. 31, 350–354.

[Biryukov et al., 1995]Biryukov1995 
    Birykov, S.V., Gulyaev, Y.V., Krylov, V.V., Plessky, V.P., 1995.
    Surface acoustic waves in inhomogeneous media.
    Springer Berlin, Heidelberg.

[Boas and Paulli, 1908]Boas:1908:BP Boas J.E.V., Paulli S., 1908. The Elephant's Head: Studies in the Comparative Anatomy of the Organs of the Head of the Indian Elephant and Other Mammals, published at the cost of the Carlsberg-Fund, Copenhagen.

[Dagenais et al., 2021]Dagenais:2021:etal Dagenais P., Hensman S., Haechler V., 2021. Elephants evolved strategies reducing the biomechanical complexity of their trunk. Current Biology 31(21), 4727–-4737 (doi: 10.1016/j.cub.2021.08.029).

[Dai and Liu, 2014]Dai2014EPL Dai, H.-H., Liu, Y., 2014. Critical thickness ratio for buckled and wrinkled fruits and vegetables. Europhysics Letters 108(4), 44003  (doi: 10.1209/0295-5075/108/44003).

[De Pascalis et al., 2010]DePascalis:2011:PDG De Pascalis R., Destrade M., Goriely A., 2010. Nonlinear correction to the Euler buckling formula for compressed cylinders with guided-guided end conditions. Journal of Elasticity 102 (2), 191–200 (doi: 10.1007/s10659-010-9265-6).

[Destrade et al., 2014]Destrade:2014:DOSV Destrade M., Ogden R.W., Sgura I., Vergori L., 2014. Straightening wrinkles, Journal of the Mechanics and Physics of Solids 65, 1–11 (doi: 10.1016/j.jmps.2014.01.001).

[Duka et al., 1993]Duka1993AM Duka E.D., England A.H., Spencer A.J.M., 1993. Bifurcation of a solid circular cylinder under finite extension and torsion. Acta Mechanica 98, 107–121 (doi: 10.1007/BF01174297).

[Eales, 1926]Eales:1926 Eales NB. 1926. XI—The Anatomy of the head of a foetal African elephant, Elephas africanus (Loxodonta africana). Transactions of the Royal Society of Edinburgh 54(3), 491–551. (doi: 10.1017/S0080456800016082).

[Emery and Fu, 2021]Emery:2021:EF Emery D., Fu Y., 2021. Elasto-capillary circumferential buckling of soft tubes under axial loading: existence and competition with localised beading and periodic axial modes. Mechanics of Soft Materials 3, 3 (doi: 10.1007/s42558-021-00034-x).

[Flügge, 1962]Flugge:1962 Flügge W., 1962. Statik und Dynamik der Schalen. Berlin, Springer. 

[Flügge, 1973]Flugge:1973 Flügge W., 1963. Stresses in Shells. Berlin, Heidelberg, New York, Springer-Verlag.

[Fosdick and Shield, 1963]Fosdick:1963:FS Fosdick R.A., Shield R.T., 1963. Small bending of a circular bar superposed on finite extension or compression. Archive for Rational Mechanics and Analysis 12, 223–248 (doi: 10.1007/BF00281227).

[Fu et al., 2016]Fu:2016:FLF Fu Y., Liu J.L., Francisco G.S., 2016. Localized bulging in an inflated cylindrical tube of arbitrary thickness - the effect of bending stiffness. Journal of the Mechanics and Physics of Solids 90, 45–60 (doi: 10.1016/j.jmps.2016.02.027).

[Fu et al., 2021]Fu:2021:FLG Fu Y., Jin L., Goriely A., 2021. Necking, beading, and bulging in soft elastic cylinders. Journal of the Mechanics and Physics of Solids 147, 104250 (doi: 10.1016/j.jmps.2020.104250).
    
[Goriely, 2017]goriely17 Goriely A., 2017.  The Mathematics and Mechanics of Biological Growth, Springer-Verlag, New York.

[Goriely et al., 2008]Goriely2008PRSA Goriely A., Vandiver R., Destrade M., 2008. Nonlinear Euler buckling. Proceedings of the Royal Society A 464(2099), 3003–3019 (doi:10.1098/rspa.2008.0184).

[Green et al., 1952]Green:1952:GRS Green A.E., Rivlin R.S., Shield R.T., 1952. General theory of small elastic deformations superposed on finite elastic deformations. Proceedings of the Royal Society A 211, 128–154 (doi: 10.1098/rspa.1952.0030).

[Gustafson and Abe, 1998]Gustafson:1998:GA Gustafson K., Abe T., 1998. The third boundary condition – Was it Robin's?. The Mathematical Intelligencer 20, 63–71 (doi: 10.1007/BF03024402).

[Hildebrandt et al., 2007]Hildebrandt:2007:etal Hildebrandt T., Drews B., Gaeth A.P., Goeritz F., Hermes R., Schmitt D., Gray C., Rich P., Streich W.J., Short R.V., Renfree M.B., 2007. Foetal age determination and development in elephants. Proceedings of the Royal Society B 274, 323–331 (doi: 10.1098/rspb.2006.3738).

[Hinch, 1991]Hinch1991 Hinch, E.J., 1991. Perturbation Methods. Cambridge University Press.

[Ilichev and Fu, 2014]Ilichev:2014:IF Il'ichev A.T., Fu Y., 2014. Stability of an inflated hyperelastic membrane tube with localized wall thinning. International Journal of Engineering Science 80, 53–61 (doi: 10.1016/j.ijengsci.2014.02.031).

[Jia et al., 2014]Jia:2014:etal Jia F., Cao Y.P., Zhao Y., Feng X.Q., 2014. Buckling and surface wrinkling of an elastic graded cylinder with elastic modulus arbitrarily varying along radial direction. International Journal of Applied Mechanics 6(1), 1450003 (doi: 10.1142/S1758825114500033).

[Jia et al., 2015]Jia2015PRE Jia F., Li B., Cao Y.P., Xie W.H., Feng X.Q., 2015. Wrinkling pattern evolution of cylindrical biological tissues with differential growth. Physical Review E 91(1), 012403 (doi: 10.1103/PhysRevE.91.012403).

[Jia et al., 2018]Jia:2018:etal Jia F., Pearce S.P., Goriely A., Curvature delays growth-induced wrinkling. Physical Review E 98, 033003 (doi: 10.1103/PhysRevE.98.033003).

[Jin et al., 2018]Jin:2018:etal Jin L., Liu Y., Cai Z., Asymptotic solutions on the circumferential wrinkling of growing tubular tissues. International Journal of Engineering Science  128, 31–43 (doi: 10.1016/j.ijengsci.2018.03.005).

[Kaczmarski et al., 2024a]Kaczmarski:2024a:etal Kaczmarski B., Leanza S., Zhao R., Kuhl E., Moulton D.E., Goriely A., 2024. Minimal design of the elephant trunk as an active filament. Physical Review Letters 132, 248402 (doi: 10.1103/PhysRevLett.132.248402).

[Kaczmarski et al., 2024b]Kaczmarski:2024b:etal Kaczmarski B., Moulton D.E., Goriely A., Kuhl E., 2024. Minimal activation with maximal reach: Reachability clouds of bio-inspired slender manipulators. Extreme Mechanics Letters Volume 71, 102207 (doi: 10.1016/j.eml.2024.102207).

[Kier and Smith, 1985]Kier:1985:KS Kier W.M., Smith K.K., 1985. Tongues, tentacles and trunks: The biomechanics of movement in muscular-hydrostats. Zoological Journal of the Linnean Society 83(4), 307–-324 (doi: 10.1111/j.1096-3642.1985.tb01178.x).

[Leanza et al., 2024]Leanza:2024:etal Leanza S., Lu-Yang J., Kaczmarski B., Wu S., Kuhl E., Zhao R.R., 2024. Elephant trunk inspired multimodal deformations and movements of soft robotic arms. Advanced Functional Materials, 202400396 (doi: 10.1002/adfm.202400396).

[Li and Yen, 1972]Li:1972:LY Li R.C.M., Yen K.H., 1972. Elastic waves guided by a solid layer between adjacent substrates. IEEE Transactions on Microwave Theory and Techniques 20(7), 477–486 (doi: 10.1109/TMTT.1972.1127788).

[Liu and Dai, 2014]Liu2014IJES Liu Y., Dai H.-H., 2014. Compression of a hyperelastic layer-substrate structure: Transitions between buckling and surface modes. International Journal of Engineering Science 80, 74–89 (doi: 10.1016/j.ijengsci.2014.02.020).

[Liu and Dorfmann, 2024]Liu2024MMS Liu Y., Dorfmann L., 2024. Localized necking and bulging of finitely deformed residually stressed solid cylinder. Mathematics and Mechanics of Solids 29(6), 1153–1175 (doi: 10.1177/10812865231186951).

[Liu et al., 2022]Liu:2022:etal Liu C., Du Y., Li K., Zhang Y., Han Z., Zhang Y., Qu S., Lü C., 2022. Geometrical incompatibility guides pattern selection in growing bilayer tubes. Journal of the Mechanics and Physics of Solids 169, 105087 (doi: 10.1016/j.jmps.2022.105087).

[Liu et al., 2024]Liu:2024:etal Liu R.-C., Liu Y., Goriely A., 2024. Surface wrinkling of a film coated to a graded substrate. Journal of the Mechanics and Physics of Solids, 186,  105603 (doi: 10.1016/j.jmps.2024.105603).

[Pearce and Fu, 2010]Pearce:2010:PF Pearce S.P., Fu Y., 2010. Characterization and stability of localized bulging/necking in inflated membrane tubes. IMA Journal of Applied Mathematics 75(4), 581–602 (doi: 10.1093/imamat/hxq026).

[Schulz et al., 2020]Schulz:2020:SFSSH Schulz A.K., Fourney E., Sordilla S., Sukhwani A., Hu D.L., 2020. Elephant trunk skin: Nature's flexible kevlar. International Conference on Intelligent Robots and Systems (IROS), 10.

[Schulz et al., 2022a]Schulz:2022a:SBBSRHARHH Schulz A.K., Boyle M., Boyle C., Sordilla S., Rincon C., Hooper S., Aubuchon C., Reidenberg J.S., Higgins. C, Hu D.L., 2022. Skin wrinkles and folds enable asymmetric stretch in the elephant trunk. PNAS Biophysics and Computational Biology 119(31), e2122563119 (doi: 10.1073/pnas.2122563119).

[Schulz et al., 2022b]Schulz:2022b:SRWTSMEH Schulz A.K., Reidenberg J.S., Wu J.N., Tang C.Y., Seleb B., Mancebo J., Elgart N., Hu D.L., 2022. Elephants trunks use an adaptable prehensile grip. Bioinspiration & Biomimetics 18, 026228 (doi: 10.1088/1748-3190/acb477).

[Schulz et al., 2023]Schulz:2023:SSZS Schulz A.K., Schneider N., Zhang M., Singal K., 2023. A year at the forefront of hydrostatic motion. Biology Open 12, bio059834 (doi: 10.1242/bio.059834).

[Schulz et al., 2024]Schulz:2024:SRKRHB Schulz A.K., Reveyaz N., Kaufmann L., Ritter C., Hildebrandt T., Brecht M., 2023. Elephants develop wrinkles through both form & function. Society of Integrative and Comparative Biology, Seattle, USA, 1 (doi: 10.1101/2023.08.24.554618).

[Shuvalov, 2003a]Shuvalov2003PRSA 
    Shuvalov, A.L., 2003.
    A sextic formalism for three-dimensional elastodynamics of cylindrically anisotropic radially inhomogeneous materials.
    Proc. R. Soc. Lond. A 459, 1611–1639.
    
    [Shuvalov, 2003b]Shuvalov2003QJMAM 
    Shuvalov, A.L., 2003.
    The frobenius power series solution for cylindrically anisotropic radially inhomogeneous elastic materials.
    Q. J. Mech. Appl. Math, 56(3), 327–345.

[Sigaeva et al., 2018]Sigaeva:2018:etal Sigaeva T., Mangan R., Vergori L., Destrade M., Sudak L., 2018. Wrinkles and creases in the bending, unbending and eversion of soft sectors. Proceedings of the Royal Society A 474, 20170827 (doi: 0.1098/rspa.2017.0827). 

[Springhetti et al., 2023]Springhetti:2023:SRB Springhetti R., Rossetto G., Bigoni D., 2023. Buckling of thin-walled cylinders from three dimensional nonlinear elasticity. Journal of Elasticity 154(1-4), 297–323 (doi: 10.1007/s10659-022-09905-4).

[Stroh, 1962]Stroh:1962 Stroh A.N., 1962. Steady state problems in anisotropic elasticity. Journal of Mathematics and Physics 41, 77–103 (doi: 10.1002/sapm196241177).

[Thompson, 1942]Thompson:1942 Thompson D.W., 1942. On Growth and Form, Cambridge University Press, Cambridge.

[Trivedi et al., 2008] Trivedi:2008:etal Trivedi D., Rahn C.D., Kier W.M., Walker I.D., 2008. Soft robotics: biological inspiration, state of the art, and future research. Applied Bionics and Biomechanics 5(3), 99–117 (doi: 10.1080/11762320802557865).

[Wang and Fu, 2021]Wang:2021:WF Wang M., Fu Y., 2021. Necking of a hyperelastic solid cylinder under axial stretching: Evaluation of the infinite-length approximation. International Journal of Engineering Science 159, 103432 (doi: 10.1016/j.ijengsci.2020.103432).

[Wang et al., 2019]Wang:2019:etal Wang S.B., Guo G.M., Zhou L., Li L.A., Fu Y., 2019. An experimental study of localized bulging in inflated cylindrical tubes guided by newly emerged analytical results. Journal of the Mechanics and Physics of Solids 124, 536–554 (doi: 10.1016/j.jmps.2018.11.011).

[Wilkes, 1955]Wilkes1955 Wilkes E.W., 1955. On the stability of a circular tube under end thrust. The Quarterly Journal of Mechanics and Applied Mathematics 8, 88–100 (doi:10.1093/qjmam/8.1.88).

[Wilson et al., 1991]Wilson:1991:etal Wilson JF, Mahajan U, Wainwright SA, Croner LJ. 1991. A continuum model of elephant trunks. ASME Journal of Biomechanical Engineering 113(1), 79–84 (doi: 10.1115/1.2894088).

[Woo and Shield, 1962]Woo:1962:WS Woo T.C., Shield R.T., 1962. Fundamental solutions for small deformations superposed on finite biaxial extension of an elastic body. Archive for Rational Mechanics and Analysis 9, 196–224 (doi: 10.1007/BF00253345).

[Wu et al., 2024]Wu:2024:etal Wu W., Yin Y., Li Y., Fan X., 2024. Theoretical analysis of inflated tube wrinkling behavior under pure bending. International Journal of Mechanical Sciences 273, 109166 (doi: 10.1016/j.ijmecsci.2024.109166).

[Ye et al., 2020]Ye:2020:YLF Ye Y., Liu Y., Fu Y., 2020. Weakly nonlinear analysis of localized bulging of an inflated hyperelastic tube of arbitrary wall thickness. Journal of the Mechanics and Physics of Solids 135, 103804 (doi: 10.1016/j.jmps.2019.103804).

[Ye et al., 2019]Ye:2019:etal Ye S., Yin S.F., Li B., Feng X.Q., 2019. Torsion instability of anisotropic cylindrical tissues with growth. Acta Mechanica Solida Sinica 32(5), 621–632 (doi: 10.1007/s10338-019-00087-6).

[Zhang et al., 2023]Zhang:2023:etal Zhang J., Li Y., Kan Z., Yuan Q., Rajabi H., Wu Z., Peng H., Wu J., 2023. A preprogrammable continuum robot inspired by elephant trunk for dexterous manipulation. Soft Robotics 10(3), 636-646 (doi: 10.1089/soro.2022.0048).

[Zhao et al., 2014]Zhao2014JMPS Zhao Y., Cao Y.P., Feng X.Q., Ma K., 2014. Axial compression-induced wrinkles on a core–shell soft cylinder: Theoretical analysis, simulations and experiments. Journal of the Mechanics and Physics of Solids 73, 212–227 (doi: 10.1016/j.jmps.2014.09.005). 
 
[Zhu et al., 2008]Zhu:2008:ZLO Zhu Y., Luo X.Y., Ogden R.W., 2008. Asymmetric bifurcations of thick-walled circular cylindrical elastic tubes under axial loading and external pressure. International Journal of Solids and Structures 45, 3410–3429 (doi: 10.1016/j.ijsolstr.2008.02.005).
 
 
 


